Source code for pygimli.frameworks.modelling

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""pyGIMLi - Inversion Frameworks.

These are basic modelling proxies.
"""
import numpy as np
import pygimli as pg

DEFAULT_STYLES = {
    'Default': {
        'color': 'C0',
        'lw': 1.5,
        'linestyle': '-'
    },
    'Data': {
        'color': 'C0',  # blueish
        'lw': 1.5,
        'linestyle': ':',
        'marker': 'o',
        'ms': 6
    },
    'Response': {
        'color': 'C0',  # blueish
        'lw': 2.0,
        'linestyle': '-',
        'marker': 'None',
        'alpha': 0.4
    },
    'Error': {
        'color': 'C3',  # reddish
        'lw': 0,
        'linestyle': '-',
        'elinewidth': 2,
        'alpha': 0.5
    },
}


[docs] class Modelling(pg.core.ModellingBase): """Abstract Forward Operator. Abstract Forward Operator that is or can use a Modelling instance. Can be seen as some kind of proxy Forward Operator. Todo ---- * Docu: - describe members (model transformation, dictionary of region properties) * think about splitting all mesh related into MeshModelling * clarify difference: setData(array|DC), setDataContainer(DC), setDataValues(array) * clarify dataSpace(comp. ModelSpace): The unique spatial or temporal origin of a datapoint (time, coordinates, receiver/transmitter indices, counter) - Every inversion needs, dataValues and dataSpace - DataContainer contain, dataValues and dataSpace - initialize both with initDataSpace(), initModelSpace * createJacobian should also return J """
[docs] def __init__(self, **kwargs): """Initialize. Attributes ---------- fop : pg.frameworks.Modelling data : pg.DataContainer modelTrans : [pg.trans.TransLog()] Parameters ---------- **kwargs : fop : Modelling """ self._fop = None # pg.frameworks.Modelling .. not needed .. remove it self._data = None # dataContainer self._modelTrans = None self.fop = kwargs.pop('fop', None) # super(Modelling, self).__init__(**kwargs) super().__init__(**kwargs) self._regionProperties = {} self._interRegionCouplings = [] self._regionsNeedUpdate = False self._regionChanged = True self._regionManagerInUse = False self.modelTrans = pg.trans.TransLog() # Model transformation operator
def __hash__(self): """Create a hash for Method Manager.""" # ^ pg.utils.dirHash(self._regionProperties) if self._data is not None: return pg.utils.strHash(str(type(self))) ^ hash(self._data) else: return pg.utils.strHash(str(type(self))) def __call__(self, *args, **kwargs): """Call forward operator.""" return self.response(*args, **kwargs) @property def fop(self): """Forward operator.""" return self._fop @fop.setter def fop(self, fop): """Set forward operator.""" if fop is not None: if not isinstance(fop, pg.frameworks.Modelling): pg.critical('Forward operator needs to be an instance of ' 'pg.modelling.Modelling but is of type:', fop) self._fop = fop @property def data(self): """Return data.""" return self._data @data.setter def data(self, d): """Set data (short property setter).""" self.setData(d)
[docs] def setData(self, data): """Set data (actual version).""" if isinstance(data, pg.DataContainer): self.setDataContainer(data) else: self._data = data
[docs] def setDataPost(self, data): """Called when the dataContainer has been set sucessfully.""" pass
[docs] def setDataContainer(self, data): """Set Data container.""" if self.fop is not None: pg.critical('in use?') self.fop.setData(data) else: super().setData(data) self._data = data self.setDataPost(self.data)
@property def modelTrans(self): """Return model transformation.""" self._applyRegionProperties() if self.regionManager().haveLocalTrans(): return self.regionManager().transModel() return self._modelTrans @modelTrans.setter def modelTrans(self, tm): """Set model transformation.""" if isinstance(tm, str): if tm.lower() == "log": tm = pg.trans.TransLog() elif tm.lower() == "linear" or tm.lower() == "lin": tm = pg.trans.Trans() else: # something like "10-1000" raise Exception("Could not use transformation" + tm) self._modelTrans = tm
[docs] def regionManager(self): """Region manager.""" self._regionManagerInUse = True # initialize RM if necessary super().regionManager() # set all local properties self._applyRegionProperties() return super().regionManager()
@property def parameterCount(self): """Return parameter count.""" pC = self.regionManager().parameterCount() if pC == 0: pg.warn("Parameter count is 0") return pC
[docs] def ensureContent(self): """Whatever this is.""" pass
[docs] def initModelSpace(self, **kwargs): """TODO.""" pass
[docs] def createDefaultStartModel(self, dataVals): """Create the default startmodel as the median of the data values.""" pg.critical("'don't use me")
[docs] def createStartModel(self, dataVals=None): """Create the default startmodel as the median of the data values. Overwriting might be a good idea. Its used by inversion to create a valid startmodel if there are no starting values from the regions. """ if dataVals is not None: mv = pg.math.median(dataVals) pg.info("Use median(data values)={0}".format(mv)) sm = pg.Vector(self.parameterCount, mv) else: sm = self.regionManager().createStartModel() return sm
[docs] def clearRegionProperties(self): """Clear all region parameter.""" self._regionChanged = True self._regionProperties = {}
[docs] def regionProperties(self, regionNr=None): """Return dictionary of all properties for region number regionNr.""" if regionNr is None: return self._regionProperties try: return self._regionProperties[regionNr] except KeyError: print(self._regionProperties) pg.error("no region for region #:", regionNr)
[docs] def setRegionProperties(self, regionNr, **kwargs): """Set region properties. regionNr can be '*' for all regions. startModel=None, limits=None, trans=None, cType=None, zWeight=None, modelControl=None, background=None, fix=None, single=None, correlationLengths=None, dip=None, strike=None Parameters ---------- regionNr : int, [ints], '*' Region number, list of numbers, or wildcard "*" for all. startModel : float starting model value limits : [float, float] lower and upper limit for value using a barrier transform trans : str transformation for model barrier: "log", "cot", "lin" cType : int constraint (regularization) type zWeight : float relative weight for vertical boundaries background : bool exclude region from inversion completely (prolongation) fix : float exclude region from inversion completely (fix to value) single : bool reduce region to one unknown correlationLengths : [floats] correlation lengths for geostatistical inversion (x', y', z') dip : float [0] angle between x and x' (first correlation length) strike : float [0] angle between y and y' (second correlation length) """ if regionNr == '*': for regionNr in self.regionManager().regionIdxs(): self.setRegionProperties(regionNr, **kwargs) return elif isinstance(regionNr, (list, tuple)): for r in regionNr: self.setRegionProperties(r, **kwargs) return pg.verbose(f'Set property for region: {regionNr}: {kwargs}') if regionNr not in self._regionProperties: self._regionProperties[regionNr] = {'startModel': None, 'modelControl': 1.0, 'zWeight': 1.0, 'cType': None, # RM defaults 'limits': [0, 0], 'trans': 'Log', # RM defauts 'background': None, 'single': None, 'fix': None, 'correlationLengths': None, 'dip': None, 'strike': None, } for key in list(kwargs.keys()): val = kwargs.pop(key) if val is not None: if self._regionProperties[regionNr][key] != val: self._regionsNeedUpdate = True self._regionProperties[regionNr][key] = val if len(kwargs) > 0: pg.warn('Unhandled region properties:', kwargs)
[docs] def setInterRegionCoupling(self, region1, region2, weight=1.0): """Set the weighting for constraints across regions.""" if region1 == "*": region1 = self.regionManager().regionIdxs() else: region1 = [region1] if region2 == "*": region2 = self.regionManager().regionIdxs() else: region2 = [region2] for reg1 in region1: for reg2 in region2: if reg1 != reg2 and \ (not self._regionProperties[reg1]['background'] and not self._regionProperties[reg2]['background']): self._interRegionCouplings.append([reg1, reg2, weight]) self._regionsNeedUpdate = True
def _applyRegionProperties(self): """Apply the region properties from dictionary into the region man.""" if self._regionsNeedUpdate is False: return # call super class her because self.regionManager() calls allways # __applyRegionProperies itself rMgr = super().regionManager() for rID, vals in self._regionProperties.items(): if vals['fix'] is not None: if rMgr.region(rID).fixValue() != vals['fix']: vals['background'] = True rMgr.region(rID).setFixValue(vals['fix']) self._regionChanged = True if vals['background'] is not None: if rMgr.region(rID).isBackground() != vals['background']: rMgr.region(rID).setBackground(vals['background']) self._regionChanged = True if vals['single'] is not None: if rMgr.region(rID).isSingle() != vals['single']: rMgr.region(rID).setSingle(vals['single']) self._regionChanged = True if vals['startModel'] is not None: rMgr.region(rID).setStartModel(vals['startModel']) if vals['trans'] is not None: rMgr.region(rID).setModelTransStr_(vals['trans']) if vals['cType'] is not None: if rMgr.region(rID).constraintType() != vals['cType']: self.clearConstraints() rMgr.region(rID).setConstraintType(vals['cType']) if vals['zWeight'] is not None: rMgr.region(rID).setZWeight(vals['zWeight']) self.clearConstraints() rMgr.region(rID).setModelControl(vals['modelControl']) if vals['limits'][0] != 0: rMgr.region(rID).setLowerBound(vals['limits'][0]) if vals['limits'][1] > vals['limits'][0]: rMgr.region(rID).setUpperBound(vals['limits'][1]) if vals['correlationLengths'] is not None: self.clearConstraints() if vals['dip'] is not None: self.clearConstraints() if vals['strike'] is not None: self.clearConstraints() for r1, r2, w in self._interRegionCouplings: rMgr.setInterRegionConstraint(r1, r2, w) self._regionsNeedUpdate = False
[docs] def setDataSpace(self, **kwargs): """Set data space, e.g., DataContainer, times, coordinates.""" if self.fop is not None: pg.critical('in use?') self.fop.setDataSpace(**kwargs) else: data = kwargs.pop('dataContainer', None) if isinstance(data, pg.DataContainer): self.setDataContainer(data) else: print(data) pg.critical("nothing known to do? " "Implement me in derived classes")
[docs] def estimateError(self, data, **kwargs): """Create data error fallback when the data error is not known. Should be implemented method-specific. """ raise Exception("Needed?? Implement me in derived classes")
# data = data * (pg.randn(len(data)) * errPerc / 100. + 1.) # return data
[docs] def drawModel(self, ax, model, **kwargs): """Draw a model into a given axis.""" if self.fop is not None: pg.critical('in use?') self.fop.drawModel(ax, model, **kwargs) else: print(kwargs) raise Exception("No yet implemented")
[docs] def drawData(self, ax, data, **kwargs): """Draw data into a given axis.""" if self.fop is not None: self.fop.drawData(ax, data, **kwargs) else: print(kwargs) raise Exception("No yet implemented")
[docs] class LinearModelling(Modelling): """Modelling class for linearized problems with a given matrix."""
[docs] def __init__(self, A): """Initialize by storing the (reference to the) matrix.""" super().__init__() self.A = A self.setJacobian(self.A)
[docs] def response(self, model): """Linearized forward modelling by matrix-vector product.""" return self.A * model
[docs] def createJacobian(self, model): """Do not compute a jacobian (linear).""" pass
@property def parameterCount(self): """Define the number of parameters from the matrix size.""" return self.A.cols()
[docs] class Block1DModelling(Modelling): """General forward operator for 1D layered models. Model space: [thickness_i, parameter_jk], with i = 0 - nLayers-1, j = (0 .. nLayers), k=(0 .. nPara) """
[docs] def __init__(self, nPara=1, nLayers=4, **kwargs): """Constructor. Parameters ---------- nLayers : int [4] Number of layers. nPara : int [1] Number of parameters per layer (e.g. nPara=2 for resistivity and phase) """ self._nLayers = 0 super(Block1DModelling, self).__init__(**kwargs) self._withMultiThread = True self._nPara = nPara # number of parameters per layer self.initModelSpace(nLayers)
@property def nPara(self): """Number of parameters.""" return self._nPara @property def nLayers(self): """Number of layers.""" return self._nLayers @nLayers.setter def nLayers(self, nLayers): """Set number of layers.""" return self.initModelSpace(nLayers)
[docs] def initModelSpace(self, nLayers): """Set number of layers for the 1D block model.""" if nLayers == self._nLayers: return self._nLayers = nLayers if nLayers < 2: pg.critical("Number of layers need to be at least 2") mesh = pg.meshtools.createMesh1DBlock(nLayers, self._nPara) self.clearRegionProperties() self.setMesh(mesh) # setting region 0 (layers) and 1..nPara (values) for i in range(1 + self._nPara): self.setRegionProperties(i, trans='log') if self._withMultiThread: self.setMultiThreadJacobian(2*nLayers - 1)
# self._applyRegionProperties()
[docs] def drawModel(self, ax, model, **kwargs): """Draw model into a given axis.""" pg.viewer.mpl.drawModel1D(ax=ax, model=model, plot='loglog', xlabel=kwargs.pop('xlabel', 'Model parameter'), **kwargs) return ax
[docs] def drawData(self, ax, data, err=None, label=None, **kwargs): """Default data view. Modelling creates the data and should know best how to draw them. Probably ugly and you should overwrite it in your derived forward operator. """ nData = len(data) yVals = range(1, nData+1) ax.loglog(data, yVals, label=label, **DEFAULT_STYLES.get(label, DEFAULT_STYLES['Default']) ) if err is not None: ax.errorbar(data, yVals, xerr=err*data, label='Error', **DEFAULT_STYLES.get('Error', DEFAULT_STYLES['Default']) ) ax.set_ylim(max(yVals), min(yVals)) ax.set_xlabel('Data') ax.set_ylabel('Data Number') return ax
[docs] class MeshModelling(Modelling): """Modelling class with a mesh discretization."""
[docs] def __init__(self, **kwargs): super().__init__(**kwargs) self._axs = None self._meshNeedsUpdate = True self._baseMesh = None # optional p2 refinement for forward task self._refineP2 = False self._refineH2 = True self._pd = None self._C = None # custom Constraints matrix
def __hash__(self): """Unique hash for caching.""" return super().__hash__() ^ hash(self.mesh()) @property def mesh(self): """Return mesh.""" pg._r("inuse ?") if self._fop is not None: pg._r("inuse ?") return self._fop.mesh else: return self.mesh() @property def paraDomain(self): """Return parameter (inverse) mesh.""" # We need our own copy here because its possible that we want to use # the mesh after the fop was deleted if not self.mesh(): pg.critical('paraDomain needs a mesh') self._pd = pg.Mesh(self.regionManager().paraDomain()) return self._pd
[docs] def setCustomConstraints(self, C): """ Set custom constraints matrix for lazy evaluation. To remove them set it to 'None' again. """ self._C = C
[docs] def createConstraints(self): """Create constraint matrix.""" # just ensure there is valid mesh # self.mesh() foundGeoStat = False for reg, props in self.regionProperties().items(): if not props['background'] and \ props['correlationLengths'] is not None or \ props['dip'] is not None or props['strike'] is not None: cL = props.get('correlationLengths') or 5 dip = props.get('dip') or 0 strike = props.get('strike') or 0 pg.info('Creating GeostatisticConstraintsMatrix for region' + f' {reg} with: I={cL}, dip={dip}, strike={strike}') if foundGeoStat is True: pg.critical('Only one global GeostatisticConstraintsMatrix' 'possible at the moment.') # keep a copy of C until refcounting in the core works self._C = pg.matrix.GeostatisticConstraintsMatrix( mesh=self.paraDomain, I=cL, dip=dip, strike=strike, ) foundGeoStat = True self.setConstraints(self._C) if foundGeoStat is False: super().createConstraints() return self.constraints()
[docs] def paraModel(self, model): """Return parameter model, i.e. model mapped back with cell markers.""" mod = model[self.paraDomain.cellMarkers()] if isinstance(mod, np.ndarray): mod = pg.Vector(mod) # Else next line fails as np.array does not allow set attributes. mod.isParaModel = True return mod
[docs] def ensureContent(self): """Internal function to ensure there is a valid initialized mesh. Initialization means the cell marker are recounted and/or there was a mesh refinement or boundary enlargement, all to fit the needs for the method-depending forward problem. """ # Need to call this once to be sure the mesh is initialized when needed self.mesh()
[docs] def setMeshPost(self, mesh): """Interface to be called when the mesh has been set successfully. Might be overwritten by child classes. """ pass
[docs] def createRefinedFwdMesh(self, mesh): """Refine the current mesh for higher accuracy. This is called automatic when accessing self.mesh() so it ensures any effect of changing region properties (background, single). """ m = pg.Mesh(mesh) if self._refineH2: pg.info("Creating refined mesh (H2) to solve forward task.") m = m.createH2() if self._refineP2: pg.info("Creating refined mesh (P2) to solve forward task.") m = m.createP2() pg.info("Mesh for forward task:", m) return m
[docs] def createFwdMesh_(self): """Create forward mesh.""" pg.info("Creating forward mesh from region infos.") m = pg.Mesh(self.regionManager().mesh()) regionIds = self.regionManager().regionIdxs() for iId in regionIds: pg.verbose("\tRegion: {0}, Parameter: {1}, PD: {2}," " Single: {3}, Background: {4}, Fixed: {5}".format( iId, self.regionManager().region(iId).parameterCount(), self.regionManager().region(iId).isInParaDomain(), self.regionManager().region(iId).isSingle(), self.regionManager().region(iId).isBackground(), self.regionManager().region(iId).fixValue())) m = self.createRefinedFwdMesh(m) self.setMeshPost(m) self._regionChanged = False super().setMesh(m, ignoreRegionManager=True) if self._C is not None: pg.info('Set custom constraints matrix.') self.setConstraints(self._C)
[docs] def mesh(self): """Returns the currently used mesh.""" self._applyRegionProperties() if self._regionManagerInUse and self._regionChanged is True: self.createFwdMesh_() return super().mesh()
[docs] def setMesh(self, mesh, ignoreRegionManager=False): """Set mesh and specify whether region manager can be ignored.""" # keep a copy, just in case self._baseMesh = mesh if ignoreRegionManager is False: self._regionManagerInUse = True # Modelling without region manager if ignoreRegionManager is True or not self._regionManagerInUse: self._regionManagerInUse = False if self.fop is not None: pg.critical('in use?') self.fop.setMesh(mesh, ignoreRegionManager=True) else: super(Modelling, self).setMesh(mesh, ignoreRegionManager=True) pass self.setMeshPost(mesh) return # copy the mesh to the region manager who renumber cell markers self.clearRegionProperties() self.regionManager().setMesh(mesh) self.setDefaultBackground()
[docs] def setDefaultBackground(self): """Set the lowest region to background if several exist.""" regionIds = self.regionManager().regionIdxs() pg.info("Found {} regions.".format(len(regionIds))) if len(regionIds) > 1: bk = pg.sort(regionIds)[0] pg.info("Region with smallest marker set to background " "(marker={0})".format(bk)) self.setRegionProperties(bk, background=True)
[docs] def drawModel(self, ax, model, **kwargs): """Draw the model as mesh-based distribution.""" mod = None # TODO needs to be checked if mapping is always ok (region example) # is (len(model) == self.paraDomain.cellCount() or \ if hasattr(model, "isParaModel") and model.isParaModel is False: # pg._y(model.isParaModel) mod = self.paraModel(model) elif hasattr(model, "isParaModel") and model.isParaModel is True: # pg._g(model.isParaModel) mod = model elif len(model) == self.paraDomain.nodeCount(): # why nodeCount? a field as model result? # pg._b('node count') mod = model elif len(model) == self.paraDomain.cellCount(): # pg._b('cell count') mod = model else: mod = self.paraModel(model) if ax is None: if self._axs is None: self._axs, _ = pg.show() ax = self._axs if hasattr(ax, '__cBar__'): # we assume the axes already holds a valid mappable and we only # update the model data cBar = ax.__cBar__ kwargs.pop('label', None) kwargs.pop('cMap', None) pg.viewer.mpl.setMappableData(cBar.mappable, mod, **kwargs) else: diam = kwargs.pop('diam', None) ax, cBar = pg.show(mesh=self.paraDomain, data=mod, label=kwargs.pop('label', 'Model parameter'), logScale=kwargs.pop('logScale', False), ax=ax, **kwargs ) if diam is not None: pg.viewer.mpl.drawSensors(ax, self.data.sensors(), diam=diam, edgecolor='black', facecolor='white') return ax, cBar
[docs] class PetroModelling(MeshModelling): """Combine petrophysical relation with the modelling class f(p). Combine petrophysical relation :math:`p(m)` with a modelling class :math:`f(p)` to invert for the petrophysical model :math:`p` instead of the geophysical model :math:`m`. :math:`p` be the petrophysical model, e.g., porosity, saturation, ... :math:`m` be the geophysical model, e.g., slowness, resistivity, ... """
[docs] def __init__(self, fop, petro, **kwargs): """Save forward class and transformation, create Jacobian matrix.""" self._f = fop # self._f createStartModel might be called and depends on the regionMgr self._f.regionManager = self.regionManager # self.createRefinedFwdMesh depends on refinement strategy of self._f self.createRefinedFwdMesh = self._f.createRefinedFwdMesh super(PetroModelling, self).__init__(fop=None, **kwargs) # petroTrans.fwd(): p(m), petroTrans.inv(): m(p) self._petroTrans = petro # class defining p(m) self._jac = pg.matrix.MultRightMatrix(self._f.jacobian()) self.setJacobian(self._jac)
@property def petro(self): """Petrophysical model transformation.""" return self._petroTrans
[docs] def setMeshPost(self, mesh): """Set mesh after init.""" self._f.setMesh(mesh, ignoreRegionManager=True)
[docs] def setDataPost(self, data): """Set data after init.""" self._f.setData(data)
[docs] def createStartModel(self, data): """Use inverse transformation to get m(p) for the starting model.""" sm = self._f.createStartModel(data) pModel = self._petroTrans.inv(sm) return pModel
[docs] def response(self, model): """Use transformation to get p(m) and compute response f(p).""" tModel = self._petroTrans.fwd(model) ret = self._f.response(tModel) return ret
[docs] def createJacobian(self, model): r"""Fill the individual jacobian matrices. J = dF(m) / dm = dF(m) / dp * dp / dm """ tModel = self._petroTrans.fwd(model) self._f.createJacobian(tModel) self._jac.A = self._f.jacobian() self._jac.r = self._petroTrans.deriv(model) # set inner derivative # print(self._jac.A.rows(), self._jac.A.cols()) # print(self._jac.r) # pg._r("create Jacobian", self, self._jac) self.setJacobian(self._jac) # to be sure .. test if necessary
# 220817 to be changed later!! # class JointModelling(Modelling):
[docs] class JointModelling(MeshModelling): """Cumulative (joint) forward operator."""
[docs] def __init__(self, fopList): """Initialize with lists of forward operators.""" super().__init__() self.fops = fopList self.jac = pg.matrix.BlockMatrix() # self.modelTrans = self.fops[0].modelTrans self.modelTrans = pg.trans.TransLogLU() self.fops[0].regionManager() self.setRegionManager(self.fops[0].regionManagerRef())
[docs] def createStartModel(self, data): """Use inverse transformation to get m(p) for the starting model.""" sm = self._f.createStartModel(data) pModel = self._petroTrans.inv(sm) return pModel
[docs] def response(self, model): """Concatenate responses for all fops.""" resp = [] for f in self.fops: resp.extend(f.response(model)) return resp
[docs] def createJacobian(self, model): """Fill the individual Jacobian matrices.""" self.initJacobian() for f in self.fops: f.createJacobian(model)
[docs] def setData(self, data): """Distribute list of data to the forward operators.""" if len(data) != len(self.fops): pg.critical("Please provide data for all forward operators") self._data = data nData = 0 for i, fi in enumerate(self.fops): fi.setData(data[i]) self.jac.addMatrix(fi.jacobian(), nData, 0) nData += data[i].size() # update total vector length self.setJacobian(self.jac)
[docs] def setMesh(self, mesh, **kwargs): # to be removed from here """Set the parameter mesh to all fops.""" for fi in self.fops: fi.setMesh(mesh)
# 220817 to be implemented!! # class JointMeshModelling(JointModelling): # def __init__(self, fopList): # super().__init__(self, fopList) # self.setRegionManager(self.fops[0].regionManagerRef())
[docs] class LCModelling(Modelling): """2D Laterally constrained (LC) modelling. 2D Laterally constrained (LC) modelling based on BlockMatrices. """
[docs] def __init__(self, fop, **kwargs): """Parameters: fop class .""" super(LCModelling, self).__init__() 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
[docs] 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)
[docs] def initModelSpace(self, nLayers): """Initialize model space.""" for i, f in enumerate(self._fops1D): f.initModelSpace(nLayers)
[docs] 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
[docs] 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]) # print("i:", i, mods[i], r) resp = pg.cat(resp, r) return resp
[docs] 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])
[docs] def createParametrization(self, nSoundings, nLayers=4, nPar=1): """Create LCI mesh and suitable constraints informations. Parameters ---------- nLayers : int Numbers of depth layers nSoundings : int Numbers of 1D measurements to laterally constrain nPar : int Numbers of independent parameter types, e.g., nPar = 1 for VES (invert for resisitivies), nPar = 2 for VESC (invert for resisitivies and phases) """ nCols = (nPar+1) * nLayers - 1 # fail for VES-C self._parPerSounding = nCols self._nSoundings = nSoundings self._mesh = pg.meshtools.createMesh2D(range(nCols + 1), range(nSoundings + 1)) self._mesh.rotate(pg.RVector3(0, 0, -np.pi/2)) cm = np.ones(nCols * nSoundings) * 1 if not self._singleRegion: for i in range(nSoundings): for j in range(nPar): cm[i * self._parPerSounding + (j+1) * nLayers-1: i * self._parPerSounding + (j+2) * nLayers-1] += (j+1) self._mesh.setCellMarkers(cm) self.setMesh(self._mesh) pID = self.regionManager().paraDomain().cellMarkers() cID = [c.id() for c in self._mesh.cells()] # print(np.array(pID)) # print(np.array(cID)) # print(self.parameterCount perm = [0]*self.parameterCount for i in range(len(perm)): perm[pID[i]] = cID[i] # print(perm) self.regionManager().permuteParameterMarker(perm)
# print(self.regionManager().paraDomain().cellMarkers())
[docs] def initJacobian(self, dataVals, nLayers, nPar=None): """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) if nPar is None: # TODO get nPar Infos from fop._fopTemplate nPar = 1 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__'): if len(val) == nSoundings: 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.J.rows(), self.J.cols(), nData) self.setJacobian(self._jac)
[docs] 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)
[docs] class ParameterModelling(Modelling): """Model with symbolic parameter names instead of numbers."""
[docs] def __init__(self, funct=None, **kwargs): """Initialize, optionally with given function.""" self.function = None self._params = {} self.dataSpace = None # x, t freqs, or whatever self.defaultModelTrans = 'lin' super(ParameterModelling, self).__init__(**kwargs) if funct is not None: self._initFunction(funct)
@property def params(self): """Return number of parameters.""" return self._params def _initFunction(self, funct): """Init any function and interpret possible args and kwargs.""" self.function = funct # the first varname is suposed to be f or freqs self.dataSpaceName = funct.__code__.co_varnames[0] pg.debug('data space:', self.dataSpaceName) args = funct.__code__.co_varnames[1:funct.__code__.co_argcount] for varname in args: if varname != 'verbose': pg.debug('add parameter:', varname) self._params[varname] = 0.0 # nPara = len(self._params.keys()) # not used! for i, [k, p] in enumerate(self._params.items()): self.addParameter(k, id=i, cType=0, single=True, trans=self.defaultModelTrans, startModel=1)
[docs] def response(self, params): """Compute and return model response.""" if np.isnan([*params]).any(): print(params) pg.critical('invalid params for response') if self.dataSpace is None: pg.critical('no data space given') ret = self.function(self.dataSpace, *params) return ret
[docs] def setRegionProperties(self, k, **kwargs): """Set Region Properties by parameter name.""" if isinstance(k, int) or (k == '*'): super(ParameterModelling, self).setRegionProperties(k, **kwargs) else: self.setRegionProperties(self._params[k], **kwargs)
[docs] def addParameter(self, name, id=None, **kwargs): """Add a parameter.""" if id is None: id = len(self._params) self._params[name] = id self.regionManager().addRegion(id) self.setRegionProperties(name, **kwargs) return id
[docs] def drawModel(self, ax, model): """Draw model.""" label = '' for k, p in self._params.items(): label += k + "={0} ".format(pg.utils.prettyFloat(model[p])) pg.info("Model: ", label)
[docs] class PriorModelling(MeshModelling): """Forward operator for grabbing values out of a mesh (prior data)."""
[docs] def __init__(self, mesh=None, pos=None, **kwargs): """Init with mesh and some positions that are converted into ids.""" super().__init__(**kwargs) self.pos = pos if mesh is not None: self.setMesh(mesh)
[docs] def setMesh(self, mesh): """Set mesh, save index vector and compute Jacobian.""" super().setMesh(mesh) self.ind = np.zeros(len(self.pos), dtype=np.int32) for i, po in enumerate(self.pos): cell = mesh.findCell(po) if cell is None: raise IndexError(f"Could not find cell at position {po}!") else: self.ind[i] = cell.id() # self.ind = np.array([mesh.findCell(po).id() for po in self.pos]) self.J = pg.SparseMapMatrix() self.J.resize(len(self.ind), mesh.cellCount()) for i, ii in enumerate(self.ind): self.J.setVal(i, ii, 1.0) self.setJacobian(self.J)
[docs] def response(self, model): """Return values at the indexed cells.""" return model[self.ind]
[docs] def createJacobian(self, model): """Do nothing (linear).""" pass
[docs] def createRefinedFwdMesh(self, mesh): """Create refined forward mesh: do nothing here to prevent this.""" return mesh