Errors with setRegularization in Traveltime inversion

[carmsignes]

Hello! Recently I raised a similar issue. I tried implementing the region treatment from the region example, but adding another structure to the parametric mesh for the inversion causes the inversion to crash. My objective is to fix the sound speed in a specific obstacle in the medium for the inversion, as it is known a priori. This will hopefully also give me a better resolution in the rest of the medium. Here is my code:

import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
from pygimli import meshtools as mt
from pygimli.physics import ert
from pygimli.physics import traveltime as tt
from pygimli.physics.petro import transFwdArchieS as ArchieTrans
from pygimli.physics.petro import transFwdWyllieS as WyllieTrans
from pygimli.frameworks import (PetroInversionManager,
                                JointPetroInversionManager)



N_sensors = 256
tank_radius = 1
n_Segments = 1024

### Create initial geometry and nodes

tank = mt.createCircle(radius=tank_radius,nSegments=n_Segments, marker=1, BoundaryMarker=-1)

water_layer_1 = mt.createCircle(radius=0.25,nSegments=n_Segments, marker=2)
water_layer_2 = mt.createCircle(radius=0.50,nSegments=n_Segments, marker=3)
#pipe = mt.createCircle(pos=[0.26,0.26] ,radius=0.1, nSgments=n_Segments*0.5, marker=4, markerPosition=[0.32,0.32]) + mt.createCircle(pos=[0.26,0.26] ,radius=0.07, nSgments=n_Segments*0.5, marker=5, markerPosition=[0.25,0.25]) 
pipe= mt.createCircle(pos=[0.26,0.26] ,radius=0.1, nSgments=n_Segments*2, marker=4, markerPosition=[0.32,0.32])
geom = tank + water_layer_2  + pipe + water_layer_1 


c = mt.createCircle(radius = 0.99 * tank_radius, nSegments=N_sensors,start=np.pi,end=np.pi*3)
for p in c.nodes():
    geom.createNode(p.pos(),-99)
    
mMesh = mt.createMesh(geom, q=44, area=0.0005, smooth=[1, 10])



# Create the parametric mesh that only reflects the domain geometry.
world1 = mt.createCircle(nSegments=n_Segments, area=0.0005, marker=2)
world=world1 + mt.createCircle(pos=[0.26,0.26] ,radius=0.1, nSgments=n_Segments*2, marker=1, markerPosition=[0.32,0.32])
    
pMesh = pg.meshtools.createMesh(world, q=44, smooth=[1, 10])


speed_of_sound = pg.solver.parseArgToArray([[1,1482],[2,1550],[3,1528], [4,4000]], mMesh.cellCount(),mMesh)
speed_of_sound_p = pg.solver.parseArgToArray([[1,1482],[2,1550],[3,1528],[4,4000]], pMesh.cellCount(),pMesh)

##Show the initial geometry and the mesh for the inversion

velKW = dict(cMin=1400, cMax=1600, logScale=False, cMap="plasma_r")


ax, cb = pg.show(mMesh,speed_of_sound,**velKW, showMesh=False, label='Speed of sound [m/s]')

ax, cb = pg.show(pMesh,speed_of_sound_p,**velKW, showMesh=True, label='Speed of sound [m/s]')
sensors = mMesh.positions()[mMesh.findNodesIdxByMarker(-99)]



pg.info("Simulate Traveltime")
ttScheme = tt.createRAData(sensors)
ttData = tt.simulate(mMesh, scheme=ttScheme, vel=speed_of_sound,
                     noiseLevel=0.001, noiseAbs=4e-7)

ax, cb = tt.showVA(ttData, **velKW)

pg.show(pMesh)


pg.info("Traveltime Inversion")
TT = tt.TravelTimeManager(verbose=True)
TT.setMesh(pMesh)
TT.inv.setRegularization('*', background=False)
TT.inv.setRegularization(1, fix=4000, single=True)
velInv = TT.invert(ttData, mesh=pMesh, lam=10000, useGradient=0, zWeight=1.0, verbose=True)
TT.inv.echoStatus()
ax, _ = pg.show(pMesh, velInv, **velKW, showMesh=False, label='Reconstructed speed of sound [m/s]')




velKW = dict(cMin=1450, cMax=1800, logScale=False, cMap="plasma_r")
ax, _ = pg.show(pMesh, velInv, **velKW, showMesh=False, label='Reconstructed speed of sound [m/s]')

And this is the resulting error:

