Possible bug in BaseRx.eval() when using resistivity

DC and IP
1

Hi,
Using SimPEG 0.21.1 on Ubuntu.

I’m trying to write a 3D resistivity inversion script, basing my code off the various tutorials. When I run it, it starts the actual inversion then stops with an error message as below

                   SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same Solver, and solver_opts as the Simulation3DNodal problem***

model has any nan: 0
============================ Inexact Gauss Newton ============================

beta phi_d phi_m f |proj(x-g)-x| LS Comment

x0 has any nan: 0

ValueError: shape mismatch: value array of shape (289,) could not be broadcast to indexing result of shape (17,)

My first Source electrode pair has 17 receivers. Note that 17*17=289.

An abridged stack trace is below:

File ~/.virtualenvs/simpeg/lib/python3.10/site-packages/SimPEG/simulation.py:291, in BaseSimulation.dpred(self, m, f)
289 for src in self.survey.source_list:
290 for rx in src.receiver_list:
→ 291 data[src, rx] = rx.eval(src, self.mesh, f)
292 return mkvc(data)

File ~/.virtualenvs/simpeg/lib/python3.10/site-packages/SimPEG/data.py:332, in Data.setitem(self, key, value)
330 def setitem(self, key, value):
331 index = self.index_dictionary[key[0]][key[1]]
→ 332 self.dobs[index] = mkvc(value)

Tracing back, the error is thrown when an attempt is made to assign 289 values into 17 array locations in dobs. The 289 values are generated in BaseRx.eval() by the line

return v / self.geometric_factor[src]

‘v’ is a 2 dimesional array of shape (17, 1) in this case, while geometric_factor is a 1 dimensional array of shape (17,). I’m not a python expert, but I assume that each row of v is being divided by the entire geometric_factor array, therefore producing 289 values.

I’ve debugged the tutorial scripts and examined survey and dobs to see if I’ve done something wrong with assembling my sources and receivers to create a survey or with the dobs, but the format looks to be identical.

The source_list is being assembled manually, then used to create a Survey, then set_geometric_factor() is called as below:

source_list = []
dobs_res_list = []

dfgg = dfg.groupby(['C1X', 'C1Y', 'C2X', 'C2Y'])

for group, frame in dfgg:
    A_location = np.r_[group[0], group[1], 0.0]
    B_location = np.r_[group[2], group[3], 0.0]
    M_electrodes = frame[['P1X', 'P1Y', 'P1Z']]
    N_electrodes = frame[['P2X', 'P2Y', 'P2Z']]
    
    fsc.sampleRaster(topo_path, A_location, 2)
    fsc.sampleRaster(topo_path, B_location, 2)
    
    M_locations = M_electrodes.to_numpy()
    N_locations = N_electrodes.to_numpy()
    fsc.sampleRaster(topo_path, M_locations, 2)
    fsc.sampleRaster(topo_path, N_locations, 2)
    
    dobs_res_list.extend(frame.Res.to_list())
    
    receiver_list = dc.receivers.Dipole(M_locations, N_locations, data_type=data_type)
    src = dc.sources.Dipole(receiver_list, A_location, B_location)
    source_list.append(src)

dobs_res = np.array(dobs_res_list, dtype='float64')
dobs_cha = np.array(dobs_cha_list, dtype='float64')


# Define survey
survey_res = dc.Survey(source_list)

# Set geometric factor for survey (to be handled internally in future versions)
survey_res.set_geometric_factor()

data_res = data.Data(survey_res, dobs=dobs_res)

Thanks for your help and please let me know if it’s something I’ve done.

Regards,
Roland

2

Using flatten() on the ‘v’ array appears to fix it. I’ll send a PR when I get everything else working to confirm.

return v.flatten() / self.geometric_factor[src]

Regards,
Roland