Extract the gradient used for the optimisation on each iteration

[lfrincond]

import numpy as np
import pygimli as pg
from pygimli.physics import ert

res = np.load('/syn_mod_squared.npy') #a resistivity model of shape (15,104)
real_x, real_y = np.shape(res)[1]*2, np.shape(res)[0]*2
elec = 27
ac_array = 'dd'

mesh = pg.createGrid(np.linspace(0, real_x, np.shape(res)[1] + 1), np.linspace(-real_y, 0, np.shape(res)[0] + 1))
mesh.cellMarkers = res.flatten()


scheme = ert.createData(elecs=np.linspace(start=0, stop=res.shape[1]*2, num=elec), schemeName=ac_array);
d_obs = ert.simulate(mesh, scheme=scheme, res=mesh.cellMarkers, noiseLevel=1, noiseAbs=1e-6, seed=1337)


#starting_model
st_model = np.ones(np.shape(res))*np.mean(res)
model_reg= 3.6

mgr = ert.ERTManager(d_obs)
inv = mgr.invert(d_obs, mesh=mesh, lam=model_reg, verbose=True, startModel=st_model, maxIter=0)

I would like to get from the inv object the computed gradient used for the optimisation.. I would like to see how the model and data regularisations are changing this gradient on each iteration.

Thank you!

[halbmy]

In a Gauss-Newton inversion, the gradient is not explicitly used (like in other methods like NLCG).

For doing the computations, you don't have to transform the Jacobian into a numpy matrix, but you can use fop.jacobian().transMult(residual). However, this does not yet take the error weighting and data and model transformations into account. I will add this functionality in the code and share it here.

[halbmy]

For accessing any information in the course of inversion, you can define a function to be called before (setPreStep) or after (setPostStep) each iteration:

def preStep(iter, inv):
    # do something, e.g.
    inv.fop.jacobian()...
    # save some stuff

inv = pg.Inversion(fop=fop)
inv.setPreStep(preStep)  # and postStep alike

[halbmy]

The inversion instance has now the methods gradient()
https://dev.pygimli.org/pygimliapi/_generated/pygimli.frameworks.html#pygimli.frameworks.Inversion.gradient
but also modelGradient() and dataGradient() as well as jacobianMatrix().

[halbmy]

This is all included in the last version v1.5.2. We will also prepare a tutorial about this, based on the notebook https://github.com/gimli-org/notebooks/blob/main/inversion/gradientVES.ipynb
I am closing this. In case of problems or questions, reopen or raise a new issue.