question: How can we get the constraint matrix 'C' in inversion and which the inverse subproblem is used?

[makeabhishek]

1. How can we get constraint matrix while performing traveltime inversion. As shown below in the equation ''C' can be a smoothening operator. How to show that in pyGimli? $$\Delta m^k = \frac{J^T D^T D (d-f(m^k))-\lambda C^TC(m^k - m^0)}{(J^T D^TD J + \lambda C^T C))}$$

2. What is reference model $m^0$ in the above equation and how do we choose that?

3. Also, for inversion pygimli use a Gauss-Newton method and for solving the inverse subproblem, either the conjugate gradient based CGLS (default) or LSQR algorithm is used. How do we know which subproblem solver is used out of CGLS and LSQR?

Reference: Günther, Thomas, Carsten Rücker, and Klaus Spitzer. "Three-dimensional modelling and inversion of DC resistivity data incorporating topography—II. Inversion." Geophysical Journal International 166.2 (2006): 506-517.

 --------------------------------------------------------------------------------
  Date: Thu Mar 28 16:58:01 2024 Mountain Daylight Time

                OS : Windows
            CPU(s) : 64
           Machine : AMD64
      Architecture : 64bit
               RAM : 511.9 GiB
       Environment : Jupyter

  Python 3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:29:04) [MSC
  v.1929 64 bit (AMD64)]

           pygimli : 1.4.3
            pgcore : 1.4.0
             numpy : 1.25.0
        matplotlib : 3.7.2
             scipy : 1.11.2
           IPython : 8.12.0
           pyvista : 0.43.2
--------------------------------------------------------------------------------

Operating system: e.g. Windows? Python version: Python 3.9.18 pyGIMLi version: 1.4.6 Way of installation: Conda package

velInv = mgr.invert(data, mesh=mesh, lam=10, useGradient=1, zWeight=15.0, secNodes=5, absoluteError=error,
                    blockyModel = True, robustData = False,
                    lambdaFactor=1,
                    cType=2,
                    verbose=True
                   )

[halbmy]

1. The matrix can be accessed (through the forward operator), printed (in your case a SparseMapMatrix) and shown by

   C=mgr.fop.constraintMatrix()
   print(C)
   pg.show(C, markersize=1)

2. By default, the reference model does not exist (or is zero) but can be set by using invert(startModel=model, isReference=True)
3. By default, GCLS is used but can be switched to LSQR by using the LSQRInversion framework. The solver yields the same model update but opens some new possibilities, e.g. to add parameter constraints, and will be the default in the future.

[makeabhishek]

I tried to obtain constraint matrix but its giving me following error AttributeError: 'TravelTimeDijkstraModelling' object has no attribute 'constraintMatrix'

Similarly for LSQRInversion I did

from pygimli.frameworks.lsqrinversion import LSQRInversion
pg.info("Traveltime Inversion")
velInv = mgr.LSQRInversion(data, mesh=mesh, lam=20, useGradient=1, zWeight=15.0, secNodes=5, absoluteError=error,
                    blockyModel = True, robustData = True,
                    lambdaFactor=1,
                    cType=1,
                    startModel=1/2200, isReference=True,
                    #verbose=True
                   )
mgr.inv.echoStatus()
print(mgr.inv.chi2History )
plt.plot(mgr.inv.chi2History)

[halbmy]

Sorry, I made a mistake, it is fop.constraints() but constraintMatrix is the better name, so we will create an alias for it.

[makeabhishek]

Thanks, Does x and y-axisof constraint matrix shows the velocity values?

Other than that when I run inversion it prompt following error because there is no error column in data. Can I specify error in percent while doing inversion. As default it is 3%, I did define absolute error absoluteError=error

ERROR - <class 'pygimli.physics.traveltime.TravelTimeManager.TravelTimeManager'>.checkError(C:\Users\376189\Documents\pyGimli\pygimli\gimli\pygimli\physics\traveltime\TravelTimeManager.py:97)

Is this the correct way of doing that? Can we specify both error and absolute error? data['err'] = mgr.estimateError(data['t'], errLevel=0.03, absError=1e-3)

[halbmy]

Please read some notes about inverse theory to learn about the constraint matrix. The model vector contains the (log-transformed) velocity vector and the constraint matrix is a smoothness operator.

Yes, the line you specified for the error estimation is correct. Instead, you can pass either relativeError or absoluteError to the inversion.