# Polyfit¶

This tutorial shows a flexible inversion with an own forward calculation that includes an own jacobian. We start with fitting a polynomial of degree $$P$$

$f(x) = p_0 + p_1 x + \ldots + p_P x^P = \sum\limits_{i=0}^{P} p_i x^i$

to given data $$y$$. The unknown model is the coefficient vector $${\bf m}=[p_0,\ldots,p_P]$$. The vectorized function for a vector $${\bf x}=[x_1,\ldots,x_N]^T$$ can be written as matrix-vector product

$\begin{split}{\bf f} ({\bf x}) = {\bf A} {\bf x} \quad\mbox{with}\quad {\bf A}= \left[ \begin{array}{cccc} 1 & x_1 & \ldots & x_1^P \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_N & \ldots & x_N^P \end{array} \right] = [ {\bf 1}\quad {\bf x} \quad {\bf x}^2 \ldots {\bf x}^P ] \;.\end{split}$

We set up the modelling operator, i.e. to return $${\bf f}({\bf x})$$ for given $$p_i$$, as a class derived from the modelling base class. The latter holds the main mimic of generating Jacobian, gradients by brute force. The only function to overwrite is cw{response()}.

Python is a very flexible language for programming and scripting and has many packages for numerical computing and graphical visualization. For this reason, we built Python bindings and compiled the library pygimli. As a main advantage, all classes can be used and derived. This makes the use of GIMLi very easy for non-programmers. All existing modelling classes can be used, but it is also easy to create new modelling classes.

We exemplify this by the preceding example.

First, the library must be imported.

To avoid name clashes with other libraries we suggest to import pygimli and alias it to an easy name (as usually done for numpy or matplotlib), e.g. by

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

The modelling class is derived from ModellingBase, a constructor is defined and the response function is defined. Due to the linearity of the problem we store the matrix $${\bf A}$$, which is also the Jacobian matrix and use it for the forward calculation. A second function is just added as reference. We overwrite the method createJacobian as we know it but do nothing in the actual computation. If $${\bf J}$$ depends on $${\bf m}$$ this function must be filled.

class FunctionModelling(pg.ModellingBase):
def __init__(self, nc, xvec, verbose=False):
pg.ModellingBase.__init__(self, verbose)
self.x_ = xvec
self.nc_ = nc
nx = len(xvec)
self.regionManager().setParameterCount(nc)
self.jacobian().resize(nx, nc)
for i in range(self.nc_):
self.jacobian().setCol(i, pg.pow(self.x_, i))

def response(self, model):
return self.jacobian() * model

def responseDirect(self, model):
y = pg.RVector(len(self.x_), model[0])

for i in range(1, self.nc_):
y += pg.pow(self.x_, i) * model[i]

return y

def createJacobian(self, model):
pass  # if J depends on the model you should work here

def startModel(self):
return pg.RVector(self.nc_, 0.5)

Let us create some synthetic data for some x values

x = np.arange(0., 10., 0.5)
y = 1.1 + 2.1 * x - 0.2 * x**2
noise = 0.1
y += np.random.randn(len(y)) * noise

We now start by setting up the modelling operator, and inversion and run it.

fop = FunctionModelling(3, x)

# initialize inversion with data and forward operator and set options
inv = pg.RInversion(y, fop)

# constant absolute error of 0.01 is 1% (not necessary, only for chi^2)
inv.setAbsoluteError(noise)

# the problem is well-posed and does not need any regularization
inv.setLambda(0)

# actual inversion run yielding coefficient model
coeff = inv.run()

inv.echoStatus()

print(coeff)

Out:

3 [1.1936214920549455, 2.050850357252683, -0.19500213573790504]

The result is easily plotted by

plt.plot(x, y, 'rx', x, inv.response(), 'b-')
plt.show()

Gallery generated by Sphinx-Gallery

Created using Bootstrap and Sphinx. Last updated on Sep 23, 2018.