Source code for pygimli.frameworks.timelapse

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import pygimli as pg
from .modelling import MeshModelling

[docs] class MultiFrameModelling(MeshModelling): """Full frame (multiple fop parallel) forward modelling."""
[docs] def __init__(self, modellingOperator, scalef=1.0, **ini): """Init class and jacobian matrix.""" super().__init__() self.ini = ini self.modellingOperator = modellingOperator self.jac = pg.matrix.BlockMatrix() self.scalef = scalef self.fops = [] = 0 self.nm = 0
[docs] def setData(self, data, modellingOperator=None, **ini): """Distribute the data containers amongst the fops.""" modellingOperator = modellingOperator or self.modellingOperator ini = self.ini or ini self.fops = [] for datai in data: fopi = modellingOperator(**ini) fopi.setData(datai) # fopi.setData(pg.Vector(datai)) self.fops.append(fopi)
[docs] def setMeshPost(self, mesh): """Set mesh to all forward operators.""" for fop in self.fops: fop.setMesh(mesh, ignoreRegionManager=True) self.prepareJacobian()
[docs] def setDefaultBackground(self): """Set the default background behaviour.""" regionIds = self.regionManager().regionIdxs()"Found {} regions.".format(len(regionIds))) if len(regionIds) > 1: bk = pg.sort(regionIds)[0]"Region with smallest marker ({0}) " "set to background".format(bk)) self.setRegionProperties(bk, background=True)
@property def parameterCount(self): """Return number of parameters.""" pc = self.regionManager().parameterCount() if pc == 0: pc = self.fops[0].parameterCount return pc * len(self.fops)
[docs] def prepareJacobian(self): """Build up Jacobian block matrix (once the sizes are known).""" self.jac.clear() self.nm = self.regionManager().parameterCount() if self.nm == 0: self.nm = self.fops[0].parameterCount = len(self.fops) print(self.nm, "model cells") nd = 0 for i, fop in enumerate(self.fops): self.jac.addMatrix(fop.jacobian(), nd, i*self.nm) nd += self.jac.recalcMatrixSize() self.setJacobian(self.jac)
[docs] def createConstraints(self, C=None): """Create constraint matrix (special type for this).""" if C is not None: self.C1 = C elif isinstance(super().createConstraints(), pg.SparseMapMatrix): self.C1 = pg.SparseMapMatrix(self.constraintsRef()) # make a copy because it will be overwritten else: self.C1 = self.constraints() self.C = pg.matrix.FrameConstraintMatrix(self.C1, len(self.fops), self.scalef) self.setConstraints(self.C) # cw = self.regionManager().constraintWeights() # self.regionManager().setConstraintsWeights(np.tile(cw, # switch off automagic inside core.inversion which checks for # local modeltransform of the regionManager self.regionManager().setLocalTransFlag(False)
[docs] def response(self, model): """Forward response.""" mod = np.reshape(model, [len(self.fops), -1]) return np.concatenate([fop.response(mo) for fop, mo in zip(self.fops, mod)])
[docs] def createJacobian(self, model): """Create Jacobian matrix.""" mod = np.reshape(model, [len(self.fops), -1]) for i, fop in enumerate(self.fops): fop.createJacobian(mod[i])
[docs] def createDefaultStartModel(self): # , dataVals): """Create standard starting model.""" return pg.Vector(self.nm*, 10.0) # look up in fop
[docs] def createStartModel(self, dataVals): """Create a starting model from mean data values.""" # return np.concatenate([f.createStartModel() for f in self.fops]) self.nm = self.regionManager().parameterCount() return pg.Vector(self.nm*len(self.fops), np.median(dataVals))