Error message with fop.modelTrans in LCInversion

[mariacarrizo]

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.1

Operating 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 fop1D is the EM forward function in empymod in 1D
• The second class FDEM1D is the EM forward modelling framework in 1D
• The third class is LCModelling to perform the modelling using FDEM1D
• The third class is LCInversion to 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)

[mariacarrizo]

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!

[ahsan435]

Hello @mariacarrizo I run your code and it is showing these results

[halbmy]

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.

[ahsan435]

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.