pygimli.physics.sNMR

Surface nuclear magnetic resonance (NMR) data inversion

Overview

Classes

MRS([name, verbose]) Magnetic resonance sounding (MRS) manager class.
MRS1dBlockQTModelling(nlay, K, zvec, t[, …]) MRS1dBlockQTModelling - pygimli modelling class for block-mono QT inversion
MRSprofile([filename, x, dx, x0]) manager class for several MRS data along a profile for joint inversion

Classes

MRS

class pygimli.physics.sNMR.MRS(name=None, verbose=True, **kwargs)[source]

Magnetic resonance sounding (MRS) manager class.

Attributes:
t, q : ndarray - time and pulse moment vectors
data, error : 2d ndarray - data and error cubes
K, z : ndarray - (complex) kernel and its vertical discretization
model, modelL, modelU : vectors - model vector and lower/upper bound to it

Methods

loadMRSI - load MRSI (MRSmatlab format) data  
showCube - show any data/error/misfit as data cube (over q and t)  
showDataAndError - show data and error cubes  
showKernel - show Kernel matrix  
createFOP - create forward operator  
createInv - create pygimli Inversion instance  
run - run block-mono (alternatively smooth-mono) inversion (with bootstrap)  
calcMCM - compute model covariance matrix and thus uncertainties  
splitModel - return thickness, water content and T2* time from vector  
showResult/showResultAndFit - show inversion result (with fit)  
runEA - run evolutionary algorithm (GA, PSO etc.) using inspyred  
plotPopulation - plot final population of an EA run  
__init__(name=None, verbose=True, **kwargs)[source]

MRS init with optional data load from mrsi file

Parameters:
name : string

Filename with load data and kernel (.mrsi) or just data (.mrsd)

verbose : bool

be verbose

kwargs - see :func:`MRS.loadMRSI`.
calcMCM()[source]

Compute linear model covariance matrix.

calcMCMbounds()[source]

Compute model bounds using covariance matrix diagonals.

checkData(**kwargs)[source]

Check data and retrieve data and error vector.

static createFOP(nlay, K, z, t)[source]

Create forward operator instance.

createInv(nlay=3, lam=100.0, verbose=True, **kwargs)[source]

Create inversion instance (and fop if necessary with nlay).

genMod(individual)[source]

Generate (GA) model from random vector (0-1) using model bounds.

invert(nlay=3, lam=100.0, startvec=None, verbose=True, uncertainty=False, **kwargs)[source]

Easiest variant doing all (create fop and inv) in one call.

loadDataCube(filename='datacube.dat')[source]

Load data cube from single ascii file (old stuff)

loadDataNPZ(filename, **kwargs)[source]

Load data and kernel from numpy gzip packed file.

The npz file contains the fields: q, t, D, (E), z, K

loadDir(dirname)[source]

Load several standard files from dir (old Borkum stage).

loadErrorCube(filename='errorcube.dat')[source]

Load error cube from a single ascii file (old stuff).

loadKernel(name='')[source]

Load kernel matrix from mrsk or two bmat files.

loadKernelNPZ(filename, **kwargs)[source]

Load data and kernel from numpy gzip packed file.

The npz file contains the fields: q, t, D, (E), z, K

loadMRSD(filename, usereal=False, mint=0.0, maxt=2.0)[source]

Load mrsd (MRS data) file: not really used as in MRSD.

loadMRSI(filename, **kwargs)[source]

Load data, error and kernel from mrsi or mrsd file

Parameters:
usereal : bool [False]

use real parts (after data rotation) instead of amplitudes

mint/maxt : float [0.0/2.0]

minimum/maximum time to restrict time series

loadResult(filename)[source]

Load inversion result from column file.

loadZVector(filename='zkernel.vec')[source]

Load the kernel vertical discretisation (z) vector.

plotEAstatistics(fname=None)[source]

Plot EA statistics (best, worst, …) over time.

plotPopulation(maxfitness=None, fitratio=1.05, savefile=True)[source]

Plot fittest individuals (fitness<maxfitness) as 1d models

Parameters:
maxfitness : float

maximum fitness value (absolute) OR

fitratio : float [1.05]

maximum ratio to minimum fitness

result()[source]

Return block model results (thk, wc and T2 vectors).

run(verbose=True, uncertainty=False, **kwargs)[source]

Easiest variant doing all (create fop and inv) in one call.