%runfile C:/Users/carlos/Desktop/pyekfmm/gimli/pygimli/testing/tank_nopipe.py --wdir
28/04/25 - 10:30:18 - pyGIMLi - WARNING - parseMapToCellArray: cannot find marker 3 within mesh.
28/04/25 - 10:30:18 - pyGIMLi - WARNING - parseMapToCellArray: cannot find marker 4 within mesh.
28/04/25 - 10:30:19 - pyGIMLi - INFO - Simulate Traveltime
28/04/25 - 10:30:19 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
R3Vector: n=256
Reset region parameter
RegionManager copying mesh ...0.0552822 s 
create NeighborInfos ... 2.62e-05 s 
analysing mesh ... 2 regions.
creating para domain ... 0.025667 s
28/04/25 - 10:33:26 - pyGIMLi - INFO - Traveltime Inversion
28/04/25 - 10:33:26 - pyGIMLi - INFO - Found 2 regions.
creating para domain ... 0.0343351 s
creating para domain ... 0.0344297 s
creating para domain ... 0.0336631 s
Reset region parameter
RegionManager copying mesh ...0.0514978 s 
create NeighborInfos ... 2.35e-05 s 
analysing mesh ... 2 regions.
creating para domain ... 0.0244273 s
28/04/25 - 10:33:26 - pyGIMLi - INFO - Region with smallest marker set to background (marker=1)
28/04/25 - 10:33:26 - pyGIMLi - INFO - Found 2 regions.
28/04/25 - 10:33:26 - pyGIMLi - INFO - Region with smallest marker set to background (marker=1)
creating para domain ... 0.035943 s
28/04/25 - 10:33:26 - pyGIMLi - INFO - Creating forward mesh from region infos.
28/04/25 - 10:33:27 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
ModellingBase::setMesh() copying new mesh ... 0.292803 s
FOP updating mesh dependencies ... 2.49e-05 s
min/max(dweight) = 580499/2.40156e+06
28/04/25 - 10:33:27 - pyGIMLi - INFO - Create constant starting model: 0.0006789466395558284
28/04/25 - 10:33:27 - pyGIMLi - INFO - Created startmodel from forward operator:14858, min/max=0.000679/0.000679
28/04/25 - 10:33:27 - pyGIMLi - INFO - Starting inversion.
use model trans from RegionManager
min/max(dweight) = 580499/2.40156e+06
Building constraints matrix
constraint matrix of size(nBounds x nModel) 21763 x 14858
check Jacobian: wrong dimensions: (0x0) should be (65280x14858)  force: 1
jacobian size invalid, forced recalc: 1
Calculating Jacobian matrix (forced=1)...*** ./core/src/ttdijkstramodelling.cpp:500	
Warning: There are cells with marker -1. Did you define a boundary region? (This is not needed).
fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x000002A9ABD534C0>
Data transformation: Identity transform
Model transformation (cumulative):
	 0 Logarithmic LU transform, lower bound 0.0, upper bound 0.0
min/max (data): 1.5e-05/0.0013
min/max (error): 0.13%/2.69%
min/max (start model): 6.8e-04/6.8e-04
--------------------------------------------------------------------------------
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
File c:\users\carlos\desktop\pyekfmm\gimli\pygimli\testing\tank_nopipe.py:99
     96 TT.inv.setRegularization('*', background=False)
     97 #TT.inv.setRegularization(4, limits=[3999,4001], fix=NUMBER)
     98 #TT.inv.setRegularization(5, limits=[339,341], fix=NUMBER)
---> 99 velInv = TT.invert(ttData, mesh=pMesh, lam=10000, useGradient=0, zWeight=1.0, verbose=True)
    100 TT.inv.echoStatus()
    101 ax, _ = pg.show(pMesh, velInv, **velKW, showMesh=False, label='Reconstructed speed of sound [m/s]')

File ~\.conda\envs\pg\Lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py:268, in TravelTimeManager.invert(self, data, useGradient, vTop, vBottom, secNodes, **kwargs)
    265     self.fop._useGradient = None
    267 ### invert return mapped models
--> 268 slowness = super().invert(data, mesh, **kwargs)
    269 velocity = 1.0 / slowness
    270 velocity.isParaModel = slowness.isParaModel

File ~\.conda\envs\pg\Lib\site-packages\pygimli\frameworks\methodManager.py:852, in MeshMethodManager.invert(self, data, mesh, startModel, **kwargs)
    846 # if "blockyModel" in kwargs:  # redundant, now in fw
    847 #     self.fw.blockyModel = kwargs["blockyModel"]
    848 # if "stopAtChi1" in kwargs:
    849 #     self.fw._stopAtChi1 = kwargs["stopAtChi1"]
    851 self.preRun(**kwargs)
--> 852 self.fw.run(dataVals, errorVals, **kwargs)
    853 self.postRun(**kwargs)
    854 return self.paraModel(self.fw.model)

