Dual moment fails at interpolation

[leonfoks]

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.

[leonfoks]

Example 2

Here is a different example for WalkTEM data.

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

mesh1D = set_mesh_1d(np.asarray([30.0, np.inf]))

depth = -mesh1D.gridN[:-1]
LocSigZ = -mesh1D.gridCC

time_LM = np.asarray([1.149E-05, 1.350E-05, 1.549E-05, 1.750E-05, 2.000E-05, 2.299E-05, 2.649E-05, 3.099E-05, 3.700E-05, 4.450E-05, 5.350E-05, 6.499E-05, 7.949E-05, 9.799E-05, 1.215E-04, 1.505E-04, 1.875E-04, 2.340E-04, 2.920E-04, 3.655E-04, 4.580E-04, 5.745E-04, 7.210E-04])
time_input_currents_LM = np.asarray([-8.333e-03, -8.033e-03, 0.0, 5.6000e-06])
input_currents_LM = np.asarray([0.0, 1.0, 1.0, 0.0])

time_HM = np.asarray([9.810E-05, 1.216E-04, 1.506E-04, 1.876E-04, 2.341E-04, 2.921E-04, 3.656E-04, 4.581E-04, 5.746E-04, 7.211E-04, 9.056E-04, 1.138E-03, 1.431E-03, 1.799E-03, 2.262E-03, 2.846E-03, 3.580E-03, 4.505E-03, 5.670E-03, 7.135E-03 ])
time_input_currents_HM = np.asarray([ -1.041e-03, -9.850e-04, 0.0,4.0000e-06])
input_currents_HM = np.asarray([0.0, 1.0, 1.0, 0.0])

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=22.5675833419,
            I=1.,
            time=time_HM,
            time_input_currents=time_input_currents_HM,
            input_currents=input_currents_HM,
            n_pulse=1,
            base_frequency=30.,
            use_lowpass_filter=True,
            high_cut_frequency=450000,
            moment_type='dual',
            time_dual_moment=time_LM,
            time_input_currents_dual_moment=time_input_currents_LM,
            input_currents_dual_moment=input_currents_LM,
            base_frequency_dual_moment=240,
          #  half_switch = True
        )

expmap = Maps.ExpMap(mesh1D)
m_1D = np.log(np.asarray([10.0, 1.0]))
chi = np.zeros(TDsurvey.n_layer)

prob = EM1D(
    mesh1D, sigmaMap=expmap, chi=chi
)

prob.pair(TDsurvey)
dBzdtTD = TDsurvey.dpred(m_1D)

[lheagy]

Hi @leonfoks: thanks for these issues! @sgkang is away at the moment (he gets back next week), and he is the most well-suited to address these. Sorry for the delay!

[sgkang]

I'll reply this asap.

[sgkang]

This was due to a bug when computing time range for interpolation. It is fixed now in #22