simpegEM1D
simpegEM1D copied to clipboard
Dual moment fails at interpolation
Ill attach two examples causing this to fail.
The first is the example in simpegEM1D but with the order of the HM and LM information switched during instantiation of the TDsurvey. The second is a custom walkTEM system, I simply specified the time gate center times and waveforms using numpy arrays.
Each causes an out of bounds error in Scipy's interpolate.
Example 1
Notebook cell
from SimPEG import Mesh, Maps, Utils
%matplotlib notebook
%load_ext autoreload
%autoreload 2
import matplotlib.pyplot as plt
from simpegEM1D import (
EM1D, EM1DSurveyTD, EM1DAnalytics, piecewise_pulse, set_mesh_1d,
skytem_HM_2015, skytem_LM_2015, get_vertical_discretization_time
)
import numpy as np
from scipy import io
from scipy.interpolate import interp1d
wave_HM = skytem_HM_2015()
wave_LM = skytem_LM_2015()
time_HM = wave_HM.time_gate_center[0::2]
time_LM = wave_LM.time_gate_center[0::2]
hz = get_vertical_discretization_time(
np.unique(np.r_[time_LM, time_HM]), facter_tmax=0.5, factor_tmin=10.,
n_layer=2
)
mesh1D = set_mesh_1d(hz)
depth = -mesh1D.gridN[:-1]
LocSigZ = -mesh1D.gridCC
time_input_currents_HM = wave_HM.current_times[-7:]
input_currents_HM = wave_HM.currents[-7:]
time_input_currents_LM = wave_LM.current_times[-13:]
input_currents_LM = wave_LM.currents[-13:]
TDsurvey = EM1DSurveyTD(
rx_location=np.array([0., 0., 0.]),
src_location=np.array([0., 0., 0.]),
topo=np.r_[0., 0., 0.],
depth=depth,
rx_type='dBzdt',
wave_type='general',
src_type='CircularLoop',
a=13.,
I=1.,
time=time_LM,
time_input_currents=time_input_currents_LM,
input_currents=input_currents_LM,
n_pulse=2,
base_frequency=210.,
use_lowpass_filter=True,
high_cut_frequency=7e4,
moment_type='dual',
time_dual_moment=time_HM,
time_input_currents_dual_moment=time_input_currents_HM,
input_currents_dual_moment=input_currents_HM,
base_frequency_dual_moment=25,
# half_switch = True
)
sig_half=1e-2
chi_half=0.
expmap = Maps.ExpMap(mesh1D)
m_1D = np.log((np.arange(TDsurvey.n_layer)+1.0)*sig_half)
chi = np.zeros(TDsurvey.n_layer)
prob = EM1D(
mesh1D, sigmaMap=expmap, chi=chi
)
prob.pair(TDsurvey)
dBzdtTD = TDsurvey.dpred(m_1D)
Error message
Produces the following with not much to go by.
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-7-6d55b4039c6a> in <module>()
1 prob.unpair()
2 prob.pair(TDsurvey)
----> 3 dBzdtTD = TDsurvey.dpred(m_1D)
4 err = np.linalg.norm(p_saved-dBzdtTD)
5 print(err, err < 1e-12)
/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Utils/codeutils.py in requiresVarWrapper(self, *args, **kwargs)
214 if getattr(self, var, None) is None:
215 raise Exception(extra)
--> 216 return f(self, *args, **kwargs)
217
218 doc = requiresVarWrapper.__doc__
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/Survey.py in dpred(self, m, f)
535 """
536 if f is None:
--> 537 f = self.prob.fields(m)
538
539 return self._pred
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/EM1D.py in fields(self, m)
423 def fields(self, m):
424 f = self.forward(m, output_type='response')
--> 425 self.survey._pred = Utils.mkvc(self.survey.projectFields(f))
426 return f
427
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/Survey.py in projectFields(self, u)
471 self.input_currents_dual_moment,
472 self.period_dual_moment,
--> 473 n_pulse=self.n_pulse
474 )
475 # concatenate dual moment response
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/Waveforms.py in piecewise_pulse(step_func, t_off, t_currents, currents, T, n, n_pulse)
249 response = (
250 piecewise_ramp(
--> 251 step_func, t_off, t_currents, currents, n=n
252 ) -
253 piecewise_ramp(
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/Waveforms.py in piecewise_ramp(step_func, t_off, t_currents, currents, n, eps)
230 if abs(const) > eps:
231 response += np.array(
--> 232 [fixed_quad(step_func, t, t+t0, n=n)[0] for t in time]
233 ) * const
234 return response
/Users/nfoks/Codes/gitRepos/simpegEM1D_original/simpegEM1D/Waveforms.py in <listcomp>(.0)
230 if abs(const) > eps:
231 response += np.array(
--> 232 [fixed_quad(step_func, t, t+t0, n=n)[0] for t in time]
233 ) * const
234 return response
/Users/nfoks/Codes/gitRepos/scipy/scipy/integrate/quadrature.py in fixed_quad(func, a, b, args, n)
93 # print((b-a)/2.0 * np.sum(w*func(y, *args), axis=-1))
94
---> 95 return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
96
97
/Users/nfoks/Codes/gitRepos/scipy/scipy/interpolate/polyint.py in __call__(self, x)
81 print('__call__')
82 print(x)
---> 83 y = self._evaluate(x)
84 return self._finish_y(y, x_shape)
85
/Users/nfoks/Codes/gitRepos/scipy/scipy/interpolate/interpolate.py in _evaluate(self, x_new)
662 y_new = self._call(self, x_new)
663 if not self._extrapolate:
--> 664 below_bounds, above_bounds = self._check_bounds(x_new)
665 if len(y_new) > 0:
666 # Note fill_value must be broadcast up to the proper size
/Users/nfoks/Codes/gitRepos/scipy/scipy/interpolate/interpolate.py in _check_bounds(self, x_new)
697 print('_check_bounds')
698 print(x_new)
--> 699 raise ValueError("A value in x_new is above the interpolation "
700 "range.")
701
ValueError: A value in x_new is above the interpolation range.