Simple fitΒΆ

This tutorial shows how to do the simplest inversion case, a curve fit, by setting up an own forward operator. The function to be fitted is`

\[f(x) = A * e^{-x/X}\]

with the two unknown coefficients A (a signal amplitude) and X (a decay rate). Both A and X are assumed to be positive which is often the case for physical properties. The easiest way to do this is via a logarithmic transformation of the model vector (containing A and X) which is very easily done in pyGIMLi.

First we import the pygimli library under a short name pg and the numerics library numpy. Additionally we load the python plotting module of the library matplotlib. Both are contained in most python distributions and systems.

import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt

We set up the modelling operator, i.e. to return \({\bf f}({\bf x})\) for given model parameters A and X subsumed in a vector. In order to be able to use operator in inversion, we derive from the abstract modelling base class. The latter holds the main mimic of generating Jacobian and adminstrating the model, the regularization and so on. The only function to overwrite is response(). If no function createJacobian is provided, they are computed by brute force (forward calculations with altered parameters).

class ExpModelling(pg.ModellingBase):
    def __init__(self, xvec, verbose=False):
        pg.ModellingBase.__init__(self, verbose)
        self.x = xvec
        self.setMesh(pg.createMesh1D(1, 2))
        # self.regionManager().setParameterCount(2)

    def response(self, model):
        return model[0] * pg.exp(-self.x / model[1])

    def startModel(self):
        return [1.0, 0.3]

The init function saves the x vector and defines the parameterization, i.e. two independent parameters (a 1D mesh with 1 cell and 2 properties). The response function computes the function using A=model[0] and X=model[1] The function startModel defines a meaningful starting vector. There are other methods to set the starting model as inv.setModel() but this one is a default one for people who use the class and forget about a starting model.

# We first create an abscissa vector using numpy (note that pygimli also
# provides an exp function and generate synthetic data with two arbitrary A and
# X values.

x = np.arange(0, 1, 1e-2)
data = 10.5 * np.exp(- x / 550e-3)

We define an (absolute) error level and add Gaussian noise to the data.

error = 0.1
data += np.random.randn(*data.shape)*error

Next, an instance of the forward operator is created. We could use it for calculating the synthetic data using f.response([10.5, 0.55]) or just f([10.5, 0.55]). We create a real-valued (R) inversion passing the forward operator, the data. A verbose boolean flag could be added to provide some output the inversion, another one prints more and saves files for debugging.

f = ExpModelling(x)
inv = pg.RInversion(data, f)

We create a real-valued logarithmid transformation and appy it to the model. Similar could be done for the data which are by default treated linearly. We then set the error level that is used for data weighting. It can be a float number or a vector of data length. One can also set a relative error. Finally, we define the inversion style as Marquard scheme (pure local damping with decreasing the regularization parameter subsequently) and start with a relatively large regularization strength to avoid overshoot. Finally run yields the coefficient vector and we plot some statistics.

tLog = pg.RTransLog()
inv.setTransModel(tLog)
inv.setAbsoluteError(error)
inv.setMarquardtScheme()
inv.setLambda(100)
coeff = inv.run()
inv.echoStatus()  # result not shown on live docu
print(inv.absrms(), inv.relrms(), inv.chi2())
print(coeff)

Out:

0.10949324795365084 3.2249010780300535 1.198877134743966
<class 'pygimli.core._pygimli_.RVector'> 2 [10.552048285467066, 0.5463869572556201]

We see that after 5 iterations the absolute rms value equals the noise level corresponding to a chi-squared misfit value of 1 as it should be the case for synthetic data. The relative rms (in %) is less relevant here but can be for other applications. Additionally the ranges for model and model response are given and the objective function consisting of data misfit and model roughness times lambda. Note that due to the local regularization the second term does not contribute to Phi. Set verbose to True to see the whole course of inversion. The values of the coefficient vector (a GIMLi real vector) are as expected close (i.e. equivalent) to the synthetic model.

We finally create a plotting figure and plot both data and model response.

plt.figure()
plt.plot(x, data, 'rx', x, inv.response(), 'b-')
../../_images/sphx_glr_plot-0-expfit_001.png

The createMesh1D automatically attributed the markers 0 and 1 to the two model parameters A and X, respectively. Each marker leads to a region that can be individually treated, e.g. the starting value, lower or upper bounds, or all three at the same time (setParameters). This changes the model transformation which can of course be region-specific.

f.region(0).setLowerBound(0.1)
f.region(0).setStartValue(3)
f.region(1).setParameters(0.3, 0.01, 1.0)

If these are set before the inversion is used, they are used automatically. We set the model by hand using the new starting model

inv.setVerbose(True)
inv.setModel(f.createStartVector())
print(inv.run())
inv.echoStatus()

Out:

<class 'pygimli.core._pygimli_.RVector'> 2 [10.552151487039865, 0.5464002620046111]

The result is pretty much the same as before but for stronger equivalence or smoothness-constrained regularization prior information might help a lot.

Total running time of the script: ( 0 minutes 0.075 seconds)

Generated by Sphinx-Gallery

News (Aug 2017):
pyGIMLi Paper out now! Download PDF here.


© 2017 - GIMLi Development Team
Created using Bootstrap and Sphinx. Last updated on Aug 15, 2017.