comet.snmr.modelling package

Submodules

comet.snmr.modelling.errors module

Part of comet/snmr/modelling

comet.snmr.modelling.errors.DepricationWarning(msg=None)

exception comet.snmr.modelling.errors.InputError(file=None, msg=None)

Bases: Exception

exception comet.snmr.modelling.errors.KernImportError(value)

Bases: Exception

comet.snmr.modelling.mrs module

Part of comet/snmr/modelling

Magnetic resonance sounding module.

class comet.snmr.modelling.mrs.MRS(survey=None, fid=0, kernel=None, mtype='smooth', dtype='rotatedAmplitudes', **kwargs)

Bases: comet.snmr.modelling.nmr_base.SNMRBase

Magnetic resonance sounding (MRS) manager class.

searchForLambda(startLam=20000)

Runs several inversion runs to find the highest lambda which is able to fit the data within its errors.

showDataAndError(ax=None, figsize=(10, 8), as_log=False)

Show data cube along with error cube.

showDataAndFit(compare_to=None, figsize=(8, 6), savename=None, clim=None, suptitle=None, separated=False, savematrices=False)

data and error weighted misfit. 1,1 or 2,2 for complex

showKernel(ax=None, save=None, **kwargs)

Show the kernel as matrix (Q over z). If Kernel is a class object, the plotting order is redirected to Kernel.show(**kwargs)

To see more about the plotting options type this in the console:

>>> import kernel as k
>>> help(k.Kernel.show)

showResult(figsize=(10, 8), save='', fig=None, ax=None, syn=None, wclabel=None, t2label=None, color=None)

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

showResultAndFit(figsize=(12, 10), save='', maxdep=0.0, clim=None, suptitle=None, syn=None, wclabel=None, t2label=None)

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

splitModel(model=None)

Split model vector into d, theta and T2*.

class comet.snmr.modelling.mrs.MRSGenetic(*args, **kwargs)

Bases: comet.snmr.modelling.mrs.MRS

MRS class derivation using a genetic algorithm for inversion.

genMod(individual)

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

plotEAstatistics(fname=None)

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

plotPopulation(maxfitness=None, fitratio=1.05, savefile=True)

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

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

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

comet.snmr.modelling.mrs.showErrorBars(ax, thk, val, thkL, thkU, valL, valU, *args, **kwargs)

Plot wc and t2 models with error bars.

comet.snmr.modelling.mrs.showT2(ax, thk, t2, maxdep=0.0, label=None, color='g')

Show T2 function nicely.

comet.snmr.modelling.mrs.showWC(ax, thk, wc, maxdep=0.0, dw=0.1, label=None, color='g')

Show water content function nicely.

comet.snmr.modelling.mrs_survey module

Part of comet/snmr/modelling

class comet.snmr.modelling.mrs_survey.MRT(survey=None, dim=2, dtype='complex', mtype='smooth')

Bases: object

create1DKernelMesh(verbose=False)

createFOP(kernelmesh=None, secondary=False, para_mesh_2d=None, **kwargs)

kwargs: order (h refinement order for kernel mesh)

createFOPMesh()

createINV(lam=1000, verbose=True, debug=False, **kwargs)

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

lam: float [100]

Lambda factor for inversion.

verbose: bool [True]

Additional verbose decission, can be True, even if the rest of the Manager should remain silent. Most information of the different iterations is printed in the console. It’s recommended to set verbose in this case to True (default).

lambdaFactor, float [0.8]

Sets lambda factor for Marquardt scheme.

robust, bool [False]

Sets the robust flag for the data. See pg.RInversion for more details

logTrans, bool [True]

Applies a logarithmic transformation to the data. Its recommended to do so (default), due to the dealing with water contents, which can’t be negative. Logarithmic transformation is the easiest way to archieve that.

blockyModel, bool [False]

Instead of the standard L2-Norm a L1 Norm can be used to allow for more blocky models. Heavy changes in watercontent and relaxation times can sometimes be fitted better this way.

data

Concatenated data vectors of sounds.

dataIndices

data_slices

Slices to get single data from self.data.