File ~\.conda\envs\pg\Lib\site-packages\pygimli\frameworks\inversion.py:1422, in ClassicInversion.run(self, dataVals, errorVals, **kwargs)
   1420 # self.inv.start()  # start is reset() and run() so better run?
   1421 self.inv.setMaxIter(0)
-> 1422 self.inv.start()
   1423 self.maxIter = maxIterTmp
   1424 if self.verbose:

File ~\.conda\envs\pg\Lib\site-packages\pygimli\physics\traveltime\modelling.py:93, in TravelTimeDijkstraModelling.createJacobian(self, par)
     90 if not self.mesh():
     91     pg.critical("no mesh")
---> 93 return self._core.createJacobian(par)

RuntimeError: ./core/src/sparsematrix.h:342		GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement GIMLI::SparseMapMatrix<ValueType, IndexType>::Aux::operator[](IndexType) [with ValueType = double; IndexType = long long unsigned int; GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement = GIMLI::MatrixElement<double, long long unsigned int, std::map<std::pair<long long unsigned int, long long unsigned int>, double, std::less<std::pair<long long unsigned int, long long unsigned int> >, std::allocator<std::pair<const std::pair<long long unsigned int, long long unsigned int>, double> > > >]  idx = 18446744073709551615, 116 maxcol = 14858 stype: 0

I have tried multiple fixes based on the region example and other discussions in the GitHub, but I always get this error when introducing another structure to the inversion mesh.

Thanks for your help!

[halbmy]

I tried it with a reduced number of sensors to save time, but changed a few things. You set a mesh, set region properties, but then you again pass a mesh to the invert() call. If you avoid this, it should work:

TT = tt.TravelTimeManager(ttData, verbose=True)
TT.setMesh(pMesh)
TT.inv.setRegularization('*', background=False)
TT.inv.setRegularization(1, fix=4000)
velInv = TT.invert(lam=10000, useGradient=0, zWeight=1.0, verbose=True)

[carmsignes]

Thanks for the reply! Indeed I was overwriting the region properties by passing the mesh again. I introduced the changes you suggested which makes the inversion run and take into account the region properties, but I still encounter some errors:

The changes:

pg.info("Traveltime Inversion")
TT = tt.TravelTimeManager(ttData, verbose=True)
TT.setMesh(pMesh)
TT.inv.setRegularization('*', background=False)
TT.inv.setRegularization(1, fix=4000)
velInv = TT.invert(lam=1000, useGradient=0, zWeight=1.0, verbose=True)
TT.inv.echoStatus()
ax, _ = pg.show(pMesh, velInv, **velKW, showMesh=False, label='Reconstructed speed of sound [m/s]')

The error:

%runfile C:/Users/carlos/Desktop/tank_nopipe.py --wdir
30/04/25 - 16:52:56 - pyGIMLi - WARNING - parseMapToCellArray: cannot find marker 3 within mesh.
30/04/25 - 16:52:56 - pyGIMLi - WARNING - parseMapToCellArray: cannot find marker 4 within mesh.
30/04/25 - 16:52:57 - pyGIMLi - INFO - Simulate Traveltime
30/04/25 - 16:52:57 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
Reset region parameter
30/04/25 - 16:54:32 - pyGIMLi - INFO - Traveltime Inversion
RegionManager copying mesh ...0.045093 s 
create NeighborInfos ... 4.01e-05 s 
analysing mesh ... 2 regions.
creating para domain ... 0.0184423 s
30/04/25 - 16:54:32 - pyGIMLi - INFO - Found 2 regions.
30/04/25 - 16:54:32 - pyGIMLi - INFO - Region with smallest marker set to background (marker=1)
creating para domain ... 0.0261893 s
creating para domain ... 0.0266568 s
30/04/25 - 16:54:32 - pyGIMLi - INFO - Creating forward mesh from region infos.
30/04/25 - 16:54:32 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
ModellingBase::setMesh() copying new mesh ... 0.10002 s
FOP updating mesh dependencies ... 5.71e-05 s
min/max(dweight) = 581273/1.51038e+06
30/04/25 - 16:54:33 - pyGIMLi - INFO - Create constant starting model: 0.000679645063577222
30/04/25 - 16:54:33 - pyGIMLi - INFO - Created startmodel from forward operator:7406, min/max=0.000680/0.000680
30/04/25 - 16:54:33 - pyGIMLi - INFO - Starting inversion.
use model trans from RegionManager
Info: fixing region:  1  to:  4000
min/max(dweight) = 581273/1.51038e+06
Building constraints matrix
constraint matrix of size(nBounds x nModel) 10591 x 7406
check Jacobian: wrong dimensions: (0x0) should be (240x7406)  force: 1
jacobian size invalid, forced recalc: 1
Calculating Jacobian matrix (forced=1)...*** ./core/src/ttdijkstramodelling.cpp:500	
Warning: There are cells with marker -1. Did you define a boundary region? (This is not needed).
Info: fixing region:  1  to:  4000
11.1916 s
min data = 0.000260136 max data = 0.00132266 (240)
min error = 0.00130052 max error = 0.00254515 (240)
min response = 0.000264226 max response = 0.00136374 (240)
calc without reference model
0: rms/rrms(data, response) = 4.4521e-05/3.73422%
0: chi^2(data, response, error, log) = 775.906
0: Phi = 186217 + 0 * 1000 = 186217
solve CGLSCDWWtrans with lambda = 1000
fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x0000021A6CAFE6B0>
Data transformation: Identity transform
Model transformation (cumulative):
	 0 Logarithmic LU transform, lower bound 0.0, upper bound 0.0