runEA(nlay=None, eatype='GA', pop_size=100, num_gen=100, runs=1, mp_num_cpus=8, **kwargs)[source]

Run evolutionary algorithm using the inspyred library

Parameters:
nlay : int [taken from classic fop if not given]

number of layers

pop_size : int [100]

population size

num_gen : int [100]

number of generations

runs : int [pop_size*num_gen]

number of independent runs (with random population)

eatype : string [‘GA’]
algorithm, choose among:

‘GA’ - Genetic Algorithm [default] ‘SA’ - Simulated Annealing ‘DEA’ - Discrete Evolutionary Algorithm ‘PSO’ - Particle Swarm Optimization ‘ACS’ - Ant Colony Strategy ‘ES’ - Evolutionary Strategy

saveFigs(basename=None, extension='pdf')[source]

Save all figures to (pdf) files.

saveResult(filename)[source]

Save inversion result to column text file for later use.

setBoundaries()[source]

Set parameter boundaries for inversion.

showCube(ax=None, vec=None, islog=None, clim=None, clab=None)[source]

Plot any data (or response, error, misfit) cube nicely.

showDataAndError(figsize=(10, 8), show=False)[source]

Show data cube along with error cube.

showKernel(ax=None)[source]

Show the kernel as matrix (Q over z).

showResult(figsize=(10, 8), save='', fig=None, ax=None)[source]

Show theta(z) and T2*(z) (+uncertainties if there).

showResultAndFit(figsize=(12, 10), save='', plotmisfit=False, maxdep=0, clim=None)[source]

Show ec(z), T2*(z), data and model response.

static simulate(model, K, z, t)[source]

Do synthetic modelling.

splitModel(model=None)[source]

Split model vector into d, theta and T2*.

MRS1dBlockQTModelling

class pygimli.physics.sNMR.MRS1dBlockQTModelling(nlay, K, zvec, t, verbose=False)[source]

MRS1dBlockQTModelling - pygimli modelling class for block-mono QT inversion

f=MRS1dBlockQTModelling(lay, KR, KI, zvec, t, verbose = False )

Methods

