Object structure

Study

class mse.study.study(aName, nbParam=0)

getParam(): Return the list of parameter values

restoreState(): Can be overload by user

saveState(): Can be overload by user

setGeometry(geomName, tol=-1, aMeshCallBack=None)

This function builds the spatial discretisation (meshes) of the study zone

Parameters:: string (geomName) – the geomzone of the mse env geometry directory

:param tol double : the number of nodes by lenght unit

setParam(aParam, index=None)

Sets new parameter for the study.

The parameters corresponding to p1,p2,p3,… in the source terms. The index of param p1 is 0, p2 is 1,… pn is n-1. This function updates the list of study parameters with the new aParam parameters in specific index (by default index = [0,1,…,len(aParam)-1) : study.Param[index[i]] = aParam[i] Then calls the setParamDolf function for all the dynamicalSystem in the study.

Parameters:

aParam (list(float)) – list of new parameters to change. The list must be equal length of index. If aParam is longer than index, the rest of the list will be ignored.
index (list(int)/int) – if index is list, index of parameters to modify, starting at 0. Default=list(range(self.nbParam)). If index[i] > number of parameters, prints an error message, stop updating of study.param and gswitchies to updating param of dynamicalSystem.

if index is int, create new list index = [index, index+1, …, index+len(aParam)-1], so that len(index)==len(aParam).

setSimuDir(aSimuDir): This fonction allows to set the simulation traj output It builds the directory if necessary param aSimuDir string: add to the study directory

Diffusion process

Constant Diffusion

class mse.DISPERSIONS.diffusionModel.diffusionModel(aDs, aDiffCoef=1)

enableAdjoint(): Function to check if derivate file to param exist, else create it

setDiffCoef(aDiffCoef)

Set a new const coef for the diffusion

Parameters:: aDiffCoef (float) – new coeffficient for diffusion

Diffusion with parameter

class mse.DISPERSIONS.diffusionModel.txtDiffusionModel(aDs, aTxtFile)

Class used to add diffusion in model with parameters.

The diffusion model is stored in a txt file. NOTE: only avaible for compment in 2d The txt file must be this form: fax=… fay=… with ‘a’ the compement. We are in 2 dim, each compement contain 2 diffusion

Parameters:

aDs (dynamicalSystem) – class allows applying a PDE to a geometry
strModel (str) – file contain diffModel.

enableAdjoint(): Function to check if derivate file to param exist, else create it

setDiffCoef(aTxtFile)

Set a new diffusion model.

The new model is store in a txt file.

Parameters:: aTxtFile (str) – the file wich contain the new diffusion model

Source Term Model

class mse.MODELS.models.txtModel(aDs, aTxtFile)

Class used to add source term in model.

The model is stored in a txt file and applied to a dynamical system.

Parameters:

aDs (dynamicalSystem) – class allows applying a PDE to a geometry

aTxtFile (str) – file which contain source term, file = aTxtFile+”.txt”. An other file with partial derivate can be create in the same directory, file = aTxtFile+”Dv.txt”.

members:

noindex:

covariable

builderCovariableFromRaster

class mse.COVARIABLE.covariable.builderCovariableFromRaster(aName, aDs, pathToRaster)

This class write files containing a covariable from a raster input.

Parameters:

string (pathToRaster) – the name of the covariable
dynamicalSystem (aDs) – the dynamical system concerned
string – path to de raster (”…/rasterT”). pathToRaster and pathToRaster”Info.txt” should exits

build(): output: write file containing covariable

covariable

class mse.COVARIABLE.covariable.covariable(aName, aDS, aSubDir='.')

A covariable can appear in the model of the dynamical system. It can be space and time heterogen.

Parameters:

string (aSubDir) – the name used in the model
dynamicalSystem (aDS) – the dynamical system concerned
string – sub directory containing files defining the values of the covariable

computeMass(): compute the meseaure of the covariable on the geometrical space

load(require=False): set the value of the covariable from file.

setValue(): call load.

covariableFromExpression

class mse.COVARIABLE.covariable.covariableFromExpression(aName, aDS)

Allow to build a covariable from a Expression

setValue(aExpression): call load.

covariableCst

class mse.COVARIABLE.covariable.covariableCst(aName, aDS, pathToTimeValueFile='')

It defines a covariable space homogen. It could be time variable using pathToTimeValueFile.

Parameters:

string (pathToTimeValueFile) – the name of the covariable aeppearing in the model
dynamicalSystem (aDS) – dynamicalSystem concerned by the covariable
string – if time variable, self.pathToTimeValueFile is a numpy save “xxx.npy” containing tow collums (time|value)

axpy(alpha, anOtherCovCst): make self+=alpha*anOtherCovCst

buildTimeVal(aTimeFct, t0, t1, stepT)

build a cov file pathToTimeValueFile with value returned by aTimeFct

Parameters:

function (aTimeFct) – returning value for any time
float (stepT) – init time
float – end time
float – step time

getValue(t)

get the value at time t

Parameters:: float (t) – the time
Retruns:: the value interpolated in file pathToTimeValueFile

setValue(value): set the value of the dolfinx Constant with value

setValueAtTime(t)

set the dolcfinx Constant with value at time t

Parameters:: float (t) – the time

setValueAverage(aTimeBegin, aTimeEnd): set the value using the pathTimeValueFile

covariableFromCovFile

class mse.COVARIABLE.covariable.covariableFromCovFile(aName, aDs, aSubDir='.')

It build a covariable from file(s)

beginTimeBuild(): It build an empty time file in view of build a time-variable covariable.

getFirstTimeInterval(): read the first time

saveStepTime(step, time): Add a value associated to time

setValue(aTime): set the covariable value at time aTime. The values are interpolated at t=aTime and the dolfinx function is updated. This fonction is used by the time scheme(implict, theta …)

setValue(aTime): set the covariable value at time aTime. The values are interpolated at t=aTime and the dolfinx function is updated. This fonction is used by the time scheme(implict, theta …)

setValueAverage(aTimeBegin, aTimeEnd): set the value of the covariable using average value on the time interval

setValueFromExpression(aExpression): set the value from expression

Dynamical System

the module dynamicalSystem allow to apply a model on a geom patch (2d or 1d)

class mse.dynamicalSystem.dynamicalSystem(aStudy, aGeoDim, aMeshFile, nbComps=1, obsProc=None)

addSecondOrderTerm(aModel)

set the second order term of the pde

Parameters:: aModel (txtDiffusionModel / diffusionModel) – the new diffusion model

addSourcesTerm(aModel)

Add a model to the list of source term of the dynamical system

Parameters:: aModel (txtModel) – the source term model to be added

exportStateVtk(nStep)

Export current state to vtk

Parameters:: nStep (int) – num export vtk

mse.dynamicalSystem.computeMass(ds, mass)

It compute the mass of the current state of the dynamicalcSystem

Parameters:

ds (dynamicalSystem)
mass (array) – array[0 to 0+number of compatiment] will be filled

State

A state is the mermory representation of the solution of the PDE at a fixed time. It allows to build the FEM représentation of the simulation.

class mse.state.state(*args)

axpy(aOtherState, alpha): self=self+alpha*otherState

closeFilesToExportVtk(): close the file containing vtk animation

exportVtk(nStep)

build a vtk file with self and add in animation file

Parameters:: nStep (real) – the number of step or time

initFilesToExportVtk(filePrefix)

prepar vtk export in pvd file for animation

:param string filePrefix : prefix name for vtk file

scale(alpha): self=alpha*self

zero(): self=0

observation Processus

class mse.obsProcess.obsProcess(aDS, aDirURef=None, enableSysAdjoint=False)

Virtual class for observation process making the link between the observations and the simulation.: For exemple, an observation process is introduced measuring the likelihood of the data for a set of fixed parameters.

Parameters:

aDS (dynamicalSystem) – the dynamical system involved
aDirURef (string) – the directory contaning the observations
enableSysAdjoint (bool) – a boolean to anable the adjoint system computation

computeGradV()

By default, it builds the adjoint state via systemAdjoint class and build the grad of the obsprocess using \(<P.\partial_aF(a,u_a)e_i>\).

read file from adjoint simu directory and create gradV.

quadratic difference

class mse.obsProcess.quadraticDiff(aDS, aDirURef=None, enableSysAdjoint=False)

Subclass of obsProcess

This class implements the observation process as a quadratic difference.

Parameters:

aDS (dynamicalStructure)
aDirURef (str) – an absolute pasth to the directory where observation files are. By default in dir study.absoulteSutyDir+’UREF/’.

computeTildevPrim(temps)

It computes \(\tilde v\) described in Computing variations using adjoint method

Here, we have : \(\tilde v'(u_k) = [0.5*(u_k-u_{ref_k})]\), with k = nbComps

Parameters:: temps (float) – time of simulation to calculate \(\tilde v'(u_k)(.,t=temps)\)

computeV()

Method to compute the quadratic difference for the current dynamical system parameters. Assume simulation has been running. The returned value is

\[\sum_k \int_Q |u_k-u_{ref_k}|^2.\]

:return

V (float): value of the observation process, the quadratic difference.

discrets observations

class mse.obsProcess.obsDiscret(obsPointsCoord, obsTimePerPoints)

Defines the spatial and temporal discretization of the observations.

Parameters:

obsPos (numpy.ndarray) – 2D array containing the coordinates of the observation points. For example, np.array([[1, 0], [0, 1]]) specifies two observation points located at (1, 0) and (0, 1).
timePerObs (numpy.ndarray) –
Array describing, for each observation time, the indices of the corresponding observation points. For example, np.array([[0.1, [0, 1]], [0.2, [1]]]) means that:
- at t = 0.1, two observations are taken at points with indices 0 and 1;
- at t = 0.2, one observation is taken at point with index 1.
Thus, timePerObs associates each observation time with the indices of the observation points where measurements are performed.

discret observation process

class mse.obsProcess.discretObsProcess(aDS, aDirURef=None, enableSysAdjoint=False, obsPointsCoord=None, obsTimePerPoints=None)

This class provides a generic interface to define a discrete observation process. It must be subclassed to implement a specific observation model by overloading the methods computeVtildePrim() and computeV(), which together are responsible for filling vPrimObs (see for example discretQuadraticObs).

Note about the numerical aspects and implementation: - Because the observation operator must support both interpolation and pointwise assignment on a dolfin function, the implementation uses a Dirac mass located at the degree of freedom (DoF) nearest to each observation point.

For each observation point:

identify the nearest DoF,

compute the normalization factor ensuring that the pointwise value is consistent with the finite-element representation.

The theta method of the adjoint system is enforced to 0 (self.timeSchemeBackward.theta = 0), because observations taken

at the final time of the simulation must be treated as fully explicit.

computeTildevPrim(temps): compute the state tilte v described in Computing variations using adjoint method

discret quadratic observation process

class mse.obsProcess.discretQuadraticObsProcess(aDS, aDirURef=None, enableSysAdjoint=False, obsPointsCoord=None, obsTimePerPoints=None)

Specialized observation process based on discretObsProcess.

This class overrides computeV() to define the quadratic difference with respect to a given reference (quadratic misfit). In particular, it sets tildeVPrimObs, which is then used to assemble and compute the adjoint system.

computeV(): Method to compute the observation process for the current study parameters.

Adjoint system

class mse.systemAdjoint.systemAdjoint(aStudy)

An attribute of study for compute adjoint state

To calculate the gradient of J, we have to solve an adjoint system such as :

\[< \tilde v'(u_\theta)+\partial _u F(\theta,u_\theta)^*P,w >_H=0\]

with \(\tilde v'\) wich depends of the observation processus

and \(\partial_uF(\theta,u_\theta)^*\) such that :

\[<\partial_uF(\theta,u_\theta)\epsilon;w> = <\epsilon;\partial_uF(\theta,u_\theta)^* w>\]

buildAdjointState()

Build adjoint state, It computes the solution of the adjoint system contained in the study and depending in the obsProcess via the obsProcess menber tildev.

A fisrts implementation consists in building the adjoint system, using dolfinx, for a system using 1 compartiments without interactions. A second step is to write if this computation holds for 2 compartiments. A third step is to integrate 2D1D interactions

valJAndGradJ(param)

return val of J(param) and gradJ(param)

The simulation bakward will be run in computeGradV()

Parameters:: array (param) – parameter to evaluate \(J\) and \(\nabla J\)

Explore trajectory

class mse.explorTrajectory.explorTrajectory(aStudy, trajDirName, varName='ds')

Parameters:

aStude (study) – contains the system of dynamical system
trajDirName – the directory containing the states writting in files

interpolToTime(aTime)

Function to interpolate the state in time aTime.

It update the state ds.trajState :param float time: time to interpolate value

replay()

function to update the state ds.utm1 at the next time

use it in a loop

rewind()

prepar loop in trajectory steps to update ds.utm1

This function must be called before self.replay

stopReplay(): Function call in self.replay to stop the rewind

Interactions

Interaction 2D1D

class mse.INTERACTIONS.interaction.interaction2D1D(aStudy, aDS2d, aDS1d, aCompsInvolved=None)

This Class implement a 2D1D interaction

Parameters:

aStudy (study) – the study
aDS2d (dynamicalSystem) – the first dynamical system is the 2d dynamical system
aDS1d (dynamicalSystem) – the second dynamical system is the 1d dynamical system
aCompsInvolved=NONE (numpy.array) – the list of the compement involved, if all are involved the value is None

Interaction 1D1D

class mse.INTERACTIONS.interaction.interaction1D1D(aStudy, aDS1d1, aDS1d2, aCompsInvolved=None)

This Class implement a 1D1D interaction

Parameters:

aStudy (study) – the study
aDS1d1 (dynamicalSystem) – the first 1d dynamical system
aDS1d2 (dynamicalSystem) – the second 1d dynamical system
aCompsInvolved=NONE (numpy.array) – the list of the compement involved, if all are involved the value is None

Interaction to Edges connect

class mse.INTERACTIONS.interaction.interactionEdgesConnect(aStudy, aDs1, aDs2, aCompsInvolved=None)

This class implement an interaction between 2 edges sharing one vertex

Parameters:

aStudy (study) – the study
aDs1 (dynamicalSystem) – the first 1d dynamical system
aDs2 (dynamicalSystem) – the second 1d dynamical system
aCompsInvolved=NONE (numpy.array) – the list of the compement involved, if all are involved the value is None

Interaction 2D2D

class mse.INTERACTIONS.interaction.interaction2D2D(aStudy, aDS2d1, aDS2d2, aDS1d=None, aCompsInvolved=None, tol=None)

Presentation of class

more details

Parameters:

aStudy (study) – the study
aDS2d1 (dynamicalSystem) – the first 2d dynamical system
aDS2d2 (dynamicalSystem) – the second 2d dynamical system
aDS1d (dynamicalSystem) – the common 1d dynamical system
aCompsInvolved=NONE (numpy.array) – the list of the compement involved, if all are involved the value is None

Metropolis-Hasting

class mse.bayesianMethods.metro(aStudy, inf, sup, cptmax, initGuess=None)

Simulated Metropolis-Hastings

Parameters:

aStudy (study) – the mse study
inf (numpy.array(float)) – inf boundary of parameter
sup (numpy.array(float)) – sup boundary of parameter
cptmax (int) – max iteration, if cptmax<0, infinite loop
initGuess (array(float)) – the initial guess for the method.

If None, statrt with a random initial guess.

acceptProp(i, w, thetab, LLb)

Function to accept or not the proposal We create this function for 2 reasons 1- user can overload it 2- it is easier to test (w is defined as a parameter to facilitate testing)

Parameters:

int (i) – local index of tab
float (LLb) – float in ]0,1[. Used to accept or decline the proposal
array (thetab) – theta proposal
float – likelihood of theta proposal

likelihood(): Likelihood calculation

priorF(thetaToTest)

Check if the proposals are within the bounds

Parameters:: array (thetaToTest) – the array to test

proposal(oldTheta)

Return a proposal for parameter Theta

Parameters:: array (oldTheta) – the old value of theta

restartMetro()

restart metro

It will first search for the files, then the lastest values and finally call _runMetro

runMetro(): Run the metro method.

saveMetro(tablen)

Method to save ll, theta and ratio in binary file.

The files are copied in case a programme interruption damages a file during writing.

mse.bayesianMethods.readFileMetro(path)

function to read files create by metro.saveMetro()

Parameters:: path (str) – path to LL, Theta and Ratio files

:return

LL, (numpy.array) - array containing LL
Theta, (numpy.array) - array containing Theta
Ratio, (numpy.array) - array containing Ratio

Linear algebra

modulde mse.linAlgTools contains an programming interface on algebra tools based on petsc.

class mse.linAlgTools.matrix(aMat=None)

The matrix classe define the programming interface on a matrix implementation, here petsc.

Parameters:: aMat (obj) – A petsc mat

assembleMat(aMat, beginN, beginM)

Block addition of aMat starting at row beginN, column beginM

:param matrix aMat :param int beginN :param int beginM

setDim(aN, aM)

It builds an empty sparce matrix.

Warning

be carful

Parameters:

aN (int) – number of row
aM (int) – number of column