min/max (data): 2.6e-04/0.0013
min/max (error): 0.13%/0.25%
min/max (start model): 6.8e-04/6.8e-04
--------------------------------------------------------------------------------
(ALL THE ITERATIONS, REMOVED FOR SPACE)
--------------------------------------------------------------------------------
inv.iter 9 ... tau = 0. Trying parabolic line search with step length 0.3Info: fixing region:  1  to:  4000
 ==> tau = -0.0408534
 tau < 0.03 ==> tau = 0.03
9: LS newModel: min = 0.000437427; max = 0.000729762
9: LS newResponse: min = 0.000260565; max = 0.00129869
9: rms/rrms(data, LS newResponse) = 7.81813e-06/0.650484%
9: chi^2(data, LS newResponse, error, log) = 23.5852
9: Phi = 5660.45+0.796987*1000=6457.43
Performing line search with tau = 0.03
Info: fixing region:  1  to:  4000
9: Model: min = 0.000449173; max = 0.000716937
9: Response: min = 0.000260597; max = 0.00130877
9: rms/rrms(data, Response) = 4.32945e-06/0.374919%
9: chi^2(data, Response, error, log) = 7.57795
9: Phi = 1818.71+0.603371*1000=2422.08
9: Model: min = 0.000449173; max = 0.000716937
9: Response: min = 0.000260597; max = 0.00130877
9: rms/rrms(data, Response) = 4.32945e-06/0.374919%
9: chi^2(data, Response, error, log) = 7.57795
9: Phi = 1818.71+0.603371*1000=2422.08
30/04/25 - 17:01:07 - pyGIMLi - ERROR - None._drawField(C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py:399)
Data size invalid
chi² =    7.58 (dPhi = 1.07%) lam: 1000.0
################################################################################
#                Abort criterion reached: dPhi = 1.07 (< 2.0%)                 #
################################################################################
Data:  7406 1394.823022465174 2226.3119928348174 False
Mesh:  Mesh: Nodes: 4226 Cells: 7426 Boundaries: 11651
30/04/25 - 17:01:08 - pyGIMLi - ERROR - None._drawField(C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py:399)
Data size invalid
cannot access local variable 'gci' where it is not associated with a value
Traceback (most recent call last):
  File "C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py", line 426, in showMesh
    gci, validData = _drawField(ax, mesh, data, kwargs)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py", line 406, in _drawField
    return gci, validData
           ^^^
UnboundLocalError: cannot access local variable 'gci' where it is not associated with a value
Mesh: Nodes: 4226 Cells: 7426 Boundaries: 11651
Data:  7406 1394.823022465174 2226.3119928348174 False
Mesh:  Mesh: Nodes: 4226 Cells: 7426 Boundaries: 11651
cannot access local variable 'gci' where it is not associated with a value
Traceback (most recent call last):
  File "C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py", line 426, in showMesh
    gci, validData = _drawField(ax, mesh, data, kwargs)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\carlos\.conda\envs\pg\Lib\site-packages\pygimli\viewer\showmesh.py", line 406, in _drawField
    return gci, validData
           ^^^
UnboundLocalError: cannot access local variable 'gci' where it is not associated with a value

I cannot map the produced inversion data. Any pointers for this error? Trying TT.inv.setRegularization(1, limits=[3999,4001]) produces an output I can plot by my previous code, but the region in which I set the limits is not really constrained (there are cells with values outside the limits). Did my code run for you @halbmy ?

[halbmy]

Since you use regions and a fixed region particularly (which is not part of the inversion) you cannot use the original mesh for plotting. Please use mgr.showResult() instead.

[halbmy]

Did this work? Can we close the issue?

[carmsignes]

Yes, this was the issue! For further reference, as I am simulating a medium where I have a priori knowledge of certain regions, providing them with a high slowness value provides a much clearer picture of the surrounding medium. In a experimental setting, of course, it also means you lose some of the possible travel paths.