Set measurements in crosshole tomography

[GGDRriedel]

Problem description

Loaded in my own S/G Geometry for a crosshole tomography Forward simulation works fairly well.

Trying to set picked times in crosshole tomography

Your environment

Operating system: Win 7 Python version: 3.6 pyGIMLi version: 0.11 Way of installation: Conda Wheel

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:

# Most stuff taken from example

rays=measurements[['s','g']].values

# Empty container
scheme = pg.DataContainer()

for sen in sensorarray:
    scheme.createSensor(sen)

scheme.resize(len(rays))
scheme.add("s", rays[:, 0])
scheme.add("g", rays[:, 1])
scheme.add("valid", np.ones(len(rays)))
scheme.registerSensorIndex("s")
scheme.registerSensorIndex("g")

tt = Refraction()
mesh_fwd.createSecondaryNodes(5)
data = tt.simulate(mesh=mesh_fwd, scheme=scheme, slowness=1. / model,
                   noiseLevel=0.001, noiseAbs=1e-5)

########### CRUCIAL PART:
data.set('t',list(dt['t'].values)) # dt contains basically the sgt file positions and picks like s g t as pandas df

ttinv = Refraction()
ttinv.setData(data)
ttinv.setMesh(mesh, secNodes=5)
invmodel = ttinv.invert(lam=1100, vtop=2000, vbottom=2000, zWeight=1.0,maxIter=1,verbose=True)

ttinv.showResult()
...

Expected behavior

Inversion running with new t values

Actual behavior

Throws an error

C:/msys64/home/Guenther.T/src/gimli/gimli36/src/inversion.h: 638        double GIMLI::Inversion<ModelValType>::getPhiD(const Vec&) const [with ModelValType = double; GIMLI::Inversion<ModelValType>::Vec = GIMLI::Vector<double>]  getPhiD == inf

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

I guess the error occurs in the ttinv.invert() step? So we can focus on that part. Can you attach the data file (data.save(filename, 's g t')) and the mesh (or code for generating it), so that we can try to reproduce it? Maybe it has already been solved in a more recent version.

[GGDRriedel]

Indeed! crosshole.py.txt datafile.#dt.txt dtreader.py.txt thefile.sgt.txt

I am generating the data from the #dt file, the .sgt is the one generated by dat.save.

Code is a bit messy since I am iteratively working from your example

[halbmy]

I cannot run the code due to missing dependencies (dtreader), and it is indeed too messy. Can you please just attach the mesh (`mesh.save() or the x and y vectors) so that we can just use the last 5 lines of the code and not the messy part at the beginning.

[GGDRriedel]

dtreader is in the files but nevertheless here is the mesh themesh.bms.txt

[GGDRriedel]

I noticed that in my setup the scheme has 804 data points while the data coming out of the simulation has 777

[halbmy]

can you add a screenshot of how the synthetic model looks like?

[halbmy]

No idea whether there are some invalid data reducing the number of data. Anyway, the following script

import matplotlib as mpl
import pygimli as pg
from pygimli.physics.traveltime import TravelTimeManager as Refraction  # pg 1.1
# from pygimli.physics.traveltime import Refraction # pg1.0
mpl.rcParams['image.cmap'] = 'inferno_r'

data = pg.DataContainer('thefile.sgt', "s g")
mesh = pg.Mesh('themesh.bms')

ttinv = Refraction()
ttinv.setData(data)
ttinv.setMesh(mesh, secNodes=5)
invmodel = ttinv.invert(lam=1100, vtop=2000, vbottom=2000,
                        zWeight=1.0, maxIter=1, verbose=True)

ax, _ = ttinv.showResult()

works (using Python 3.7 from Anaconda on Windows) using pyGIMLi 1.1: but not in pyGIMLi 1.0.12.1 (leading to the mentioned PhiD error) I'll check that in detail.

[GGDRriedel]

Mesh and sensors.

It's probably fixed in a later version since most errors I get seem to be of numerical nature.

I'll try to update to the branch and come back with info.

EDIT: Just noticed that it's not released yet and I do not have the time to go through compilation and so on so I'll have to wait it out. Thanks anyway, maybe there's still something funky going on with the data.

[florian-wagner]

@GGDRriedel Release candidates of pygimli 1.1 are available in the test channel: https://anaconda.org/gimli/pygimli/files?version=1.1.0rc5

conda install -c gimli/label/test -c conda-forge pygimli=1.1.0rc5

[halbmy]

Unlike in the crosshole traveltime example https://www.pygimli.org/_examples_auto/seismics/plot_3_crosshole_tomography.html your sensors are no nodes in the mesh. While in pyGIMLi 1.1 this seems to be accounted for by choosing the nearest node, in version 1.0 this seems to lead to inaccurate Jacobian matrix. We will not fix this for v1.0 anymore.

Anyway, I guess node-free meshes will lead to inaccuracies also in v1.1. Therefore we recommend (particularly for real-world cases) creating a triangular mesh with all sensor positions in a box instead of a rectangular mesh. This works for both v1.0 and v1.1.

A reasonable synthetic model should have some structure (otherwise an infinite regularization parameter will lead to a constant model that fits the data) and also some Gaussian noise (with its std used as error) so that the regularization parameter can be chosen accordingly. I see nothing more to do here and close the issue.

[GGDRriedel]

What do you mean by "in a box" ?

[halbmy]

A box is a rectangle bounding the domain like in the crosshole example https://www.pygimli.org/_examples_auto/seismics/plot_3_crosshole_tomography.html where all sensors are added as nodes. Please have a look at that.