gimli
gimli copied to clipboard
Published
20 hours ago •
gimli-org
gimli-org
Error message with `fop.modelTrans` in `LCInversion`
Problem description
Hello! I have a problem with the function modelTrans of the class LCModelling. It's giving me an error message but I'm not able to find the issue causing it.
Your environment
Thu Feb 08 14:06:55 2024 CET -- OS | Linux | CPU(s) | 8 | Machine | x86_64 Architecture | 64bit | RAM | 7.8 GiB | Environment | Jupyter Python 3.9.15 \| packaged by conda-forge \| (main, Nov 22 2022, 08:45:29) [GCC 10.4.0] numpy | 1.23.5 | scipy | 1.9.3 | numba | 0.56.4 empymod | 2.2.1 | IPython | 8.6.0 | matplotlib | 3.5.1Operating system: Linux
Python version: Python 3.9.15
pyGIMLi version: 1.3.0
Additional: Using empymod version '2.2.1'
Way of installation: Conda package
Steps to reproduce
The code is written in the following order: First several class are defined
- The first class
fop1Dis the EM forward function inempymodin 1D - The second class
FDEM1Dis the EM forward modelling framework in 1D - The third class is
LCModellingto perform the modelling using FDEM1D - The third class is
LCInversionto run the inversion scheme The setup are 10 1D models with 3 layers Then the classes are initialized and run
## LCI Inversion
import numpy as np
import pygimli as pg
import empymod as ep
from scipy.constants import mu_0
import matplotlib.pyplot as plt
# Example for 3 layered model
nLayers=3
# Set 1D forward function with empymod
class fop1d():
""" 1D FDEM response
Parameters
----------
sgm : array
Array of electrical conductivities [S/m], len(sigma) = nlay
thk : array
Array of thicknesses [m], len(thk) = nlay - 1
Freq : frequency of device [Hz]
coilOrient : array of strings
coil orientations: 'H' for horizontal coplanar,
'V' for vertical coplanar, 'P' for perpendicular
coilSpacing : array
separations of the transmitter and receiver coils [m]
height : float
height of the device with respect to ground [m]
Returns
-------
Array [OP, IP]
"""
def __init__(self, Freq, coilOrient, coilSpacing, height):
""" initialize class """
self.Freq = Freq
self.coilOrient = coilOrient
self.coilSpacing = coilSpacing
self.height = height
def em1d(self, sgm, thk):
# Source and receivers geometry [x, y, z]
source = [0, 0, -self.height]
receivers = [self.coilSpacing, np.zeros(len(self.coilSpacing)), -self.height]
# Depth and resistivity
res_air = 1e6
res = np.hstack((res_air, 1/sgm))
depth = np.hstack((0, np.cumsum(thk)))
# Empty array to store responses
OUT = []
if any(coilOrient == 'H'):
H_Hs = ep.dipole(source, receivers, depth, res, self.Freq, ab = 66, xdirect=None,
verb=0)*(2j * np.pi * self.Freq * mu_0)
H_Hp = ep.dipole(source, receivers, depth=[], res=[res_air], freqtime = self.Freq,
ab = 66, verb=0)*(2j * np.pi * self.Freq * mu_0)
op = (H_Hs/H_Hp).imag.amp()
ip = (H_Hs/H_Hp).real.amp()
OUT.append([op, ip])
if any(coilOrient == 'V'):
V_Hs = ep.dipole(source, receivers, depth, res, self.Freq, ab =55, xdirect=None,
verb=0)*(2j * np.pi * self.Freq * mu_0)
V_Hp = ep.dipole(source, receivers, depth=[], res=[res_air], freqtime=self.Freq,ab=55,
verb=0)*(2j * np.pi * self.Freq * mu_0)
op = (V_Hs/V_Hp).imag.amp()
ip = (V_Hs/V_Hp).real.amp()
OUT.append([op, ip])
if any(coilOrient == 'P'):
# Maybe put 0.1m in receiver offset
P_Hs = ep.dipole(source, receivers, depth, res, self.Freq, ab=46, xdirect=None,
verb=0)*(2j * np.pi * self.Freq * mu_0)
P_Hp = ep.dipole(source, receivers, depth=[], res=[res_air], freqtime= self.Freq,
ab=66, verb=0)*(2j * np.pi * self.Freq * mu_0)
op = (P_Hs/P_Hp).imag.amp()
ip = (P_Hs/P_Hp).real.amp()
OUT.append([op, ip])
return np.array(OUT).ravel()
# Class for 1D forward modelling with pygimli
class FDEM1D(pg.frameworks.Modelling):
""" Class forward modelling Frequency Domain EM data """
def __init__(self, nLayers=nLayers):
""" Initialize class and set survey configuration """
self.nLayers = nLayers
mesh = pg.meshtools.createMesh1DBlock(self.nLayers)
super().__init__()
self.setMesh(mesh)
def response(self, mod):
""" Compute response vector for a certain model [mod]
mod = [thickness_1, thickness_2, ..., thickness_n, sigma_1, sigma_2, ..., sigma_n]
"""
resp = fop1d.em1d(np.asarray(mod)[self.nLayers-1:self.nLayers*2-1], # sigma
np.asarray(mod)[:self.nLayers-1] # thickness
)
return resp
def response_mt(self, mod, i=0):
"""Multi-threaded forward response."""
return self.response(mod)
def createJacobian(self, mod, dx=1e-8):
""" compute Jacobian for a 1D model """
resp = self.response(mod)
n_rows = len(resp) # number of data values in data vector
n_cols = len(mod) # number of model parameters
J = self.jacobian() # we define first this as the jacobian
J.resize(n_rows, n_cols)
Jt = np.zeros((n_cols, n_rows))
for j in range(n_cols):
mod_plus_dx = mod.copy()
mod_plus_dx[j] += dx
Jt[j,:] = (self.response(mod_plus_dx) - resp)/dx # J.T in col j
for i in range(n_rows):
J[i] = Jt[:,i]
def createDefaultStartModel(self, mod):
""" create default start model 1D """
model_0 = pg.Vector([2, 2, 20/1000, 20/1000, 20/1000])
return model_0
# LC Modelling class
class LCModelling(pg.frameworks.LCModelling):
"""2D Laterally constrained (LC) modelling.
"""
def __init__(self, fop, **kwargs):
"""Parameters: fop class ."""
super().__init__(fop, **kwargs)
self._singleRegion = False
self._fopTemplate = fop
self._fopKwargs = kwargs
self._fops1D = []
self._mesh = None
self._nSoundings = 0
self._parPerSounding = 0
self._jac = None
self.soundingPos = None
def setDataBasis(self, **kwargs):
"""Set homogeneous data basis.
Set a common data basis to all forward operators.
If you want individual you need to set them manually.
"""
for f in self._fops1D:
f.setDataBasis(**kwargs)
def initModelSpace(self, nLayers):
"""Initialize model space."""
for i, f in enumerate(self._fops1D):
f.initModelSpace(nLayers)
def createDefaultStartModel(self, models):
"""Create default starting model."""
sm = pg.Vector()
for i, f in enumerate(self._fops1D):
sm = pg.cat(sm, f.createDefaultStartModel(models[i]))
return sm
def response(self, par):
"""Cut together forward responses of all soundings."""
mods = np.asarray(par).reshape(self._nSoundings, self._parPerSounding)
resp = pg.Vector(0)
for i in range(self._nSoundings):
r = self._fops1D[i].response(mods[i])
resp = pg.cat(resp, r)
return resp
def createJacobian(self, par):
"""Create Jacobian matrix by creating individual Jacobians."""
mods = np.asarray(par).reshape(self._nSoundings, self._parPerSounding)
for i in range(self._nSoundings):
self._fops1D[i].createJacobian(mods[i])
def initJacobian(self, dataVals, nLayers, nPar=1):
"""Initialize Jacobian matrix.
Parameters
----------
dataVals : ndarray | RMatrix | list
Data values of size (nSounding x Data per sounding).
All data per sounding need to be equal in length.
If they don't fit into a matrix use list of sounding data.
"""
nSoundings = len(dataVals)
self.createParametrization(nSoundings, nLayers=nLayers, nPar=nPar)
if self._jac is not None:
self._jac.clear()
else:
self._jac = pg.matrix.BlockMatrix()
self.fops1D = []
nData = 0
for i in range(nSoundings):
kwargs = {}
for key, val in self._fopKwargs.items():
if hasattr(val, '__iter__'):
kwargs[key] = val[i]
else:
kwargs[key] = val
f = None
if issubclass(self._fopTemplate, pg.frameworks.Modelling):
f = self._fopTemplate(**kwargs)
else:
f = type(self._fopTemplate)(self.verbose, **kwargs)
f.setMultiThreadJacobian(self._parPerSounding)
self._fops1D.append(f)
nID = self._jac.addMatrix(f.jacobian())
self._jac.addMatrixEntry(nID, nData, self._parPerSounding * i)
nData += len(dataVals[i])
self._jac.recalcMatrixSize()
#print("Jacobian size:", self._jac.rows(), self._jac.cols(), nData)
self.setJacobian(self._jac)
return self._jac
def constraint_matrix(self, dataVals, nLayers):
""" Create constraint matrix
Parameters
----------
dataVals : list of 1D pyGIMLi data vectors (nSounding x Data per sounding) [list]
nLayers : number of layers (for blocky inversion) [int]
"""
nSoundings = len(dataVals) # number of soundings
boundaries_thk = (nLayers - 1) * (nSoundings - 1) # inner mesh boundaries for thicknesses
boundaries_sig = nLayers * (nSoundings - 1) # inner mesh boundaries for conductivities
CM = np.zeros((boundaries_thk + boundaries_sig, (nLayers*2 - 1) * nSoundings))
h = -np.eye(1, nLayers * 2) + np.eye(1, nLayers * 2, k = (nLayers * 2 - 1))
for i in range(boundaries_thk + boundaries_sig):
CM[i, i:h.shape[1]+i] = h
print('Size of constraint matrix:', np.shape(CM))
self.CM = pg.utils.toSparseMatrix(CM) # convert to sparse pg matrix
return self.CM
def createWeight(self, dataVals, cWeight_1, cWeight_2, nLayers):
""" Create constraint weights (cWeights)
Blocky model : vertical constraint weights for both model parameter regions
Parameters
----------
dataVals : list of 1D pyGIMLi data vectors [list]
cWeight_1 : thickness constraint weight (blocky model) [float]
cWeight_2 : resistivity constraint weight (blocky model) [float]
nLayers : number of layers (for blocky inversion) [int]
"""
nSoundings = len(dataVals) # number of soundings
""" constraint weights for blocky model """
cWeight_thk = cWeight_1
cWeight_sig = cWeight_2
boundaries_thk = (nLayers - 1) * (nSoundings - 1)
boundaries_sig = nLayers * (nSoundings - 1)
cWeight_thk = pg.Vector(boundaries_thk, cWeight_thk)
cWeight_sig = pg.Vector(boundaries_sig, cWeight_sig)
self.cWeight = pg.cat(cWeight_thk, cWeight_sig)
print('shape of Constraint cWeight:', np.shape(self.cWeight))
def createConstraints(self):
""" create weighted constraint matrix """
self._CW = pg.matrix.LMultRMatrix(self.CM, self.cWeight) # , verbose = True)
self.setConstraints(self._CW)
def drawModel(self, ax, model, **kwargs):
"""Draw models as stitched 1D model section."""
mods = np.asarray(model).reshape(self._nSoundings,
self._parPerSounding)
pg.viewer.mpl.showStitchedModels(mods, ax=ax, useMesh=True,
x=self.soundingPos,
**kwargs)
# LC Inversion class
class LCInversion(pg.Inversion):
"""Quasi-2D Laterally constrained inversion (LCI) framework."""
def __init__(self, fop, **kwargs):
super(LCInversion, self).__init__(fop=fop, **kwargs)
def prepare(self, dataVals, errorVals, nLayers, **kwargs):
"""Prepare inversion with given data and error vectors."""
dataVec = pg.RVector3()
for d in dataVals:
dataVec = pg.cat(dataVec, d)
errVec = pg.RVector3()
for e in errorVals:
errVec = pg.cat(errVec, e)
self.fop.initJacobian(dataVals=dataVals, nLayers=nLayers, nPar=1)
# self.fop.initJacobian resets prior set startmodels
if self._startModel is not None:
self.fop.setStartModel(self._startModel)
rC = self.fop.regionManager().regionCount()
if kwargs.pop('disableLCI', False):
self.inv.setMarquardtScheme(0.7)
# self.inv.setLocalRegularization(True)
for r in self.fop.regionManager().regionIdxs():
self.fop.setRegionProperties(r, cType=0)
else:
# self.inv.stopAtChi1(False)
cType = kwargs.pop('cType', None)
if cType is None:
cType = [1] * rC
zWeight = kwargs.pop('zWeight', None)
if zWeight is None:
zWeight = [0.0] * rC
self.fop.setRegionProperties('*',
cType=cType,
zWeight=zWeight,
**kwargs)
self.inv.setReferenceModel(self.fop.startModel())
return dataVec, errVec
def run(self, dataVals, errorVals, nLayers, **kwargs):
"""Run inversion with given data and error vectors."""
lam = kwargs.pop('lam', 20)
dataVec, errVec = self.prepare(dataVals, errorVals, nLayers, **kwargs)
print('#'*50)
print(kwargs)
print('#'*50)
return super(LCInversion, self).run(dataVec, errVec, lam=lam, **kwargs)
# Set problem
Freq = 9000
coilOrient = np.array(['H', 'V', 'P'])
coilSpacing = np.array([2, 4, 8])
nLayers = 3 # careful: nlay must be equal in 1Dfop and LC
height = 0
# Initialize forward operator function with empymod
fop1d = fop1d(Freq, coilOrient, coilSpacing, height)
# Initialize forward modeller with pygimli
FOP = FDEM1D()
# True models
pos = 10 # number of soundings
x = np.linspace(0,pos,pos, endpoint=False)
sigm_1 = 10/1000
sigm_2 = 20/1000
sigm_3 = 40/1000
thk_1 = 1/5*x + 1
thk_2 = 2 - 1/7*x
models = np.zeros((pos, nLayers*2-1))
models[:,0] = thk_1
models[:,1] = thk_2
models[:,2] = sigm_1
models[:,3] = sigm_2
models[:,4] = sigm_3
# Simulate data for each 1D model
data = []
for p in range(pos):
data.append(FOP.response(models[p]))
# Define error vector
relativeError = np.ones_like(data)*1e-3
# Initialize LCModelling class
LC = LCModelling(FDEM1D)
# Initialize inversion class
inv = LCInversion(LC)
# Run inversion
models_est = inv.run(data, relativeError, nLayers=nLayers)
Expected behavior
Run the LCInversion
Actual behavior
When the inversion calls the LCModelling and initializeJacobian:
I obtain an error message related to the fop._regionProperties inside LCModelling class
##################################################
{}
##################################################
Traceback (most recent call last):
File "/home/mariacarrizo/repos/LCI-EMI/LCI_EM.py", line 401, in <module>
models_est = inv.run(data, relativeError, nLayers=nLayers)
File "/home/mariacarrizo/repos/LCI-EMI/LCI_EM.py", line 355, in run
return super(LCInversion, self).run(dataVec, errVec, lam=lam, **kwargs)
File "/home/mariacarrizo/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/frameworks/inversion.py", line 511, in run
self.inv.setTransModel(self.fop.modelTrans)
File "/home/mariacarrizo/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/frameworks/modelling.py", line 153, in modelTrans
self._applyRegionProperties()
File "/home/mariacarrizo/anaconda3/envs/pg/lib/python3.9/site-packages/pygimli/frameworks/modelling.py", line 349, in _applyRegionProperties
rMgr.region(rID).setConstraintType(vals['cType'])
Boost.Python.ArgumentError: Python argument types in
Region.setConstraintType(Region, list)
did not match C++ signature:
setConstraintType(GIMLI::Region {lvalue}, unsigned long type)
Hello!
I'm still trying to work out the the LCI Inversion, so now I'm using the pg.Inversion() framework instead of LCInversion as in the following code:
# Import libraries
import numpy as np
import pygimli as pg
import empymod as ep
from scipy.constants import mu_0
# Define forward function
def FDEM1D(sgm, thk):
""" 1D FDEM response of DUALEM842s
Parameters
----------
sgm : array
Array of electrical conductivities [S/m], len(sigma) = nlay
thk : array
Array of thicknesses [m], len(thk) = nlay - 1
Returns
-------
Array [OP, IP]
"""
# Set geometry
Freq = 9000
coilSpacing = [2, 4, 8]
pcoilSpacing = [2.1, 4.1, 8.1]
coilOrient = np.array(['H', 'V', 'P'])
height = 0
# Source and receivers geometry [x, y, z]
source = [0, 0, -height]
receivers = [coilSpacing, np.zeros(len(coilSpacing)), -height]
preceivers = [pcoilSpacing, np.zeros(len(pcoilSpacing)), -height]
# Depth and resistivity
res_air = 1e6
res = np.hstack((res_air, 1/sgm))
depth = np.hstack((0, np.cumsum(thk)))
# Empty array to store responses
OUT = []
if any(coilOrient == 'H'):
H_Hs = ep.dipole(source, receivers, depth, res, Freq, ab = 66, xdirect=None,
verb=0)*(2j * np.pi * Freq * mu_0)
H_Hp = ep.dipole(source, receivers, depth=[], res=[res_air], freqtime = Freq,
ab = 66, verb=0)*(2j * np.pi * Freq * mu_0)
op = (H_Hs/H_Hp).imag.amp()
ip = (H_Hs/H_Hp).real.amp()
OUT.append([op, ip])
if any(coilOrient == 'V'):
V_Hs = ep.dipole(source, receivers, depth, res, Freq, ab =55, xdirect=None,
verb=0)*(2j * np.pi * Freq * mu_0)
V_Hp = ep.dipole(source, receivers, depth=[], res=[res_air], freqtime=Freq,
ab=55, verb=0)*(2j * np.pi * Freq * mu_0)
op = (V_Hs/V_Hp).imag.amp()
ip = (V_Hs/V_Hp).real.amp()
OUT.append([op, ip])
if any(coilOrient == 'P'):
P_Hs = ep.dipole(source, preceivers, depth, res, Freq, ab=46, xdirect=None,
verb=0)*(2j * np.pi * Freq * mu_0)
P_Hp = ep.dipole(source, preceivers, depth=[], res=[res_air], freqtime= Freq,
ab=66, verb=0)*(2j * np.pi * Freq * mu_0)
op = (P_Hs/P_Hp).imag.amp()
ip = (P_Hs/P_Hp).real.amp()
OUT.append([op, ip])
return np.array(OUT).ravel()
# Forward Modelling in 1D
class FDEM1DModelling(pg.frameworks.Modelling):
def __init__(self, nlay=3):
self.nlay = nlay
mesh = pg.meshtools.createMesh1DBlock(nlay)
super().__init__()
self.setMesh(mesh)
def createStartModel(self, dataVals):
startThicks = 2
startSigma = 1/20
# layer thickness properties
self.setRegionProperties(0, startModel = startThicks, trans = 'log')
# electrical conductivity properties
self.setRegionProperties(1, startModel = startSigma, trans = 'log')
return super(FDEM1DModelling, self).createStartModel()
def response(self, par):
""" Compute response vector for a certain model [mod]
par = [thickness_1, thickness_2, ..., thickness_n, sigma_1, sigma_2, ..., sigma_n]
"""
resp = FDEM1D(np.asarray(par)[self.nlay-1:self.nlay*2-1], # sigma
np.asarray(par)[:self.nlay-1] # thickness
)
return resp
def response_mt(self, par, i=0):
"""Multi-threaded forward response."""
return self.response(par)
def createJacobian(self, par, dx=1e-4):
""" compute Jacobian for a 1D model """
resp = self.response(par)
n_rows = len(resp) # number of data values in data vector
n_cols = len(par) # number of model parameters
J = self.jacobian() # we define first this as the jacobian
J.resize(n_rows, n_cols)
Jt = np.zeros((n_cols, n_rows))
for j in range(n_cols):
mod_plus_dx = par.copy()
mod_plus_dx[j] += dx
Jt[j,:] = (self.response(mod_plus_dx) - resp)/dx # J.T in col j
for i in range(n_rows):
J[i] = Jt[:,i]
#print(self.jacobian())
#print(J)
#print(Jt)
def drawModel(self, ax, model):
pg.viewer.mpl.drawModel1D(ax = ax,
model = model,
plot = 'semilogx',
xlabel = 'Electrical conductivity (S/m)',
)
ax.set_ylabel('Depth in (m)')
# Define Lateral constraint modelling class
class LCModelling(pg.frameworks.LCModelling):
"""2D Laterally constrained (LC) modelling.
2D Laterally constrained (LC) modelling based on BlockMatrices.
"""
def __init__(self, fop, **kwargs):
"""Parameters: fop class ."""
super().__init__(fop, **kwargs)
def initJacobian(self, dataVals, nLay):
"""Initialize Jacobian matrix.
Parameters
----------
dataVals : ndarray | RMatrix | list
Data values of size (nSounding x Data per sounding).
All data per sounding need to be equal in length.
If they don't fit into a matrix use list of sounding data.
"""
nPar = 1
nSoundings = len(dataVals)
self.createParametrization(nSoundings, nLayers = nLay, nPar = nPar)
if self._jac is not None:
self._jac.clear()
else:
self._jac = pg.matrix.BlockMatrix()
self.fops1D = []
nData = 0
for i in range(nSoundings):
kwargs = {}
for key, val in self._fopKwargs.items():
if hasattr(val, '__iter__'):
if len(val) == nSoundings:
kwargs[key] = val[i]
else:
kwargs[key] = [val]
else:
kwargs[key] = val
f = None
if issubclass(self._fopTemplate, pg.frameworks.Modelling):
f = self._fopTemplate(**kwargs)
else:
f = type(self._fopTemplate)(self.verbose, **kwargs)
f.setMultiThreadJacobian(self._parPerSounding)
self._fops1D.append(f)
nID = self._jac.addMatrix(f.jacobian())
#print(nID)
self._jac.addMatrixEntry(nID, nData, self._parPerSounding * i)
#print(self._jac)
nData += len(dataVals[i])
self._jac.recalcMatrixSize()
#print("Jacobian size:", self._jac.rows(), self._jac.cols(), nData)
self.setJacobian(self._jac)
def createJacobian(self, par):
"""Create Jacobian matrix by creating individual Jacobians."""
mods = np.asarray(par).reshape(self._nSoundings, self._parPerSounding)
for i in range(self._nSoundings):
self._fops1D[i].createJacobian(mods[i])
#print('Jacobian shape:',np.shape(self._jac))
#return self._jac
def response(self, par):
"""Cut together forward responses of all soundings."""
mods = np.asarray(par).reshape(self._nSoundings, self._parPerSounding)
resp = pg.Vector(0)
for i in range(self._nSoundings):
r = self._fops1D[i].response(mods[i])
resp = pg.cat(resp, r)
return resp
def constraint_matrix(self, dataVals, nLay=None):
""" Create constraint matrix
Parameters
----------
dataVals : list of 1D pyGIMLi data vectors (nSounding x Data per sounding) [list]
nLayers : number of layers (for blocky inversion) [int]
"""
nSoundings = len(dataVals) # number of soundings
boundaries_thk = (nLay - 1) * (nSoundings - 1) # inner mesh boundaries for thicknesses
boundaries_sig = nLay * (nSoundings - 1) # inner mesh boundaries for conductivities
CM = np.zeros((boundaries_thk + boundaries_sig, (nLay *2 - 1) * nSoundings))
h = -np.eye(1, nLay * 2) + np.eye(1, nLay * 2, k = (nLay * 2 - 1))
for i in range(boundaries_thk + boundaries_sig):
CM[i, i:h.shape[1]+i] = h
#print('Size of constraint matrix:', np.shape(CM))
self.CM = pg.utils.toSparseMatrix(CM) # convert to sparse pg matrix
return self.CM
def createWeight(self, dataVals, cWeight_1, cWeight_2, nLay=None):
""" Create constraint weights (cWeights)
Blocky model : vertical constraint weights for both model parameter regions
Parameters
----------
dataVals : list of 1D pyGIMLi data vectors [list]
cWeight_1 : vertical constraint weight (smooth model) or thickness constraint weight
(blocky model) [float]
cWeight_2 : horizontal constraint weight (smooth model) or resistivity constraint weight
(blocky model) [float]
nLay : number of layers (for blocky inversion) [int]
"""
nSoundings = len(dataVals) # number of soundings
""" constraint weights for blocky model """
cWeight_thk = cWeight_1
cWeight_sig = cWeight_2
boundaries_thk = (nLay - 1) * (nSoundings - 1)
boundaries_sig = nLay * (nSoundings - 1)
cWeight_thk = pg.Vector(boundaries_thk, cWeight_thk)
cWeight_sig = pg.Vector(boundaries_sig, cWeight_sig)
cWeight = pg.cat(cWeight_thk, cWeight_sig)
self.cWeight = np.reshape(cWeight, (len(cWeight), 1))
# print('Constraint cWeight:', np.shape(self.cWeight))
#return self.cWeight
def createConstraints(self):
""" create weighted constraint matrix """
self._CW = pg.matrix.LMultRMatrix(self.CM, self.cWeight, verbose = True)
self.setConstraints(self._CW)
print('CW: ', self._CW)
def drawModel(self, ax, model, **kwargs):
"""Draw models as stitched 1D model section."""
mods = np.asarray(model).reshape(self._nSoundings,
self._parPerSounding)
pg.viewer.mpl.showStitchedModels(mods, ax=ax, useMesh=True,
x=self.soundingPos,
**kwargs)
# Create a test model
pos = 10 # number of soundings
nLayers = 3
# Electrical conductivities and thicknesses
x = np.linspace(0,pos,pos, endpoint=False)
sigm_1 = 100/1000
sigm_2 = 10/1000
sigm_3 = 100/1000
thk_1 = np.sin(np.pi/2*x/3) + 2
thk_2 = 1/7*x + 1.5
models = np.zeros((pos, nLayers*2-1))
models[:,0] = thk_1
models[:,1] = thk_2
models[:,2] = sigm_1
models[:,3] = sigm_2
models[:,4] = sigm_3
# Create test data
data_t = []
for p in range(pos):
data_t.append(FDEM1D(models[p,2:], models[p,:2]))
# Initialize Lateral Constraint Modelling class
LC = LCModelling(FDEM1DModelling)
LC.initJacobian(dataVals=data_t, nLay=3)
LC.createJacobian(models)
LC.constraint_matrix(data_t, nLay=3)
LC.createWeight(data_t, cWeight_1=1, cWeight_2=1, nLay=3)
LC.createConstraints()
LC.region(1).setConstraintType(1)
LC.region(2).setConstraintType(1)
# Create an Initial model
m0 = models.copy()
m0[:,:2] = 2
m0[:,2:] = 0.02
# Define inversion
inv = pg.Inversion(LC)
data_true = np.array(data_t).ravel()
relativeError = np.ones_like(data_true)*1e-3
# Run inversion
models_est = inv.run(data_true, relativeError, startModel=m0.ravel(), verbose=True)
and the inversion runs but it shows this prompt
22/02/24 - 13:33:01 - pyGIMLi - INFO - Starting model set from given array. [2. 2. 0.02 0.02 0.02 2. 2. 0.02 0.02 0.02 2. 2. 0.02 0.02
0.02 2. 2. 0.02 0.02 0.02 2. 2. 0.02 0.02 0.02 2. 2. 0.02
0.02 0.02 2. 2. 0.02 0.02 0.02 2. 2. 0.02 0.02 0.02 2. 2.
0.02 0.02 0.02 2. 2. 0.02 0.02 0.02]
22/02/24 - 13:33:01 - pyGIMLi - INFO - Starting inversion.
22/02/24 - 13:33:01 - Core - ERROR - no cWeights defined. You should create constraints matrix first.
min/max(dweight) = 17362.6/8.65792e+07
fop: <__main__.LCModelling object at 0x7fc440ef1860>
Data transformation: <pgcore._pygimli_.RTrans object at 0x7fc440c39940>
Model transformation: <pgcore._pygimli_.RTransLog object at 0x7fc440ef18b0>
min/max (data): 1.2e-05/0.06
min/max (error): 0.1%/0.1%
min/max (start model): 0.02/2
--------------------------------------------------------------------------------
min/max(dweight) = 17362.6/8.65792e+07
found valid constraints matrix. omit rebuild
It indicates that not constraint matrix is defined
22/02/24 - 13:33:01 - Core - ERROR - no cWeights defined. You should create constraints matrix first.
however it later points to
found valid constraints matrix. omit rebuild
I'm not sure if the inversion is using the constraints defined in the LCModelling class in LCModelling.createConstraints().
What does the warning message means?
Thanks!
Thank you @ahsan435 for helping with this code. I just didn't have have the time to dive deep in but I am having the same result with the newest pyGIMLi v1.5.0. So it seems everything is correct. The inversion class LCInversion is actually not needed as all the magic is (as always in pyGIMLi) in the forward operator that holds the Jacobian matrix and the regularization matrix.
@mariacarrizo you were not sure whether the constraint matrices are correct but pg.show(LC.constraints().A) shows that there is smoothness between the neighboring model conductivities and thicknesses.
Nevertheless, the framework is quite bulky and needs to be redesigned based on existing methods and documented by an example.
Dear Sir @halbmy Thank you very much for verifying the results! I'm also working on the LCI method for Time Domain Electromagnetics (TDEM) data. That's why it is also very interesting for me, and I am learning a lot from that code. I converted that code into the time domain but am still not getting satisfying results.