__call__((object)arg1, (object)model) C++ signature :
clearConstraints((object)arg1) C++ signature :
clearJacobian((object)arg1) C++ signature :
constraints((object)arg1) C++ signature :
constraintsRef((object)arg1) C++ signature :
createConstraints((object)arg1) C++ signature :
createDefaultStartModel((object)arg1) C++ signature :
createJacobian((object)arg1, (object)model) C++ signature :
createMappedModel((object)arg1, …) Read only extrapolation of model values given per cell marker to values given per cell.
createRefinedForwardMesh((object)arg1 [, …) C++ signature :
createStartModel((object)arg1) C++ signature :
createStartVector((object)arg1) DEPRECATED use createStartModel
data((object)arg1) Return the associated data container.
deleteMesh((object)arg1) Delete the actual mesh.
initConstraints((object)arg1) C++ signature :
initJacobian((object)arg1) C++ signature :
initRegionManager((object)arg1) C++ signature :
jacobian((object)arg1) Return the pointer to the Jacobian matrix associated with this forward operator.
jacobianRef((object)arg1) C++ signature :
mapModel((object)arg1, (object)model [, …) C++ signature :
mesh((object)arg1) C++ signature :
multiThreadJacobian((object)arg1) Return number of threads used for Jacobian generation.
region((object)arg1, (object)marker) Syntactic sugar for this->regionManager().region(marker).
regionManager((object)arg1) C++ signature :
regionManagerRef((object)arg1) C++ signature :
response(par) Yield model response cube as vector.
response_mt((object)arg1, (object)model [, …) C++ signature :
setConstraints((object)arg1, (object)C) C++ signature :
setData((object)arg1, (object)data) Change the associated data container
setJacobian((object)arg1, (object)J) C++ signature :
setMesh((object)arg1, (object)mesh [, …) Set new mesh to the forward operator, optionally hold region parameter for the new mesh (i.e.
setMultiThreadJacobian((object)arg1, …) Set number of threads used for brute force Jacobian generation.
setRegionManager((object)arg1, (object)reg) C++ signature :
setStartModel((object)arg1, (object)startModel) C++ signature :
setThreadCount((object)arg1, (object)nThreads) Set the maximum number of allowed threads for MT calculation.
setVerbose((object)arg1, (object)verbose) Set verbose state.
solution((object)arg1) C++ signature :
startModel((object)arg1) C++ signature :
threadCount((object)arg1) Return the maximum number of allowed threads for MT calculation
verbose((object)arg1) Get verbose state.
createJacobian_mt  
responses  
__init__(nlay, K, zvec, t, verbose=False)[source]

Constructor with number of layers, kernel, z and t vectors.

clearConstraints((object)arg1) → object :
C++ signature :
void* clearConstraints(GIMLI::ModellingBase {lvalue})

clearConstraints( (object)arg1) -> object :

C++ signature :
void* clearConstraints(ModellingBase_wrapper {lvalue})
clearJacobian((object)arg1) → object :
C++ signature :
void* clearJacobian(GIMLI::ModellingBase {lvalue})

clearJacobian( (object)arg1) -> object :

C++ signature :
void* clearJacobian(ModellingBase_wrapper {lvalue})
constraints((object)arg1) → object :
C++ signature :
GIMLI::MatrixBase* constraints(GIMLI::ModellingBase {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(ModellingBase_wrapper {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(GIMLI::ModellingBase {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(ModellingBase_wrapper {lvalue})
constraintsRef((object)arg1) → object :
C++ signature :
GIMLI::SparseMapMatrix<double, unsigned long> {lvalue} constraintsRef(GIMLI::ModellingBase {lvalue})

constraintsRef( (object)arg1) -> object :

C++ signature :
GIMLI::SparseMapMatrix<double, unsigned long> {lvalue} constraintsRef(GIMLI::ModellingBase {lvalue})
createConstraints((object)arg1) → object :
C++ signature :
void* createConstraints(GIMLI::ModellingBase {lvalue})

createConstraints( (object)arg1) -> object :

C++ signature :
void* createConstraints(ModellingBase_wrapper {lvalue})
createDefaultStartModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> createDefaultStartModel(GIMLI::ModellingBase {lvalue})

createDefaultStartModel( (object)arg1) -> object :

C++ signature :
GIMLI::Vector<double> createDefaultStartModel(ModellingBase_wrapper {lvalue})
createJacobian((object)arg1, (object)model) → object :
C++ signature :
void* createJacobian(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>)

createJacobian( (object)arg1, (object)model) -> object :

C++ signature :
void* createJacobian(ModellingBase_wrapper {lvalue},GIMLI::Vector<double>)

createJacobian( (object)arg1, (object)model, (object)resp) -> object :

C++ signature :
void* createJacobian(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>,GIMLI::Vector<double>)

createJacobian( (object)arg1, (object)model, (object)resp) -> object :

C++ signature :
void* createJacobian(ModellingBase_wrapper {lvalue},GIMLI::Vector<double>,GIMLI::Vector<double>)
createJacobian_mt(model, resp)
createMappedModel((object)arg1, (object)model[, (object)background=-1]) → object :

Read only extrapolation of model values given per cell marker to values given per cell. Exterior values will be set to background or prolongated for background = -1.

C++ signature :
GIMLI::Vector<double> createMappedModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,double=-1])
createRefinedForwardMesh((object)arg1[, (object)refine=True[, (object)pRefine=False]]) → object :
C++ signature :
void* createRefinedForwardMesh(GIMLI::ModellingBase {lvalue} [,bool=True [,bool=False]])
createStartModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> createStartModel(GIMLI::ModellingBase {lvalue})
createStartVector((object)arg1) → object :

DEPRECATED use createStartModel

C++ signature :
GIMLI::Vector<double> createStartVector(GIMLI::ModellingBase {lvalue})
data((object)arg1) → object :

Return the associated data container.

C++ signature :
GIMLI::DataContainer {lvalue} data(GIMLI::ModellingBase {lvalue})
deleteMesh((object)arg1) → object :

Delete the actual mesh.

C++ signature :
void* deleteMesh(GIMLI::ModellingBase {lvalue})
initConstraints((object)arg1) → object :
C++ signature :
void* initConstraints(GIMLI::ModellingBase {lvalue})

initConstraints( (object)arg1) -> object :

C++ signature :
void* initConstraints(ModellingBase_wrapper {lvalue})
initJacobian((object)arg1) → object :
C++ signature :
void* initJacobian(GIMLI::ModellingBase {lvalue})

initJacobian( (object)arg1) -> object :

C++ signature :
void* initJacobian(ModellingBase_wrapper {lvalue})
initRegionManager((object)arg1) → object :
C++ signature :
void* initRegionManager(GIMLI::ModellingBase {lvalue})
jacobian((object)arg1) → object :

Return the pointer to the Jacobian matrix associated with this forward operator.

C++ signature :
GIMLI::MatrixBase* jacobian(GIMLI::ModellingBase {lvalue})
jacobian( (object)arg1) -> object :

Return the pointer to the Jacobian matrix associated with this forward operator.

C++ signature :
GIMLI::MatrixBase* jacobian(GIMLI::ModellingBase {lvalue})
jacobianRef((object)arg1) → object :
C++ signature :
GIMLI::Matrix<double> {lvalue} jacobianRef(GIMLI::ModellingBase {lvalue})

jacobianRef( (object)arg1) -> object :

C++ signature :
GIMLI::Matrix<double> {lvalue} jacobianRef(GIMLI::ModellingBase {lvalue})
mapModel((object)arg1, (object)model[, (object)background=0]) → object :
C++ signature :
void* mapModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,double=0])
mesh((object)arg1) → object :
C++ signature :
GIMLI::Mesh* mesh(GIMLI::ModellingBase {lvalue})
multiThreadJacobian((object)arg1) → object :

Return number of threads used for Jacobian generation.

C++ signature :
unsigned long multiThreadJacobian(GIMLI::ModellingBase {lvalue})
region((object)arg1, (object)marker) → object :

Syntactic sugar for this->regionManager().region(marker).

C++ signature :
GIMLI::Region* region(GIMLI::ModellingBase {lvalue},int)
regionManager((object)arg1) → object :
C++ signature :
GIMLI::RegionManager regionManager(GIMLI::ModellingBase {lvalue})

regionManager( (object)arg1) -> object :

C++ signature :
GIMLI::RegionManager {lvalue} regionManager(GIMLI::ModellingBase {lvalue})
regionManagerRef((object)arg1) → object :
C++ signature :
GIMLI::RegionManager {lvalue} regionManagerRef(GIMLI::ModellingBase {lvalue})
response(par)[source]

Yield model response cube as vector.

response_mt((object)arg1, (object)model[, (object)i=0]) → object :
C++ signature :
GIMLI::Vector<double> response_mt(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,unsigned long=0])

response_mt( (object)arg1, (object)model [, (object)i=0]) -> object :

C++ signature :
GIMLI::Vector<double> response_mt(ModellingBase_wrapper {lvalue},GIMLI::Vector<double> [,unsigned long=0])
responses(models, respos)
setConstraints((object)arg1, (object)C) → object :
C++ signature :
void* setConstraints(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)

setConstraints( (object)arg1, (object)C) -> object :

C++ signature :
void* setConstraints(ModellingBase_wrapper {lvalue},GIMLI::MatrixBase*)
setData((object)arg1, (object)data) → object :

Change the associated data container

C++ signature :
void* setData(GIMLI::ModellingBase {lvalue},GIMLI::DataContainer {lvalue})
setJacobian((object)arg1, (object)J) → object :
C++ signature :
void* setJacobian(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)

setJacobian( (object)arg1, (object)J) -> object :

C++ signature :
void* setJacobian(ModellingBase_wrapper {lvalue},GIMLI::MatrixBase*)
setMesh((object)arg1, (object)mesh [, (object)ignoreRegionManager=False]) -> object : Set new mesh to the forward operator, optionally hold region parameter for the new mesh (i.e. for roll a long)

Set new mesh to the forward operator, optionally hold region parameter for the new mesh (i.e. for roll a long)

C++ signature :
void* setMesh(GIMLI::ModellingBase {lvalue},GIMLI::Mesh [,bool=False])
setMultiThreadJacobian((object)arg1, (object)nThreads) → object :

Set number of threads used for brute force Jacobian generation. 1 is default. If nThreads is greater than 1 you need to implement response_mt with a read only response function. Maybe its worth set the single setThreadCount to 1 than, that you dont find yourself in a threading overkill.

C++ signature :
void* setMultiThreadJacobian(GIMLI::ModellingBase {lvalue},unsigned long)
setRegionManager((object)arg1, (object)reg) → object :
C++ signature :
void* setRegionManager(GIMLI::ModellingBase {lvalue},GIMLI::RegionManager*)
setStartModel((object)arg1, (object)startModel) → object :
C++ signature :
void* setStartModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>)

setStartModel( (object)arg1, (object)startModel) -> object :

C++ signature :
void* setStartModel(ModellingBase_wrapper {lvalue},GIMLI::Vector<double>)
setThreadCount((object)arg1, (object)nThreads) → object :

Set the maximum number of allowed threads for MT calculation. Have to be greater than 0. Will also set ENV(OPENBLAS_NUM_THREADS) .. if used.

C++ signature :
void* setThreadCount(GIMLI::ModellingBase {lvalue},unsigned long)
setVerbose((object)arg1, (object)verbose) → object :

Set verbose state.

C++ signature :
void* setVerbose(GIMLI::ModellingBase {lvalue},bool)
solution((object)arg1) → object :
C++ signature :
GIMLI::Matrix<double> solution(GIMLI::ModellingBase {lvalue})
startModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> startModel(GIMLI::ModellingBase {lvalue})

startModel( (object)arg1) -> object :

C++ signature :
GIMLI::Vector<double> startModel(ModellingBase_wrapper {lvalue})
threadCount((object)arg1) → object :

Return the maximum number of allowed threads for MT calculation

C++ signature :
unsigned long threadCount(GIMLI::ModellingBase {lvalue})
verbose((object)arg1) → object :

Get verbose state.

C++ signature :
bool verbose(GIMLI::ModellingBase {lvalue})

verbose( (object)arg1) -> object :

C++ signature :
bool verbose(GIMLI::ModellingBase {lvalue})

MRSprofile

class pygimli.physics.sNMR.MRSprofile(filename=None, x=None, dx=1, x0=0, **kwargs)[source]

manager class for several MRS data along a profile for joint inversion

Attributes:
mrs : list of MRS objects (single soundings)
x : list of positions for the soundings

Methods

load - load mrs files from a directory  
set X - set x vector  
showData - show MRS data  
independentBlock1dInversion - perform independent 1D block inversion  
block1dInversion - 1D block inversion of all data sets together  
blockLCInversion - 1D block laterally constrained inversion of all data  
printFits - print total misfit (chi^2, rms) and individual values  
showModel - show LCI model  
__init__(filename=None, x=None, dx=1, x0=0, **kwargs)[source]

Initialize profile object by mrs objects and optional positions.

Parameters:
filename : list of str | str

list of files OR filenames(with *) OR directory to load

x : iterable

position vector of individual soundings

x0 : float [0]

starting position

dx : position [1]

position increment

block1dInversion(nlay=2, lam=100.0, show=False, verbose=True, uncertainty=False)[source]

Invert all data together by a 1D model (more general solution).

block1dInversionOld(nlay=2, startModel=None, verbose=True, uncertainty=False, **kwargs)[source]

Invert all data together by one 1D model (variant 1 - all equal).

blockLCInversion(nlay=2, startModel=None, **kwargs)[source]

Laterally constrained (piece-wise 1D) block inversion.

independentBlock1dInversion(nlay=2, lam=100, startModel=None)[source]

Independent inversion of all soundings.

Parameters:
nlay : int [2]

number of layers

lam : float

regularisation parameter

startModel : array/vector

starting model (see MRS.run parameters)

load(filenames, **kwargs)[source]

load mrs files in a list of (single) MRS handlers filename can be a list of mrsi files or a directory to search Additional parameters: usereal, mint, maxt (see MRS.load)

loadKernel(kernelfile)[source]

load one kernel file for all soundings

printFits()[source]

Show single fits and total fit.

saveFigs(basename='out', extension='pdf')[source]

Save all figures to (pdf) files.

setX(x=None, x0=0, dx=1)[source]

define positions for soundings and sort accordingly

show1dModel()[source]

Show 1D model (e.g. of joint block inversion).

showData(figsize=(15, 10), nc=0, nr=0, clim=None)[source]

show all data cubes in subplots

showFits(ax=None)[source]

Show chi-square and rms fits of individual soundings.

showInitialValues()[source]

show initial values of whole profile

showModel(showFit=0, cmap=<matplotlib.colors.LinearSegmentedColormap object>, figsize=(13, 12), wlim=(0, 0.5), tlim=(0.05, 0.5))[source]

Show 2d model as stitched 1d models along with fit.

showT2(tlim=(0.05, 0.5), ax=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, title='$T_2^*$ (s)')[source]

Show relaxation time distribution as stitched model section.

showWC(wlim=(0, 0.5), ax=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, title='$\\theta$ (-)')[source]

Show water content distribution as stitched model section.



2019 - GIMLi Development Team
Created using Bootstrap, Sphinx and pyGIMLi 1.0.11+23.g5fafebb7 on May 23, 2019.