Data[sound #2] = mrt.data[mrt.data_slices[1]]

dtype

error

Concatenated error vectors of sounds.

getSingleDataAndError(sounding_idx)

initSoundings(override=False)

Extends the sounding list for the fids in survey. Called automatically is necessary.

kernels

List with underlaying kernels from sounds.

loadResults(basename, gates=True, pulses=True)

returns (model, error, response, chi2)

mtype

saveResults(basename)

Saves orig data, model, error and forward model as well as chi2.

setDataAndErrorCube(data, error, phase, df=None)

Depricated!

Set data and error cubes using the methods of the single soundings.

Input has to be a list or iterable object of data, and error cubes (pulses x gates) a corresponding list of phase vectors for each pulse and a float defining the frequency offset per sounding.

setDataAndErrorVector(data, error=None, phi=None, df=None)

Depricated! Set Data and Error in MRT and the underlaying MRS instances.

setDataType(dtype)

setKernelMesh(mesh, order=1, **kwargs)

setKernels(basename, load_loopmesh=False, use_order_refinement=True, indices=None)

Sets the kernels for the underlaying soundings. Basename will be formatted with index. Example 5 soundings, basename = ‘kern_{}’ will result in import of kernel_0, kernel_1, …, kernel_4.

setModelTrans(thk=(10, 1, 30, 'log'), wc=(0.3, 0.0, 0.7, 'cot'), t2=(0.2, 0.005, 1.0, 'log'))

Sets model transformation for water content, relaxation times, and thickness (1D). input = (startvalue, min, max, type). Known types are cotangens (‘cot’) and logarithmic (‘log’) transformations.

setModelType(mtype)

setSurvey(survey)

Defines the survey that holds the various soundings and datasets.

setZWeight(z_weight)

showFids(to_plot='abs', rows_cols=None, ax=None, draw='data', **kwargs)

kwargs to ph.plot.drawFID(**kwargs)

showSounding(index, ax=None, to_plot='abs', draw='data', figsize=(5, 3), **kwargs)

Shows Data, Error or misfit of a site.

simulate(model, error, samplingrate=1000.0, max_time=1.0, num_gates=50, verbose=False, **kwargs)

updateData()

Update data vector in inversion instance.

comet.snmr.modelling.nmr_base module

Part of comet/snmr/modelling

Nuclear magnetic resonance base manager as used by MRS and MRT manager classes

class comet.snmr.modelling.nmr_base.SNMRBase(survey=None, fid=0, kernel=None, mtype='block', dtype='rotatedAmplitudes', update_kernel=False, dim=1, **kwargs)

Bases: object

Manager base class for MRS and MRT manager classes.

K

applyBoundsAndTrans()

Append the previously given bounaries for the model transformation to the forward operator.

calcModelCovarianceMatrix()

Compute linear model covariance matrix.

calcModelCovarianceMatrixBounds()

Compute model bounds using covariance matrix diagonals.

createFOP(nlay=3)

Creates the forward operator (FOP). Two possibilities are supported: block and smooth. The choice affects the inversion process and therefore its results.

possible keyword argument for block FOP is ‘nlay’ to define the number of layers, the FOP is calculating (default = 3).

createINV(lam=1000, **kwargs)

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

Parameters:

• lam (float [100]) – Lambda factor for inversion.
• verbose (bool [True]) – Additional verbose decission, can be True, even if the rest of the Manager should remain silent. Most information of the different iterations is printed in the console. It’s recommended to set verbose in this case to True (default).
• special kwargs for Marquardt scheme (block)
• ——————————————-
• lambdaFactor, float [0.8] – Sets lambda factor for Marquardt scheme.
• robust, bool [False] – Sets the robust flag for the data. See pg.RInversion for more details
• special kwargs for smooth scheme
• ———————————
• logTrans, bool [True] – Applies a logarithmic transformation to the data. Its recommended to do so (default), due to the dealing with water contents, which can’t be negative. Logarithmic transformation is the easiest way to archieve that.
• blockyModel, bool [False] – Instead of the standard L2-Norm a L1 Norm can be used to allow for more blocky models. Heavy changes in watercontent and relaxation times can sometimes be fitted better this way.

data

Data Vector representation with respect to self.dtype. Returns None if no sounding or data_cube in sounding is found.

dtype

error

fid

Reference to sounding (FID) class instance in survey.

invert(data=None, error=None, phase=None, lam=1000, runChi1=False, **kwargs)

# TODO!

loadMRSI(filename, verbose=True)

Load data, error and kernel from mrsi file

loadSurvey(dataname)

mtype

setBoundsAndTrans(thkBounds=[10.0, 1.0, 30.0], wcBounds=[0.3, 0.0, 0.7], t2Bounds=[0.2, 0.005, 1.0], trans=['log', 'cot', 'log'])

Sets the boundarys and transformation for the model domain.

Parameters:

• thkBounds (list of floats [ [10., 1., 30.] ]) – Startvalue, lower and upper boundary for thickness of each layer in 1D. Ignored for smooth models (or 2D).
• wcBounds (list of floats [ [0.3, 0.0, 0.7] ]) – Startvalue, lower and upper boundary for water content.
• t2Bounds (list of floats [ [0.2, 0.005, 1.0] ]) – Startvalue, lower and upper boundary for relaxation times.
• trans (list of strings [ [‘log’, ‘cot’, ‘log’] ]) – Defines the type of model transformation. logarithmic (‘log’) or cotangens (‘cot’)

setDataType(dtype)

setKernel(kernelfile=None, load_loopmesh=True, load_kernelmesh=True, use_order_refinement=True)

Load or initialize a new Kernel class instance for calcualting the NMR kernels.

setModelType(mtype)

setSurvey(survey, fid=0)

showCube(ax=None, vec=None, islog=None, clim=None, clab=None, cmap='viridis', cbar=True)

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

simulate(model, err=2.5e-07, samplingrate=1000.0, max_time=1.0, num_gates=50, verbose=False, debug=False, **kwargs)

Creates forward operator and calculates a synthetic response to a given model. Keyword arguments are passed to the function createFOP and to FOP.response. You can also define the ‘Type’ to be ‘smooth’ or ‘block’ or let the simulate function analyse the input.

returns datacube, errorcube (both complex) and phaseinformation (for rotated amplitudes)

Parameters:model (list of lists) – Given model of shape [water_content, relaxation_time] if forwarded to FOP to generate synthetic data set.

comet.snmr.modelling.nmr_base.effectiveNoise(area, noise_lvl=0.0036, sample_rate=1000.0, time=1.0)

Calculates the effective noise of a loop for simulation.

noise_lvl = 3.6e-3 nV / m² / sqrt(number_of_samples) This is a standart noise_lvl from measurements in Schillerslage, Germany. Output in Volt.

comet.snmr.modelling.nmr_base.getPhiByGridSearch(data)

comet.snmr.modelling.smooth_syn module

Part of comet/snmr/modelling

comet.snmr.modelling.smooth_syn.archie(porosity, saturation, water_resistivity, tortuosity=1.0, cementation=1.3, saturation_exponent=2.0, formation_factor=None)

porosity(z), saturation(z)

returns resititvity_bulk(z) cite{}

comet.snmr.modelling.smooth_syn.brooksCorey(z, water_table, porosity, lam=1.6, height_zero=0.12)

after Brooks and Corey (1964) cite{}

comet.snmr.modelling.smooth_syn.costabel(saturation, t2_saturation, lam=1.6)

cite{costabel2011NSG} Costabel, S., and U. Yaramanci, 2011, Relative hydraulic conductivity and effective saturation from Earth’s field nuclear magnetic resonance – a method for assessing the vadose zone: Near Surface Geophysics, 9, 155–167.

comet.snmr.modelling.smooth_syn.effectiveSaturationToWater(saturation_eff, water_saturation, water_residual=0.05)

saturation_eff = (water - water_residual)/ (water_saturated - water_residual)

comet.snmr.modelling.smooth_syn.modelVadose(z, water_table, porosity, t2_saturated, water_resistivity, height_zero=0.12, water_residual=0.05, lam=1.6, verbose=False, **kwargs)

Calculates a synthetical vadose zone on basis of a Brooks-Corey model for saturation over the vadose zone, whereas lambda is the pore size distribution index.

Also calculates the electrical resistivity(z) via Archies law, as well as the distribution of relaxation times based on Costabel and Yaramanci (2011).

returns (z, resistivity, water_content, relaxation_times)

comet.snmr.modelling.smooth_syn.test_local()

comet.snmr.modelling.snmrModelling module

Part of comet/snmr/modelling

Modelling classes for core magnetic resonance (1D, 2D)

class comet.snmr.modelling.snmrModelling.MRS1dBlockQTModelling(survey, fid=0, nlay=3, dtype='complex', kernel=None)

Bases: sphinx.ext.autodoc.importer._MockObject

MRS1dBlockQTModelling - pygimli modelling class for block-mono QT inversion

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

fid

forward(par, verbose=False, num_cpu=12)

yield model response cube as vector

iscomplex

response(par)

class comet.snmr.modelling.snmrModelling.SNMRJointModelling(mrt=None, verbose=False)

Bases: sphinx.ext.autodoc.importer._MockObject

Joint modelling operator for multiple transmitter receiver combinations

addFOP(*fops)

createJacobian(model)

forward(model)

response(model)

setFOPs(fops)

class comet.snmr.modelling.snmrModelling.SNMRModelling(survey, kernel, fid=0, dtype='complex', mesh=None, num_cpu=12, update_kernel=False)

Bases: sphinx.ext.autodoc.importer._MockObject

Modelling class for surface nuclear magnetic resonance (SNMR).

The class is based on the ModellingBase class of pygimli and therefore contains a various amount of parameters and functions as well as some protected members to ensure a generalized interface suitable for the pygimli inversion engine.

For further details about the spezifications of the modelling base, be referred to the pygimli API available from the official project website www.pygimli.org.

static amplitudeJacobian(Mcomplex, model)

calculateKernel(matrix=False, interpolate=False, forceNew=False, **kwargs)

createJacobian(model=None, **kwargs)

Caculate the Jacobian Matrix of a NMR Kernel, with or without relaxation times included (model dependancy for this).

kwargs are redirected to kernelClass.calculate()

Example

>>> # complex jacobian without relaxation time
>>> FOP = MRModelling('a valid kernel class')
>>> FOP.createJacobian()  # sets FOP.jacobian

dimension

fid

Reference to sounding (FID) class instance in survey.

forward(model, **kwargs)

Forward response of the kernel to a specific destribution of watercontent or relaxation times.

model.shape:

array.shape = 2 or 3, numLayers (watercontent only: 2, 3 with relaxation times), first entry = thickness

>>> thickness = [1, 5, 10]
>>> # first layer 0...1 m
>>> # second layer 1...6 m
>>> # third layer 6...16 m
>>> # after that homogeneous halfspace
>>> watercontent = [0.2, 0.3, 0.1, 0.2]  # 1 == 100%
>>> # one entry more than thickness, last entry for halfspace
>>> model = np.array((thickness,
>>>                  watercontent,
>>>                  [100, 200, 14, 100]))  # relaxation times
>>> measurement = mrs.response(model)

iscomplex

jshape

(data, model) == (pulses * gates, model * number of parameters)

Type:Jacobian shape

kshape

(data, model) == (pulses, model)

Type:Kernel shape

response(model)

Calculates the forward response of a SNMR measurement, and returns an 1D numpy array containing the real and imaginary parts of the response. One Voltage value for each pulse moment q and time gate g.

data type: complex

([real(V_11), ..., real(V_1Q),
  real(V_21), ..., real(V_2Q),
  ...,
  real(V_N1), ..., real(V_NQ),
  imag(V_11), ..., ...,
  ...       , ..., imag(V_NQ)]), shape: (2*N, Q)

data type: not complex

([abs(V_11), ..., abs(V_1Q),
  abs(V_21), ..., abs(V_2Q),
  ...,
  abs(V_N1), ..., abs(V_NQ)]), shape: (N, Q)

setKernel(kernel)

setModel(model)

Parameters:model (array_like) – Array that contains three array_like objects. First the thickness of the different layers (number of layers - 1). The second and third array contains the water contents and relaxation times of each layer.

Example

>>> FOP = SNMRModelling('a valid kernel class')
>>> model = [[5., 10.],  # thickness [m]
>>>          [0.1, 0.25, 0.4],  # water content [1]
>>>          [0.1, 0.1, 1.]]  # relaxation times [s]
>>> FOP.setModelVec(model)

setSurvey(survey, fid=0)

updateDataPhase()

Sets data phase for complex inversion. If no model is given the starting model is used.

vector

Module contents

Part of comet/snmr/modelling