Half space option causes TDsurvey.dpred() fail

[leonfoks]

I may not be doing this correctly...

Ive attached a notebook cell, and the error messages below.

The only difference is half_switch = True in the instantiation of TDsurvey. Looks like the ExpMap is confused about the 19 layer model, and the request for only the results due to one "layer".

Perhaps a separate issue, but related. I have not been able to create a single layer mesh to represent a half space instead of a multi-layer mesh with a half_switch logical.

It would be a good idea to be able to pass a single layer model to the forward and sensitivity methods. Inside simply use if (half_switch or n_layer == 1): instead of forcing the user to specify the half_switch argument?

Here's a 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
from geobipy import StatArray

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_HM, time_LM]), facter_tmax=0.5, factor_tmin=10.,
    n_layer=19
)

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_HM,
            time_input_currents=time_input_currents_HM,
            input_currents=input_currents_HM,
            n_pulse=2,
            base_frequency=25.,
            use_lowpass_filter=True,
            high_cut_frequency=7e4,
            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=210,
            half_switch = True
        )

sig_half=1e-2
chi_half=0.

expmap = Maps.ExpMap(mesh1D)
m_1D = np.log(np.arange(TDsurvey.n_layer)*sig_half)
chi = np.zeros(TDsurvey.n_layer)

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

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

This causes

<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/anaconda/lib/python3.5/site-packages/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/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in fields(self, m)
    422     # @profile
    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

/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in forward(self, m, output_type)
    301 
    302         # TODO: potentially store
--> 303         sig = self.sigma_cole()
    304 
    305         if output_type == 'response':

/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in sigma_cole(self)
    228         )
    229 
--> 230         sigma = np.tile(self.sigma.reshape([-1, 1]), (1, n_frequency))
    231         if np.isscalar(self.eta):
    232             eta = self.eta

/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Props.py in fget(self)
    211                     )
    212                 )
--> 213             return mapping * self.model
    214 
    215         def fset(self, value):

/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Maps.py in __mul__(self, val)
    192                 raise ValueError(
    193                     'Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(
--> 194                         str(self), str(val.shape)
    195                     )
    196                 )

ValueError: Dimension mismatch in ExpMap(19,19) and np.ndarray(1,).

[sgkang]

Hi @leonfoks, That's a fair point, and I can make some modifications so that you can basically put a single layer model without half_switch.

[sgkang]

Hi @leonfoks, error occurs because of hz, you need to set the size of this as 1 for halfspace. It seems annoying to have halfspace option. I'll remove this option later.

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
# from geobipy import StatArray

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]

mesh1D = set_mesh_1d([1])
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=np.r_[0.],
            rx_type='dBzdt',
            wave_type='general',
            src_type='CircularLoop',
            a=13.,
            I=1.,
            time=time_HM,
            time_input_currents=time_input_currents_HM,
            input_currents=input_currents_HM,
            n_pulse=2,
            base_frequency=25.,
            use_lowpass_filter=True,
            high_cut_frequency=7e4,
            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=210,
            half_switch = True
        )

sig_half=1e-2
chi_half=0.

expmap = Maps.ExpMap(mesh1D)
m_1D = np.log(np.arange(TDsurvey.n_layer)*sig_half)
chi = np.zeros(TDsurvey.n_layer)

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

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