diff --git a/.gitignore b/.gitignore index 2d8f1da0f9..cc3be5b514 100644 --- a/.gitignore +++ b/.gitignore @@ -75,6 +75,8 @@ tutorials/13-joint_inversion/cross_gradient_data/* # jupyter *.ipynb +*.ipynb_checkpoints/ +*.idea/ # results and testing files *.obs diff --git a/SimPEG/electromagnetics/time_domain/simulation_1d.py b/SimPEG/electromagnetics/time_domain/simulation_1d.py index 72f9f07d77..2102cb20bd 100644 --- a/SimPEG/electromagnetics/time_domain/simulation_1d.py +++ b/SimPEG/electromagnetics/time_domain/simulation_1d.py @@ -34,7 +34,7 @@ try: from pyMKL import mkl_set_num_threads except Exception as e: - print("No MKL support: ", e) + # print("No MKL support: ", e) mkl_set_num_threads = None diff --git a/SimPEG/electromagnetics/utils/static_instrument/base.py b/SimPEG/electromagnetics/utils/static_instrument/base.py index 4d07a56b61..a9e033dec6 100644 --- a/SimPEG/electromagnetics/utils/static_instrument/base.py +++ b/SimPEG/electromagnetics/utils/static_instrument/base.py @@ -2,6 +2,7 @@ import os from matplotlib import pyplot as plt from discretize import TensorMesh +import multiprocessing from SimPEG import maps from SimPEG.electromagnetics import time_domain as tdem @@ -149,8 +150,9 @@ def times(self): for times_full, times_filter in zip(self.times_full, self.times_filter)] - startmodel__n_layer = 30 + startmodel__n_layer : int = 30 "Number of layers in model discretization" + @property def n_layer_used(self): if "resistivity" in self.xyz.layer_data: @@ -176,15 +178,26 @@ def data_uncert_array_culled(self): return np.where(np.isnan(dobs), np.inf, self.data_uncert_array) dipole_moments = [1] - - uncertainties__std_data = 0.03 - "Noise as a factor of data" - uncertainties__std_data_override = False - "If set to true, use the std_data value instead of data std:s from stacking" - uncertainties__noise_level_1ms = 1e-9 - "Amplitude at 1ms, in V/m^2" - uncertainties__noise_exponent = -0.5 - "Slope of noise floor, t^uncertainties__noise_exponent" + + # uncertainties__std_data = 0.03 + # "Noise as a factor of data" + @property + def uncertainties__std_data(self): + channel = 1 + return self.gex.gex_dict[f'Channel{channel}']['UniformDataSTD'] + + # uncertainties__std_data_override: bool = False + # "If set to true, use the std_data value instead of data std:s from stacking" + @property + def uncertainties__std_data_override(self): + channel = 1 + return f"dbdt_std_ch{channel}gt" not in self.xyz.layer_data.keys() + + uncertainties__noise_level_1ms : float = 1e-9 + "Noise amplitude at 1e-3 seconds. Unit of value is V/m^2" + uncertainties__noise_exponent : float = -0.5 + "Slope of noise-floor. Typically expressed as: t^uncertainties__noise_exponent" + @property def uncert_array(self): n_sounding = len(self.xyz.flightlines) @@ -205,13 +218,13 @@ def uncert_array(self): return np.where(np.isnan(self.data_array_nan), np.Inf, uncertainties) - startmodel__thicknesses_type: typing.Literal['logspaced', 'geometric', 'time'] = "logspaced" + startmodel__thicknesses_type : typing.Literal['logspaced', 'geometric', 'time'] = "logspaced" "Type of model discretization" - startmodel__thicknesses_minimum_dz = 1 + startmodel__thicknesses_minimum_dz : float = 1 "Thickness of thinnest layer if using 'logspaced' or 'geometric' discretization" - startmodel__top_depth_last_layer = 400 - "Depth to the top of the last layer if using logspaced discretization" - startmodel__thicknesses_geomtric_factor = 1.15309 + startmodel__top_depth_last_layer : float = 400 + "Depth to the top of the last layer if using log-spaced discretization" + startmodel__thicknesses_geomtric_factor : float = 1.15309 "Ratio of one layer to the next if using geometric discretization" def make_thicknesses(self): @@ -255,13 +268,24 @@ def make_survey(self): def n_param(self, thicknesses): return (len(thicknesses)+1)*len(self.xyz.flightlines) - + simulation__solver : typing.Literal['LU', 'pardiso'] = 'LU' "Equation solver to use. Some solvers might require hardware support." - simulation__parallel = True - "Use multiple computation threads in parallel. Useful to set to false in a notebook for debugging." - simulation__n_cpu = 3 - "Number of threads (roughly same as CPU cores) to use" + + @property + def simulation__parallel(self): + return True + # simulation__parallel: bool = True + # "Use multiple computation threads in parallel. Useful to set to false in a notebook for debugging." + + # FIXME!!! the simulation__n_cpu will need adjusted when we start to spin up dedicated machines for inversion + # like this property will basically be "choose your inversion machine" + @property + def simulation__n_cpu(self): + return multiprocessing.cpu_count() + # simulation__n_cpu = 3 + # "Number of threads (roughly same as CPU cores) to use" + def make_simulation(self, survey, thicknesses): if 'pardiso' in self.simulation__solver.lower(): print('Using Pardiso solver') @@ -302,7 +326,7 @@ def make_misfit(self, thicknesses): dmis.W = self.make_misfit_weights() return dmis - startmodel__res=100. + startmodel__res : float = 100. "Initial resistivity (ohmm)" def make_startmodel(self, thicknesses): startmodel=np.log(np.ones(self.n_param(thicknesses)) * 1/self.startmodel__res) @@ -321,9 +345,9 @@ def make_startmodel(self, thicknesses): # 5) line spacing (if 100m then alpha_s = 1e-4, if 400m then 6.3e-6) # 6) geomean of the linespacing and sounding spacing: sqrt(line_space * sound_space) # - 25m sounding spacing, 100m line spacing: h=50, alpha_s=4e-4 - regularization__alpha_s = 1e-4 - regularization__alpha_r = 1. - regularization__alpha_z = 1. + regularization__alpha_s : float = 1e-4 + regularization__alpha_r : float = 1. + regularization__alpha_z : float = 1. def make_regularization(self, thicknesses): if False: @@ -367,18 +391,19 @@ def make_regularization(self, thicknesses): "Random seed for beta (regularization) schedule estimator" directives__beta__beta0_ratio : float = 10. "Start ratio for the beta (regularization) schedule estimator" - directives__beta__cooling_factor=2 + directives__beta__cooling_factor : float = 2 "Cooling factor for the beta (regularization) schedule" - directives__beta__cooling_rate=1 + directives__beta__cooling_rate : float = 1 "Initial cooling rate for the beta (regularization) schedule" - directives__irls__enable = False - "IRLS is used to generate a sparse model in addition to and l2 model" - directives__irls__max_iterations = 30 + directives__irls__enable : bool = False + "IRLS is used to generate a sparse [sharp] model in addition to the l2 [smooth] model" + directives__irls__max_iterations : int = 30 "Maximum number of iterations (after l2 model has converged)" - directives__irls__minGNiter = 1 - directives__irls__fix_Jmatrix = True - directives__irls__f_min_change = 1e-3 - directives__irls__coolingRate = 1 + # FIXME!!! fill out descriptions for the rest of the irls directives + directives__irls__minGNiter : int = 1 + directives__irls__fix_Jmatrix : bool = True + directives__irls__f_min_change : float = 1e-3 + directives__irls__coolingRate : float = 1 def make_directives(self): if self.directives__beta__seed: BetaEstimate = directives.BetaEstimate_ByEig(beta0_ratio=self.directives__beta__beta0_ratio, @@ -405,9 +430,11 @@ def make_directives(self): return dirs - optimizer__max_iter=40 - "Maximum number of gauss newton iterations" - optimizer__max_iter_cg=20 + optimizer__max_iter : int = 50 + "Maximum number of inversion iterations" + optimizer__max_iter_cg : int = 20 + "Maximum number of Inexact-Gauss-Newton iterations" + def make_optimizer(self): return optimization.InexactGaussNewton(maxIter = self.optimizer__max_iter, maxIterCG=self.optimizer__max_iter_cg) diff --git a/SimPEG/electromagnetics/utils/static_instrument/dual.py b/SimPEG/electromagnetics/utils/static_instrument/dual.py index 0f552a3f76..c8e9befb0f 100644 --- a/SimPEG/electromagnetics/utils/static_instrument/dual.py +++ b/SimPEG/electromagnetics/utils/static_instrument/dual.py @@ -1,23 +1,18 @@ -import numpy as np import os +import tarfile +import typing + +import matplotlib as mpl from matplotlib import pyplot as plt -from discretize import TensorMesh +from matplotlib.colors import LogNorm -from SimPEG import maps -from SimPEG.electromagnetics import time_domain as tdem -from SimPEG.electromagnetics.utils.em1d_utils import plot_layer import libaarhusxyz import pandas as pd - import numpy as np +import scipy.stats from scipy.spatial import cKDTree, Delaunay -import os, tarfile -import matplotlib as mpl -from matplotlib import pyplot as plt -from matplotlib.colors import LogNorm from discretize import TensorMesh, SimplexMesh -from SimPEG.utils import mkvc from SimPEG import ( maps, data, data_misfit, inverse_problem, regularization, optimization, directives, inversion, utils @@ -26,17 +21,18 @@ from SimPEG.utils import mkvc import SimPEG.electromagnetics.time_domain as tdem import SimPEG.electromagnetics.utils.em1d_utils -from SimPEG.electromagnetics.utils.em1d_utils import get_2d_mesh,plot_layer, get_vertical_discretization_time +from SimPEG.electromagnetics.utils.em1d_utils import ( + get_2d_mesh, plot_layer, get_vertical_discretization_time + ) from SimPEG.regularization import LaterallyConstrained, RegularizationMesh -import scipy.stats from . import base -import typing + class DualMomentTEMXYZSystem(base.XYZSystem): """Dual moment system, suitable for describing e.g. the SkyTEM instruments. This class can not be directly instantiated, but - instead, instantiable subclasses can created using the class + instead, instantiable subclasses can be created using the class method ``` @@ -49,19 +45,42 @@ class DualMomentTEMXYZSystem(base.XYZSystem): See the help for `XYZSystem` for more information on basic usage. """ - gate_filter__start_lm=5 - "Lowest used gate (zero based)" - gate_filter__end_lm=11 - "First unused gate above used ones (zero based)" - gate_filter__start_hm=12 - "Lowest used gate (zero based)" - gate_filter__end_hm=26 - "First unused gate above used ones (zero based)" - - rx_orientation : typing.Literal['x', 'y', 'z'] = 'z' - "Receiver orientation" - tx_orientation : typing.Literal['x', 'y', 'z'] = 'z' - "Transmitter orientation" + + def gt_filt_st(self, inuse_ch_key): + return np.where((self.xyz.layer_data[inuse_ch_key].sum() == 0).cumsum().diff() == 0)[0][0] + + def gt_filt_end(self, inuse_ch_key): + return len(self.xyz.layer_data[inuse_ch_key].loc[1, :]) - np.where((self.xyz.layer_data[inuse_ch_key].sum().loc[::-1] == 0).cumsum().diff() == 0)[0][0] + + @property + def gate_filter__start(self): + start_list = [] + for channel in range(1, 1 + self.gex.number_channels): + # str_ch = f"0{channel}"[-2:] + # inuse_ch_key = f"InUse_Ch{str_ch}" + inuse_ch_key = f"dbdt_inuse_ch{channel}gt" + start_list.append(self.gt_filt_st(inuse_ch_key)) + return start_list + + @property + def gate_filter__end(self): + end_list = [] + for channel in range(1, 1 + self.gex.number_channels): + # str_ch = f"0{channel}"[-2:] + # inuse_ch_key = f"InUse_Ch{str_ch}" + inuse_ch_key = f"dbdt_inuse_ch{channel}gt" + end_list.append(self.gt_filt_end(inuse_ch_key)) + return end_list + + # FIXME!!! Tx_orientation should also be broken up by moment/channel, + # but without seeing how this would look in a gex I don't know how to redo this + @property + def tx_orientation(self): + return self.gex.tx_orientation + + @property + def rx_orientation(self): + return [(self.gex.rx_orientation(channel)).lower() for channel in range(1, 1 + self.gex.number_channels)] @classmethod def load_gex(cls, gex): @@ -77,19 +96,24 @@ class GexSystem(cls): @property def sounding_filter(self): # Exclude soundings with no usable gates - return self._xyz.dbdt_inuse_ch1gt.values.sum(axis=1) + self._xyz.dbdt_inuse_ch2gt.sum(axis=1) > 0 + # return self._xyz.dbdt_inuse_ch1gt.values.sum(axis=1) + self._xyz.dbdt_inuse_ch2gt.sum(axis=1) > 0 + num_gates_per_sounding_per_moment = {} + for channel in range(self.gex.number_channels): + iu_key = f"dbdt_inuse_ch{channel + 1}gt" + num_gates_per_sounding_per_moment[channel] = self._xyz.layer_data[iu_key].values.sum(axis=1) + num_gates_per_sounding = pd.DataFrame.from_dict(num_gates_per_sounding_per_moment) + print(f"num_gates_per_sounding = {num_gates_per_sounding}") + print(f"num_gates_per_sounding.sum(axis=1) > 0 = {num_gates_per_sounding.sum(axis=1) > 0}") + return num_gates_per_sounding.sum(axis=1) > 0 @property def area(self): return self.gex.General['TxLoopArea'] - - @property - def waveform_hm(self): - return self.gex.General['WaveformHMPoint'] - + @property - def waveform_lm(self): - return self.gex.General['WaveformLMPoint'] + def waveform(self): + return [self.gex.General[f"Waveform{self.gex.transmitter_moment(channel)}Point"] for channel in range(1, 1 + self.gex.number_channels)] + @property def correct_tilt_pitch_for1Dinv(self): @@ -104,7 +128,8 @@ def lm_data(self): dbdt = np.where(self.xyz.dbdt_inuse_ch1gt == 0, np.nan, dbdt) tiltcorrection = self.correct_tilt_pitch_for1Dinv tiltcorrection = np.tile(tiltcorrection, (dbdt.shape[1], 1)).T - return - dbdt * self.xyz.model_info.get("scalefactor", 1) * self.gex.Channel1['GateFactor'] * tiltcorrection + channel = 1 + return - dbdt * self.xyz.model_info.get("scalefactor", 1) * self.gex[f"Channel{channel}"]['GateFactor'] * tiltcorrection @property def hm_data(self): @@ -113,7 +138,8 @@ def hm_data(self): dbdt = np.where(self.xyz.dbdt_inuse_ch2gt == 0, np.nan, dbdt) tiltcorrection = self.correct_tilt_pitch_for1Dinv tiltcorrection = np.tile(tiltcorrection, (dbdt.shape[1], 1)).T - return - dbdt * self.xyz.model_info.get("scalefactor", 1) * self.gex.Channel2['GateFactor'] * tiltcorrection + channel = 2 + return - dbdt * self.xyz.model_info.get("scalefactor", 1) * self.gex[f"Channel{channel}"]['GateFactor'] * tiltcorrection # NOTE: dbdt_std is a fraction, not an actual standard deviation size! @property @@ -134,31 +160,43 @@ def data_uncert_array(self): @property def dipole_moments(self): - return [self.gex.gex_dict['Channel1']['ApproxDipoleMoment'], - self.gex.gex_dict['Channel2']['ApproxDipoleMoment']] + return [self.gex.gex_dict[f'Channel{channel}']['ApproxDipoleMoment'] for channel in range(1, 1 + self.gex.number_channels)] + # return [self.gex.gex_dict['Channel1']['ApproxDipoleMoment'], + # self.gex.gex_dict['Channel2']['ApproxDipoleMoment']] @property def times_full(self): - return (np.array(self.gex.gate_times('Channel1')[:,0]), - np.array(self.gex.gate_times('Channel2')[:,0])) + return tuple(np.array(self.gex.gate_times(channel)[:, 0]) for channel in range(1, 1 + self.gex.number_channels)) + # return (np.array(self.gex.gate_times(1)[:,0]), + # np.array(self.gex.gate_times(2)[:,0])) @property def times_filter(self): times = self.times_full filts = [np.zeros(len(t), dtype=bool) for t in times] - filts[0][self.gate_filter__start_lm:self.gate_filter__end_lm] = True - filts[1][self.gate_filter__start_hm:self.gate_filter__end_hm] = True + # FIXME: loop over channel number + for channel in range(self.gex.number_channels): + filts[channel][self.gate_filter__start[channel]: self.gate_filter__end[channel]] = True return filts - + + @property + def uncertainties__std_data(self): + return [self.gex.gex_dict[f'Channel{channel}']['UniformDataSTD'] for channel in range(1, 1 + self.gex.number_channels)] + + @property + def uncertainties__std_data_override(self): + # "If set to true, use the std_data value instead of data std:s from stacking" + return [f"dbdt_std_ch{channel}gt" not in self.xyz.layer_data.keys() for channel in range(1, 1 + self.gex.number_channels)] + def make_waveforms(self): - time_input_currents_hm = self.waveform_hm[:,0] - input_currents_hm = self.waveform_hm[:,1] - time_input_currents_lm = self.waveform_lm[:,0] - input_currents_lm = self.waveform_lm[:,1] - - waveform_hm = tdem.sources.PiecewiseLinearWaveform(time_input_currents_hm, input_currents_hm) - waveform_lm = tdem.sources.PiecewiseLinearWaveform(time_input_currents_lm, input_currents_lm) - return waveform_lm, waveform_hm + time_input_currents = [] + input_currents = [] + for channel in range(self.gex.number_channels): + time_input_currents.append(self.waveform[channel][:,0]) + input_currents.append(self.waveform[channel][:,1]) + + return [tdem.sources.PiecewiseLinearWaveform(time_input_currents[channel], input_currents[channel]) for channel in range(self.gex.number_channels)] + def make_system(self, idx, location, times): # FIXME: Martin says set z to altitude, not z (subtract topo), original code from seogi doesn't work! @@ -166,22 +204,12 @@ def make_system(self, idx, location, times): receiver_location = (location[0] + self.gex.General['RxCoilPosition'][0], location[1], location[2] + np.abs(self.gex.General['RxCoilPosition'][2])) - waveform_lm, waveform_hm = self.make_waveforms() - - return [ - tdem.sources.MagDipole( - [tdem.receivers.PointMagneticFluxTimeDerivative( - receiver_location, times[0], self.rx_orientation)], - location=location, - waveform=waveform_lm, - orientation=self.tx_orientation, - i_sounding=idx), - tdem.sources.MagDipole( - [tdem.receivers.PointMagneticFluxTimeDerivative( - receiver_location, times[1], self.rx_orientation)], - location=location, - waveform=waveform_hm, - orientation=self.tx_orientation, - i_sounding=idx)] - - + waveforms = self.make_waveforms() + return [tdem.sources.MagDipole([tdem.receivers.PointMagneticFluxTimeDerivative(receiver_location, + times[channel], + self.rx_orientation[channel])], + location=location, + waveform=waveforms[channel], + orientation=self.tx_orientation, + i_sounding=idx) + for channel in range(self.gex.number_channels)] diff --git a/SimPEG/electromagnetics/utils/static_instrument/single.py b/SimPEG/electromagnetics/utils/static_instrument/single.py index 8cb4316c1c..51f03629b5 100644 --- a/SimPEG/electromagnetics/utils/static_instrument/single.py +++ b/SimPEG/electromagnetics/utils/static_instrument/single.py @@ -24,7 +24,7 @@ class SingleMomentTEMXYZSystem(base.XYZSystem): @property def dipole_moments(self): - return [self.aera * self.i_max] + return [self.area * self.i_max] def make_waveforms(self): if self.waveform is None: diff --git a/SimPEG/simulation.py b/SimPEG/simulation.py index 25e3f291a7..27f7b78097 100644 --- a/SimPEG/simulation.py +++ b/SimPEG/simulation.py @@ -21,7 +21,7 @@ print("Using Pardiso solver as default if not otherwise defined") except ImportError: from .utils.solver_utils import SolverLU as DefaultSolver - print("Cannot import Pardiso, will use SolverLU") + # print("Cannot import Pardiso, will use SolverLU") __all__ = ["LinearSimulation", "ExponentialSinusoidSimulation"]