Jacobian matrix has different size than number of elements

[mariosgeo]

I am running several forward models, and I want to extract the Jacobian for each of those. My example script is

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import numpy as np
import matplotlib.pyplot as plt
# from read_ares import read_ares
# from make_pyg import make_pyg

iter1=1
for i in range(0,iter1):
# Create forward model
    world = mt.createWorld(start=[0, 0], end=[160, -50], layers=[-1, -5],
                       worldMarker=True)

    block = mt.createCircle(pos=[80, -3.-i], radius=[4, 1], marker=4,
                        boundaryMarker=10, area=0.1)




    geom = world + block


    pg.show(geom)
   data = ert.load("line1.pyg") #some gradient array with 10801 measures

   for p in data.sensors():
       geom.createNode(p)
       geom.createNode(p - [0, 0.1])

   mesh = mt.createMesh(geom, quality=34)





    rhomap = [[1, 100.],
          [2, 75.],
          [3, 50.],
          [4, 150.]]

    # Take a look at the mesh and the resistivity distribution
    pg.show(mesh, data=rhomap, label=pg.unit('res'), showMesh=True)




    res_cell = np.array(mesh.cellMarkers())

   res_cell[res_cell==1]=100
   res_cell[res_cell==2]=75
   res_cell[res_cell==3]=50
   res_cell[res_cell==4]=150



    fop = ert.ERTModelling()
    fop.setData(data)
    fop.setMesh(mesh)


   fop.createJacobian(res_cell)
   jac=np.asarray(fop.jacobian())

   app_res =ert.simulate(mesh, scheme=data, res=rhomap, noiseLevel=1,
                    noiseAbs=1e-6, seed=1337,returnArray=True)

The issue is that the number of cells are 8158, but the Jacobian matrix has size of 10801x4838.
Could you advise what's the issue?

[mariosgeo]

Edit:
On the same code, why does the

app_res2=fop.response(res_cell)
app_res =ert.simulate(mesh, scheme=data, res=rhomap,returnArray=True)

returns different results?

[halbmy]

The forward operator sorts the mesh according to the cell markers (to ease inversion). So fop.response cannot easily be used unless you sort the values.

[mariosgeo]

Hi,
When I plot using

 pg.show(mesh, data=rhomap, label=pg.unit('res'), showMesh=True,cMap='rainbow')
 pg.show(mesh, data=res_cell, label=pg.unit('res'), showMesh=True,cMap='rainbow')

I get the same model. So I suppose they are sorted ;

[halbmy]

In addition, by default the region with the lowest marker is treated as background region, so you should disable this by fop.setRegionProperties('*', background=False). Then the matrix size should be ok.

[mariosgeo]

Hi Thomas,

Thanks for the reply, this has fixed the jacobian size.
Yet, since the forward response is still different from fop.respose and ert.simulate, how can I know that my staring model is correct?

Is there another way to get the jacobian using the ert module?

[halbmy]

No, you can only get the jacobian through the forward operator and this sorts the cells according to the marker, so you just have to create a corresponding resistivity vector.

[mariosgeo]

I think I am doing that. I get the cell markers from the mesh with the following

res_cell = np.array(mesh.cellMarkers())

res_cell[res_cell==1]=100
res_cell[res_cell==2]=75
res_cell[res_cell==3]=50
res_cell[res_cell==4]=150

and I throw it it to the fop

app_res2=fop.response(res_cell)

The difference between the fop (blue line) and ert (orange line), seems not completely off. Cloud that be related to the background resistivity value?

[halbmy]

In your case the cells will still be arbitrarily ordered. Can you try fop(res_cell[np.argsort(mesh.cellMarkers())])?

[mariosgeo]

If i do this, then it's a much-much better match (ha). in the next figure, I plot the app_res/app_res2. There is some small change 0.5%, not important for my application.

Since I am trying to understand it, you order the cells to be sorted, but only to the fop operator. This applies also to the Jacobian?
How can I now get the e.g. cell center of each cell, that refer to the jacobian?

[halbmy]

The forward operator is designed for inversion purposes. When working with a mesh using different markers/regions (besides a potential background region), the forward operator resorts the mesh according to the region number so that the regions are separated (and e.g. smoothness is only within the regions and not between them unless explicitly defined). Usually a synthetic simulation is not done with the fop.response() but with ert.simulate() which cares for different inputs like maps, simple vectors or full-length vectors.

[mariosgeo]

Thanks a lot. My question is answered. I will close this discussion.