geoist.inversion package

Subpackages

• geoist.inversion.tests package
  • Submodules
  • geoist.inversion.tests.test_geometry module
  • geoist.inversion.tests.test_mesh module
  • Module contents

Submodules

geoist.inversion.abic module

class geoist.inversion.abic.AbicLSQOperator(toep, depth_constraint=None, smooth_components={}, refer_constraint=None, weights=None)

基类：geoist.inversion.toeplitz.LSQOperator

An operator doing matrix vector multiplication. The matrix is: $lpha_g G^TG + sum lpha_i W^TB_i^TB_iW$. Where $lpha$'s are weights, $G$ is kernel matrix, $W$ is depth constraint, $B_i$'s are other constrains.

matvec(v)

If do_precond is True, apply T^T*M^{-T}*M^{-1}*T + coef*diag**2 to vector v, else apply T^T*T + coef*diag**2 to vector v

class geoist.inversion.abic.AbicLSQOperator2(op_in)

基类：object

matvec(v)

class geoist.inversion.abic.GravInvAbicModel(nzyx=[4, 4, 4], smooth_components=None, depth_constraint=None, model_density=None, refer_density=None, weights=None, source_volume=None, smooth_on='m', subtract_mean=True, data_dir='/data/gravity_inversion')

基类：object

abic_optimize()

bound_optimize(x0=None)

calc_abic()

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_abic_quiet()

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_log_obs_det()

calc_log_obs_det_quiet()

calc_log_prior_total_det()

calc_log_prior_total_det_quiet()

calc_min_u(solved=False, x=None)

calc_res()

calc_u(solved=False, x=None)

calc_u_quiet(solved=False, x=None)

property depth_constraint

do_linear_solve()

do_linear_solve_quiet()

forward(model_density=None)

gen_kernel(process=1)

gen_mesh(height=- 1)

gen_model_name()

hessp_u(x, v)

jac_u(x=None)

lasso_hessp(x, v)

lasso_jac(x)

lasso_optimize(x0=None)

lasso_target(x)

property model_density

property nx

property ny

property nz

para_grad(x)

print_summary()

property refer_density

property smooth_components

property source_volume

u_bound()

walsh_transform(keys=None)

property weights

class geoist.inversion.abic.SmoothOperator

基类：object

derivation(v, component='dx', array_order='zyx')

rderivation(v, component='dx', array_order='zyx')

geoist.inversion.abic.check_gpu()

geoist.inversion.abic.free_gpu()

free up gpu memory consumption

geoist.inversion.abic.timeit(f)

geoist.inversion.abic_fft module

geoist.inversion.base module

The base classes for inverse problem solving.

See geoist.inversion for examples, regularization, and more.

This module defines base classes that are used by the rest of the inversion package:

• MultiObjective: A "container" class that emulates a sum of different objective (goal) functions (like Misfit or some form of regularization). When two of those classes are added they generate a MultiObjective object.

• OperatorMixin: A mix-in class that defines the operators + and * (by a scalar). Used to give these properties to Misfit and the regularizing functions. Adding results in a MultiObjective. Multiplying sets the regul_param of the class (like a scalar weight factor).

• OptimizerMixin: A mix-in class that defines the fit and config methods for optimizing a Misfit or MultiObjective and fitting the model to the data.

• CachedMethod: A class that wraps a method and caches the returned value. When the same argument (an array) is passed twice in a row, the class returns the cached value instead of recomputing.

• CachedMethodPermanent: Like CachedMethod but always returns the cached value, regardless of the input. Effectively calculates only the first time the method is called. Useful for caching the Jacobian matrix in a linear problem.

---

class geoist.inversion.base.CachedMethod(instance, meth)

基类：object

Wrap a method to cache it's output based on the hash of the input array.

Store the output of calling the method on a numpy array. If the method is called in succession with the same input array, the cached result will be returned. If the method is called on a different array, the old result will be discarded and the new one stored.

Uses SHA1 hashes of the input array to tell if it is the same array.

注解

We need the object instance and method name instead of the bound method (like obj.method) because we can't pickle bound methods. We need to be able to pickle so that the solvers can be passed between processes in parallelization.

Parameters:

• instanceobject

  The instance of the object that has the method you want to cache.

• methstring

  The name of the method you want to cache.

Examples:

>>> import numpy as np
>>> class MyClass(object):
...     def __init__(self, cached=False):
...         if cached:
...             self.my_method = CachedMethod(self, 'my_method')
...     def my_method(self, p):
...         return p**2
>>> obj = MyClass(cached=False)
>>> a = obj.my_method(np.arange(0, 5))
>>> a
array([ 0,  1,  4,  9, 16])
>>> b = obj.my_method(np.arange(0, 5))
>>> a is b
False
>>> cached = MyClass(cached=True)
>>> a = cached.my_method(np.arange(0, 5))
>>> a
array([ 0,  1,  4,  9, 16])
>>> b = cached.my_method(np.arange(0, 5))
>>> a is b
True
>>> cached.my_method.hard_reset()
>>> b = cached.my_method(np.arange(0, 5))
>>> a is b
False
>>> c = cached.my_method(np.arange(0, 5))
>>> b is c
True
>>> cached.my_method(np.arange(0, 6))
array([ 0,  1,  4,  9, 16, 25])

hard_reset()

Delete the cached values.

class geoist.inversion.base.CachedMethodPermanent(instance, meth)

基类：object

Wrap a method to cache it's output and return it whenever the method is called..

This is different from CachedMethod because it will only run the method once. All other times, the result returned will be this first one. This class should be used with methods that should return always the same output (like the Jacobian matrix of a linear problem).

注解

We need the object instance and method name instead of the bound method (like obj.method) because we can't pickle bound methods. We need to be able to pickle so that the solvers can be passed between processes in parallelization.

Parameters:

• instanceobject

  The instance of the object that has the method you want to cache.

• methstring

  The name of the method you want to cache.

Examples:

>>> import numpy as np
>>> class MyClass(object):
...     def __init__(self, cached=False):
...         if cached:
...             self.my_method = CachedMethodPermanent(self, 'my_method')
...     def my_method(self, p):
...         return p**2
>>> obj = MyClass(cached=False)
>>> a = obj.my_method(np.arange(0, 5))
>>> a
array([ 0,  1,  4,  9, 16])
>>> b = obj.my_method(np.arange(0, 5))
>>> a is b
False
>>> cached = MyClass(cached=True)
>>> a = cached.my_method(np.arange(0, 5))
>>> a
array([ 0,  1,  4,  9, 16])
>>> b = cached.my_method(np.arange(0, 5))
>>> a is b
True
>>> c = cached.my_method(np.arange(10, 15))
>>> c
array([ 0,  1,  4,  9, 16])
>>> a is c
True

hard_reset()

Delete the cached values.

class geoist.inversion.base.MultiObjective(*args)

基类：geoist.inversion.base.OptimizerMixin, geoist.inversion.base.OperatorMixin

An objective (goal) function with more than one component.

This class is a linear combination of other goal functions (like Misfit and regularization classes).

It is automatically created by adding two goal functions that have the OperatorMixin as a base class.

Alternatively, you can create a MultiObjetive by passing the other goals function instances as arguments to the constructor.

The MultiObjetive behaves like any other goal function object. It has fit and config methods and can be added and multiplied by a scalar with the same effects.

Indexing a MultiObjetive will iterate over the component goal functions.

Examples:

To show how this class is generated and works, let's create a simple class that subclasses OperatorMixin.

>>> class MyGoal(OperatorMixin):
...     def __init__(self, name, nparams, islinear):
...         self.name = name
...         self.islinear = islinear
...         self.nparams = nparams
...     def value(self, p):
...         return 1
...     def gradient(self, p):
...         return 2
...     def hessian(self, p):
...         return 3
>>> a = MyGoal('A', 10, True)
>>> b = MyGoal('B', 10, True)
>>> c = a + b
>>> type(c)
<class 'geoist.inversion.base.MultiObjective'>
>>> c.size
2
>>> c.nparams
10
>>> c.islinear
True
>>> c[0].name
'A'
>>> c[1].name
'B'

Asking for the value, gradient, and Hessian of the MultiObjective will give me the sum of both components.

>>> c.value(None)
2
>>> c.gradient(None)
4
>>> c.hessian(None)
6

Multiplying the MultiObjective by a scalar will set the regularization parameter for the sum.

>>> d = 10*c
>>> d.value(None)
20
>>> d.gradient(None)
40
>>> d.hessian(None)
60

All components must have the same number of parameters. For the moment, MultiObjetive doesn't handle multiple parameter vector (one for each objective function).

>>> e = MyGoal("E", 20, True)
>>> a + e
Traceback (most recent call last):
  ...
AssertionError: Can't add goals with different number of parameters: 10, 20

The MultiObjective will automatically detect if the problem remains linear or not. For example, adding a non-linear problem to a linear one makes the sum non-linear.

>>> (a + b).islinear
True
>>> f = MyGoal('F', 10, False)
>>> (a + f).islinear
False
>>> (f + f).islinear
False

config(*args, **kwargs)

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• methodstring

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)

• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear

• newton

• levmarq

• steepest

• acor

fit()

Solve for the parameter vector that minimizes this objective function.

Uses the optimization method and parameters defined using the config method.

The estimated parameter vector can be accessed through the p_ attribute. A (possibly) formatted version (converted to a more manageable type) of the estimate can be accessed through the property estimate_.

fmt_estimate(p)

Format the current estimated parameter vector into a more useful form.

Will call the fmt_estimate method of the first component goal function (the first term in the addition that created this object).

gradient(p)

Return the gradient of the multi-objective function.

This will be the sum of all goal functions that make up this multi-objective.

Parameters:

• p1d-array

  The parameter vector.

Returns:

• result1d-array

  The sum of the gradients of the components.

hessian(p)

Return the hessian of the multi-objective function.

This will be the sum of all goal functions that make up this multi-objective.

Parameters:

• p1d-array

  The parameter vector.

Returns:

• result2d-array

  The sum of the hessians of the components.

value(p)

Return the value of the multi-objective function.

This will be the sum of all goal functions that make up this multi-objective.

Parameters:

• p1d-array

  The parameter vector.

Returns:

• resultscalar (float, int, etc)

  The sum of the values of the components.

class geoist.inversion.base.OperatorMixin

基类：object

Implements the operators + and * for the goal functions classes.

This class is not meant to be used on its own. Use it as a parent to give the child class the + and * operators.

Used in Misfit and the regularization classes in geoist.inversion.regularization.

注解

Performing A + B produces a MultiObjetive with copies of A and B.

注解

Performing scalar*A produces a copy of A with scalar set as the regul_param attribute.

copy(deep=False)

Make a copy of me.

property regul_param

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

class geoist.inversion.base.OptimizerMixin

基类：object

Defines fit and config methods plus all the optimization methods.

This class is not meant to be used on its own. Use it as a parent to give the child class the methods it implements.

Used in Misfit and geoist.inversion.base.MultiObjetive.

The config method is used to configure the optimization method that will be used.

The fit method runs the optimization method configured and stores the computed parameter vector in the p_ attribute.

Some stats about the optimization process are stored in the stats_ attribute as a dictionary.

The minimum requirement for a class to inherit from OptimizerMixin is that it must define at least a value method.

config(method, **kwargs)

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• methodstring

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)

• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear

• newton

• levmarq

• steepest

• acor

property estimate_

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a geoist.mesher object.

fit()

Solve for the parameter vector that minimizes this objective function.

Uses the optimization method and parameters defined using the config method.

The estimated parameter vector can be accessed through the p_ attribute. A (possibly) formatted version (converted to a more manageable type) of the estimate can be accessed through the property estimate_.

fmt_estimate(p)

Called when accessing the property estimate_.

Use this to convert the parameter vector (p) to a more useful form, like a geometric object, etc.

Parameters:

• p1d-array

  The parameter vector.

Returns:

• formatted

  Pretty much anything you want.

geoist.inversion.geometry module

Defines geometric primitives like prisms, spheres, etc.

class geoist.inversion.geometry.GeometricElement(props)

基类：object

Base class for all geometric elements.

addprop(prop, value)

Add a physical property to this geometric element.

If it already has the property, the given value will overwrite the existing one.

Parameters:

• propstr

  Name of the physical property.

• valuefloat

  The value of this physical property.

copy()

Return a deep copy of the current instance.

class geoist.inversion.geometry.Polygon(vertices, props=None)

基类：geoist.inversion.geometry.GeometricElement

A polygon object (2D).

注解

Most applications require the vertices to be clockwise!

Parameters:

• verticeslist of lists

  List of [x, y] pairs with the coordinates of the vertices.

• propsdict

  Physical properties assigned to the polygon. Ex: props={'density':10, 'susceptibility':10000}

Examples:

>>> poly = Polygon([[0, 0], [1, 4], [2, 5]], {'density': 500})
>>> poly.props
{'density': 500}
>>> poly.nverts
3
>>> poly.vertices
array([[0, 0],
       [1, 4],
       [2, 5]])
>>> poly.x
array([0, 1, 2])
>>> poly.y
array([0, 4, 5])

property nverts

property vertices

property x

property y

class geoist.inversion.geometry.PolygonalPrism(vertices, z1, z2, props=None)

基类：geoist.inversion.geometry.GeometricElement

A 3D prism with polygonal cross-section.

注解

The coordinate system used is x -> North, y -> East and z -> Down

注解

vertices must be CLOCKWISE or will give inverse result.

Parameters:

• verticeslist of lists

  Coordinates of the vertices. A list of [x, y] pairs.

• z1, z2float

  Top and bottom of the prism

• propsdict

  Physical properties assigned to the prism. Ex: props={'density':10, 'magnetization':10000}

实际案例

>>> verts = [[1, 1], [1, 2], [2, 2], [2, 1]]
>>> p = PolygonalPrism(verts, 0, 3, props={'temperature':25})
>>> p.props['temperature']
25
>>> print p.x
[ 1.  1.  2.  2.]
>>> print p.y
[ 1.  2.  2.  1.]
>>> print p.z1, p.z2
0.0 3.0
>>> p.addprop('density', 2670)
>>> print p.props['density']
2670

topolygon()

Get the polygon describing the prism viewed from above.

Returns:

• polygonfatiando.mesher.Polygon

  The polygon

示例

>>> verts = [[1, 1], [1, 2], [2, 2], [2, 1]]
>>> p = PolygonalPrism(verts, 0, 100)
>>> poly = p.topolygon()
>>> print poly.x
[ 1.  1.  2.  2.]
>>> print poly.y
[ 1.  2.  2.  1.]

class geoist.inversion.geometry.Prism(x1, x2, y1, y2, z1, z2, props=None)

基类：geoist.inversion.geometry.GeometricElement

A 3D right rectangular prism.

注解

The coordinate system used is x -> North, y -> East and z -> Down

Parameters:

• x1, x2float

  South and north borders of the prism

• y1, y2float

  West and east borders of the prism

• z1, z2float

  Top and bottom of the prism

• propsdict

  Physical properties assigned to the prism. Ex: props={'density':10, 'magnetization':10000}

实际案例

>>> from geoist.inversion import Prism
>>> p = Prism(1, 2, 3, 4, 5, 6, {'density':200})
>>> p.props['density']
200
>>> print p.get_bounds()
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
>>> print p
x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | density:200
>>> p = Prism(1, 2, 3, 4, 5, 6)
>>> print p
x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6
>>> p.addprop('density', 2670)
>>> print p
x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | density:2670

center()

Return the coordinates of the center of the prism.

Returns:

• coordslist = [xc, yc, zc]

  Coordinates of the center

示例

>>> prism = Prism(1, 2, 1, 3, 0, 2)
>>> print prism.center()
[ 1.5  2.   1. ]

get_bounds()

Get the bounding box of the prism (i.e., the borders of the prism).

Returns:

• boundslist

  [x1, x2, y1, y2, z1, z2], the bounds of the prism

实际案例

>>> p = Prism(1, 2, 3, 4, 5, 6)
>>> print p.get_bounds()
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]

class geoist.inversion.geometry.Sphere(x, y, z, radius, props=None)

基类：geoist.inversion.geometry.GeometricElement

A sphere.

注解

The coordinate system used is x -> North, y -> East and z -> Down

Parameters:

• x, y, zfloat

  The coordinates of the center of the sphere

• radiusfloat

  The radius of the sphere

• propsdict

  Physical properties assigned to the prism. Ex: props={'density':10, 'magnetization':10000}

实际案例

>>> s = Sphere(1, 2, 3, 10, {'magnetization':200})
>>> s.props['magnetization']
200
>>> s.addprop('density', 20)
>>> print s.props['density']
20
>>> print s
x:1 | y:2 | z:3 | radius:10 | density:20 | magnetization:200
>>> s = Sphere(1, 2, 3, 4)
>>> print s
x:1 | y:2 | z:3 | radius:4
>>> s.addprop('density', 2670)
>>> print s
x:1 | y:2 | z:3 | radius:4 | density:2670

class geoist.inversion.geometry.Square(bounds, props=None)

基类：geoist.inversion.geometry.Polygon

A square object (2D).

Parameters:

• boundslist = [x1, x2, y1, y2]

  Coordinates of the top right and bottom left corners of the square

• propsdict

  Physical properties assigned to the square. Ex: props={'density':10, 'slowness':10000}

Example:

>>> sq = Square([0, 1, 2, 4], {'density': 750})
>>> sq.bounds
[0, 1, 2, 4]
>>> sq.x1
0
>>> sq.x2
1
>>> sq.props
{'density': 750}
>>> sq.addprop('magnetization', 100)
>>> sq.props['magnetization']
100

A square can be used as a Polygon:

>>> sq.vertices
array([[0, 2],
       [1, 2],
       [1, 4],
       [0, 4]])
>>> sq.x
array([0, 1, 1, 0])
>>> sq.y
array([2, 2, 4, 4])
>>> sq.nverts
4

property bounds

The x, y boundaries of the square as [xmin, xmax, ymin, ymax]

property vertices

The vertices of the square.

class geoist.inversion.geometry.Tesseroid(w, e, s, n, top, bottom, props=None)

基类：geoist.inversion.geometry.GeometricElement

A tesseroid (spherical prism).

Parameters:

• w, efloat

  West and east borders of the tesseroid in decimal degrees

• s, nfloat

  South and north borders of the tesseroid in decimal degrees

• top, bottomfloat

  Bottom and top of the tesseroid with respect to the mean earth radius in meters. Ex: if the top is 100 meters above the mean earth radius, top=100, if 100 meters below top=-100.

• propsdict

  Physical properties assigned to the tesseroid. Ex: props={'density':10, 'magnetization':10000}

实际案例

>>> from fatiando.mesher import Tesseroid
>>> t = Tesseroid(1, 2, 3, 4, 6, 5, {'density':200})
>>> t.props['density']
200
>>> print t.get_bounds()
[1.0, 2.0, 3.0, 4.0, 6.0, 5.0]
>>> print t
w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5 | density:200
>>> t = Tesseroid(1, 2, 3, 4, 6, 5)
>>> print t
w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5
>>> t.addprop('density', 2670)
>>> print t
w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5 | density:2670

get_bounds()

Get the bounding box of the tesseroid (i.e., the borders).

Returns:

• boundslist

  [w, e, s, n, top, bottom], the bounds of the tesseroid

实际案例

>>> t = Tesseroid(1, 2, 3, 4, 6, 5)
>>> print t.get_bounds()
[1.0, 2.0, 3.0, 4.0, 6.0, 5.0]

half(lon=True, lat=True, r=True)

Divide the tesseroid in 2 halfs for each dimension (total 8)

The smaller tesseroids will share the large one's props.

Parameters:

• lon, lat, rTrue or False

  Dimensions along which the tesseroid will be split in half.

Returns:

• tesseroidslist

  A list of maximum 8 tesseroids that make up the larger one.

Examples:

>>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2})
>>> split = tess.half()
>>> print len(split)
8
>>> for t in split:
...     print t
w:-10 | e:0 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:-10 | e:0 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:-10 | e:0 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:-10 | e:0 | s:0 | n:20 | top:0 | bottom:-20 | density:2
w:0 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:0 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:0 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:0 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2
>>> tess = Tesseroid(-15, 15, -20, 20, 0, -40)
>>> split = tess.half(lat=False)
>>> print len(split)
4
>>> for t in split:
...     print t
w:-15 | e:0 | s:-20 | n:20 | top:-20 | bottom:-40
w:-15 | e:0 | s:-20 | n:20 | top:0 | bottom:-20
w:0 | e:15 | s:-20 | n:20 | top:-20 | bottom:-40
w:0 | e:15 | s:-20 | n:20 | top:0 | bottom:-20

split(nlon, nlat, nh)

Split the tesseroid into smaller ones.

The smaller tesseroids will share the large one's props.

Parameters:

• nlon, nlat, nhint

  The number of sections to split in the longitudinal, latitudinal, and vertical dimensions

Returns:

• tesseroidslist

  A list of nlon*nlat*nh tesseroids that make up the larger one.

Examples:

>>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2})
>>> split = tess.split(1, 2, 2)
>>> print len(split)
4
>>> for t in split:
...     print t
w:-10 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:-10 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:-10 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:-10 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2
>>> tess = Tesseroid(-15, 15, -20, 20, 0, -40)
>>> split = tess.split(3, 1, 1)
>>> print len(split)
3
>>> for t in split:
...     print t
w:-15 | e:-5 | s:-20 | n:20 | top:0 | bottom:-40
w:-5 | e:5 | s:-20 | n:20 | top:0 | bottom:-40
w:5 | e:15 | s:-20 | n:20 | top:0 | bottom:-40

geoist.inversion.grawavelet module

Created on Sun Mar 10 20:49:28 2019

@author: chens

geoist.inversion.grawavelet.walshMatrix(n, normalized=False)

geoist.inversion.harrwavelet module

Created on Sun Mar 10 20:49:28 2019

@author: chens

geoist.inversion.harrwavelet.haarMatrix(n, normalized=False)

geoist.inversion.harrwavelet.walshMatrix(n, normalized=False)

geoist.inversion.hyper_param module

Classes for hyper parameter estimation (like the regularizing parameter).

These classes copy the interface of the standard inversion classes based on Misfit (i.e., solver.config(...).fit().estimate_). When their fit method is called, they perform many runs of the inversion and try to select the optimal values for the hyper parameters. The class will then behave as the solver that yields the best estimate (e.g., solver[0].predicted()).

Available classes:

• LCurve: Estimate the regularizing parameter using an L-curve analysis.

---

class geoist.inversion.hyper_param.LCurve(datamisfit, regul, regul_params, loglog=True, njobs=1)

基类：geoist.inversion.base.OptimizerMixin

Use the L-curve criterion to estimate the regularization parameter.

Runs the inversion using several specified regularization parameters. The best value is the one that falls on the corner of the log-log plot of the data misfit vs regularizing function. This point is automatically found using the triangle method of Castellanos et al. (2002).

This class behaves as Misfit. To use it, simply call fit and optionally config. The estimate will be stored in estimate_ and p_. The estimated regularization parameter will be stored in regul_param_.

Parameters:

• datamisfitMisfit

  The data misfit instance for the inverse problem. Can be a sum of other misfits.

• regulA class from geoist.inversion.regularization

  The regularizing function.

• regul_paramslist

  The values of the regularization parameter that will be tested.

• loglogTrue or False

  If True, will use a log-log scale for the L-curve (recommended).

• jobsNone or int

  If not None, will use jobs processes to calculate the L-curve.

References:

Castellanos, J. L., S. Gomez, and V. Guerra (2002), The triangle method for finding the corner of the L-curve, Applied Numerical Mathematics, 43(4), 359-373, doi:10.1016/S0168-9274(01)00179-9.

Examples:

We'll use the L-curve to estimate the best regularization parameter for a smooth inversion using geoist.seismic.srtomo.

First, we'll setup some synthetic data:

>>> import numpy
>>> from geoist.mesher import SquareMesh
>>> from geoist.seismic import ttime2d, srtomo
>>> from geoist.inversion import Smoothness2D, LCurve
>>> from geoist import utils, gridder
>>> area = (0, 2, 0, 2)
>>> shape = (10, 10)
>>> model = SquareMesh(area, shape)
>>> vp = 4*numpy.ones(shape)
>>> vp[3:7,3:7] = 10
>>> vp
array([[  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.]])
>>> model.addprop('vp', vp.ravel())
>>> src_loc_x, src_loc_y = gridder.scatter(area, 30, seed=0)
>>> src_loc = numpy.transpose([src_loc_x, src_loc_y])
>>> rec_loc_x, rec_loc_y = gridder.circular_scatter(area, 20,
...                                                 random=True, seed=0)
>>> rec_loc = numpy.transpose([rec_loc_x, rec_loc_y])
>>> srcs = [src for src in src_loc for _ in rec_loc]
>>> recs = [rec for _ in src_loc for rec in rec_loc]
>>> tts = ttime2d.straight(model, 'vp', srcs, recs)
>>> tts = utils.contaminate(tts, 0.01, percent=True, seed=0)

Now we can setup a tomography by creating the necessary data misfit (SRTomo) and regularization (Smoothness2D) objects. We'll normalize the data misfit by the number of data points to make the scale of the regularization parameter more tractable.

>>> mesh = SquareMesh(area, shape)
>>> datamisfit = (1./tts.size)*srtomo.SRTomo(tts, srcs, recs, mesh)
>>> regul = Smoothness2D(mesh.shape)

The tomography solver will be the LCurve solver. It works by calling fit() and accessing estimate_, exactly like any other solver:

>>> regul_params = [10**i for i in range(-10, -2, 1)]
>>> tomo = LCurve(datamisfit, regul, regul_params)
>>> _ = tomo.fit()
>>> print(numpy.array_repr(tomo.estimate_.reshape(shape), precision=0))
array([[  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,   9.,  11.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  11.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,  10.,  11.,   9.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.]])

When fit() is called, the LCurve will run the inversion for each value of the regularization parameter, build an l-curve, and find the best solution (i.e., the corner value of the l-curve).

The LCurve object behaves like a normal multi-objective function. In fact, it will try to mirror the objective function that resulted in the best solution. You can index it to access the data misfit and regularization parts. For example, to get the residuals vector or the predicted data:

>>> predicted = tomo[0].predicted()
>>> residuals = tomo[0].residuals()
>>> print '%.4f %.4f' % (residuals.mean(), residuals.std())
-0.0000 0.0047

The estimated regularization parameter is stored in regul_param_:

>>> tomo.regul_param_
1e-05

You can run the l-curve analysis in parallel by specifying the njobs argument. This will spread the computations over njobs number of processes and give some speedup over running sequentially. Note that you should not enable any kind of multi-processes parallelism on the data misfit class. It is often better to run each inversion sequentially and run many of them in parallel. Note that you'll enough memory to run multiple inversions at the same time, so this is not suited for large, memory hungry inversions.

>>> par_tomo = LCurve(datamisfit, regul, regul_params, njobs=2)
>>> _ = par_tomo.fit()  # Will you 2 processes to run inversions
>>> par_tomo.regul_param_
1e-05
>>> print(numpy.array_repr(par_tomo.estimate_.reshape(shape), precision=0))
array([[  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,   9.,  11.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  11.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,  10.,  11.,   9.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.]])

LCurve also has a config method to configure the optimization process for non-linear problems, for example:

>>> initial = numpy.ones(mesh.size)
>>> _ = tomo.config('newton', initial=initial, tol=0.2).fit()
>>> tomo.regul_param_
1e-05
>>> print(numpy.array_repr(tomo.estimate_.reshape(shape), precision=0))
array([[  4.,   4.,   3.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  12.,   9.,  11.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,  11.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  10.,  10.,  10.,  10.,   4.,   4.,   4.],
       [  4.,   4.,   4.,  11.,  10.,  11.,   9.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   5.,   4.,   4.,   4.],
       [  4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.,   4.]])

You can view the optimization information for the run corresponding to the best estimate using the stats_ attribute:

>>> list(sorted(tomo.stats_))
['iterations', 'method', 'objective']
>>> tomo.stats_['method']
"Newton's method"
>>> tomo.stats_['iterations']
2

fit()

Solve for the parameter vector and optimum regularization parameter.

Combines the data-misfit and regularization solvers using the range of regularization parameters provided and calls fit and config on each.

The p_ and estimate_ attributes correspond to the combination that falls in the corner of the L-curve.

The regularization parameter for this corner point if stored in the regul_param_ attribute.

Returns:

• self

fmt_estimate(p)

Return the estimate_ attribute of the optimal solution.

property objective_

The objective function corresponding to the best estimate.

property p_

The estimated parameter vector obtained from the best regularization parameter.

plot_lcurve(ax=None, guides=True)

Make a plot of the data-misfit x regularization values.

The estimated corner value is shown as a blue triangle.

Parameters:

• axmatplotlib Axes

  If not None, will plot the curve on this Axes instance.

• guidesTrue or False

  Plot vertical and horizontal lines across the corner value.

property regul_param_

The regularization parameter corresponding to the best estimate.

select_corner()

Select the corner value of the L-curve formed inversion results.

This is performed automatically after calling the fit method. You can run this method separately after fit has been called to tweak the results.

You can access the estimated values by:

• The p_ and estimate_ attributes will hold the estimated parameter vector and formatted estimate, respective, corresponding to the corner value.

• The regul_param_ attribute holds the value of the regularization parameter corresponding to the corner value.

• The corner_ attribute will hold the index of the corner value in the list of computed solutions.

Uses the Triangle method of Castellanos et al. (2002).

References:

Castellanos, J. L., S. Gomez, and V. Guerra (2002), The triangle method for finding the corner of the L-curve, Applied Numerical Mathematics, 43(4), 359-373, doi:10.1016/S0168-9274(01)00179-9.

property stats_

The optimization information for the best solution found.

geoist.inversion.mesh module

Meshes (collections) of geometric objects.

Meshes behave like lists/arrays of geometric objects (they are iterables).

class geoist.inversion.mesh.PointGrid(area, z, shape, props=None)

基类：object

A regular grid of 3D point sources (spheres of unit volume).

Use this as a 1D list of Sphere.

Grid points are ordered like a C matrix, first each row in a column, then change columns. In this case, the x direction (North-South) are the rows and y (East-West) are the columns.

Parameters:

• arealist = [x1, x2, y1, y2]

  The area where the grid will be spread out

• zfloat or 1d-array

  The z coordinates of each point in the grid (remember, z is positive downward).

• shapetuple = (nx, ny)

  The number of points in the x and y directions

• propsdict

  Physical properties of each point in the grid. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property for each point in the grid.

Examples:

>>> g = PointGrid([0, 10, 2, 6], 200, (2, 3))
>>> g.shape
(2, 3)
>>> g.size
6
>>> g[0].center
array([   0.,    2.,  200.])
>>> g[-1].center
array([  10.,    6.,  200.])
>>> for p in g:
...     p.center
array([   0.,    2.,  200.])
array([   0.,    4.,  200.])
array([   0.,    6.,  200.])
array([  10.,    2.,  200.])
array([  10.,    4.,  200.])
array([  10.,    6.,  200.])
>>> g.x.reshape(g.shape)
array([[  0.,   0.,   0.],
       [ 10.,  10.,  10.]])
>>> g.y.reshape(g.shape)
array([[ 2.,  4.,  6.],
       [ 2.,  4.,  6.]])
>>> g.dx, g.dy
(10.0, 2.0)

addprop(prop, values)

Add physical property values to the points in the grid.

Different physical properties of the grid are stored in a dictionary.

Parameters:

• propstr

  Name of the physical property.

• valueslist or array

  Value of this physical property in each point of the grid

copy()

Return a deep copy of the current instance.

split(shape)

Divide the grid into subgrids.

注解

Remember that x is the North-South direction and y is East-West.

Parameters:

• shapetuple = (nx, ny)

  Number of subgrids along the x and y directions, respectively.

Returns:

• subgridslist

  List of PointGrid

Examples:

>>> import numpy as np
>>> z = np.linspace(0, 1100, 12)
>>> g = PointGrid((0, 3, 0, 2), z, (4, 3))
>>> g.addprop('bla', [1,   2,  3,
...                   4,   5,  6,
...                   7,   8,  9,
...                   10, 11, 12])
>>> grids = g.split((2, 3))
>>> for s in grids:
...     s.props['bla']
array([1, 4])
array([2, 5])
array([3, 6])
array([ 7, 10])
array([ 8, 11])
array([ 9, 12])
>>> for s in grids:
...     s.x
array([ 0.,  1.])
array([ 0.,  1.])
array([ 0.,  1.])
array([ 2.,  3.])
array([ 2.,  3.])
array([ 2.,  3.])
>>> for s in grids:
...     s.y
array([ 0.,  0.])
array([ 1.,  1.])
array([ 2.,  2.])
array([ 0.,  0.])
array([ 1.,  1.])
array([ 2.,  2.])
>>> for s in grids:
...     s.z
array([   0.,  300.])
array([ 100.,  400.])
array([ 200.,  500.])
array([ 600.,  900.])
array([  700.,  1000.])
array([  800.,  1100.])

class geoist.inversion.mesh.PrismMesh(bounds, shape, props=None)

基类：object

A 3D regular mesh of right rectangular prisms.

Prisms are ordered as follows: first layers (z coordinate), then EW rows (y) and finaly x coordinate (NS).

注解

Remember that the coordinate system is x->North, y->East and z->Down

Ex: in a mesh with shape (3,3,3) the 15th element (index 14) has z index 1 (second layer), y index 1 (second row), and x index 2 (third element in the column).

PrismMesh can used as list of prisms. It acts as an iteratior (so you can loop over prisms). It also has a __getitem__ method to access individual elements in the mesh. In practice, PrismMesh should be able to be passed to any function that asks for a list of prisms, like geoist.pfm.prism.gz.

To make the mesh incorporate a topography, use carvetopo

Parameters:

• boundslist = [xmin, xmax, ymin, ymax, zmin, zmax]

  Boundaries of the mesh.

• shapetuple = (nz, ny, nx)

  Number of prisms in the x, y, and z directions.

• propsdict

  Physical properties of each prism in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each prism of the mesh.

实际案例

>>> from geoist.inversion import PrismMesh
>>> mesh = PrismMesh((0, 1, 0, 2, 0, 3), (1, 2, 2))
>>> for p in mesh:
...     print p
x1:0 | x2:0.5 | y1:0 | y2:1 | z1:0 | z2:3
x1:0.5 | x2:1 | y1:0 | y2:1 | z1:0 | z2:3
x1:0 | x2:0.5 | y1:1 | y2:2 | z1:0 | z2:3
x1:0.5 | x2:1 | y1:1 | y2:2 | z1:0 | z2:3
>>> print mesh[0]
x1:0 | x2:0.5 | y1:0 | y2:1 | z1:0 | z2:3
>>> print mesh[-1]
x1:0.5 | x2:1 | y1:1 | y2:2 | z1:0 | z2:3

One with physical properties:

>>> props = {'density':[2670.0, 1000.0]}
>>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2), props=props)
>>> for p in mesh:
...     print p
x1:0 | x2:1 | y1:0 | y2:4 | z1:0 | z2:3 | density:2670
x1:1 | x2:2 | y1:0 | y2:4 | z1:0 | z2:3 | density:1000

or equivalently:

>>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2))
>>> mesh.addprop('density', [200, -1000.0])
>>> for p in mesh:
...     print p
x1:0 | x2:1 | y1:0 | y2:4 | z1:0 | z2:3 | density:200
x1:1 | x2:2 | y1:0 | y2:4 | z1:0 | z2:3 | density:-1000

You can use get_xs (and similar methods for y and z) to get the x coordinates of the prisms in the mesh:

>>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2))
>>> print mesh.get_xs()
[ 0.  1.  2.]
>>> print mesh.get_ys()
[ 0.  4.]
>>> print mesh.get_zs()
[ 0.  3.]

The shape of the mesh must be integer!

>>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2.5))
Traceback (most recent call last):
    ...
AttributeError: Invalid mesh shape (1, 1, 2.5). shape must be integers

addprop(prop, values)

Add physical property values to the cells in the mesh.

Different physical properties of the mesh are stored in a dictionary.

Parameters:

• propstr

  Name of the physical property.

• valueslist or array

  Value of this physical property in each prism of the mesh. For the ordering of prisms in the mesh see PrismMesh

carvetopo(x, y, height, below=False)

Mask (remove) prisms from the mesh that are above the topography.

Accessing the ith prism will return None if it was masked (above the topography). Also mask prisms outside of the topography grid provided. The topography height information does not need to be on a regular grid, it will be interpolated.

Parameters:

• x, ylists

  x and y coordinates of the grid points

• heightlist or array

  Array with the height of the topography

• belowboolean

  Will mask prisms below the input surface if set to True.

celltype

geoist.inversion.geometry.Prism 的别名

copy()

Return a deep copy of the current instance.

dump(meshfile, propfile, prop)

Dump the mesh to a file in the format required by UBC-GIF program MeshTools3D.

Parameters:

• meshfilestr or file

  Output file to save the mesh. Can be a file name or an open file.

• propfilestr or file

  Output file to save the physical properties prop. Can be a file name or an open file.

• propstr

  The name of the physical property in the mesh that will be saved to propfile.

注解

Uses -10000000 as the dummy value for plotting topography

实际案例

>>> from StringIO import StringIO
>>> meshfile = StringIO()
>>> densfile = StringIO()
>>> mesh = PrismMesh((0, 10, 0, 20, 0, 5), (1, 2, 2))
>>> mesh.addprop('density', [1, 2, 3, 4])
>>> mesh.dump(meshfile, densfile, 'density')
>>> print meshfile.getvalue().strip()
2 2 1
0 0 0
2*10
2*5
1*5
>>> print densfile.getvalue().strip()
1.0000
3.0000
2.0000
4.0000

get_layer(i)

Return the set of prisms corresponding to the ith layer of the mesh.

Parameters:

• iint

  The index of the layer

Returns:

• prismslist of Prism

  The prisms in the ith layer

Examples:

>>> mesh = PrismMesh((0, 2, 0, 2, 0, 2), (2, 2, 2))
>>> layer = mesh.get_layer(0)
>>> for p in layer:
...     print p
x1:0 | x2:1 | y1:0 | y2:1 | z1:0 | z2:1
x1:1 | x2:2 | y1:0 | y2:1 | z1:0 | z2:1
x1:0 | x2:1 | y1:1 | y2:2 | z1:0 | z2:1
x1:1 | x2:2 | y1:1 | y2:2 | z1:0 | z2:1
>>> layer = mesh.get_layer(1)
>>> for p in layer:
...     print p
x1:0 | x2:1 | y1:0 | y2:1 | z1:1 | z2:2
x1:1 | x2:2 | y1:0 | y2:1 | z1:1 | z2:2
x1:0 | x2:1 | y1:1 | y2:2 | z1:1 | z2:2
x1:1 | x2:2 | y1:1 | y2:2 | z1:1 | z2:2

get_xs()

Return an array with the x coordinates of the prisms in mesh.

get_ys()

Return an array with the y coordinates of the prisms in mesh.

get_zs()

Return an array with the z coordinates of the prisms in mesh.

layers()

Returns an iterator over the layers of the mesh.

Examples:

>>> mesh = PrismMesh((0, 2, 0, 2, 0, 2), (2, 2, 2))
>>> for layer in mesh.layers():
...     for p in layer:
...         print p
x1:0 | x2:1 | y1:0 | y2:1 | z1:0 | z2:1
x1:1 | x2:2 | y1:0 | y2:1 | z1:0 | z2:1
x1:0 | x2:1 | y1:1 | y2:2 | z1:0 | z2:1
x1:1 | x2:2 | y1:1 | y2:2 | z1:0 | z2:1
x1:0 | x2:1 | y1:0 | y2:1 | z1:1 | z2:2
x1:1 | x2:2 | y1:0 | y2:1 | z1:1 | z2:2
x1:0 | x2:1 | y1:1 | y2:2 | z1:1 | z2:2
x1:1 | x2:2 | y1:1 | y2:2 | z1:1 | z2:2

class geoist.inversion.mesh.PrismRelief(ref, dims, nodes)

基类：object

A 3D model of a relief (topography) using prisms.

Use to generate: * topographic model * basin model * Moho model * etc

PrismRelief can used as list of prisms. It acts as an iteratior (so you can loop over prisms). It also has a __getitem__ method to access individual elements in the mesh. In practice, PrismRelief should be able to be passed to any function that asks for a list of prisms, like geoist.pfm.prism.gz.

Parameters:

• reffloat

  Reference level. Prisms will have:

  • bottom on zref and top on z if z > zref;

  • bottom on z and top on zref otherwise.

• dimstuple = (dy, dx)

  Dimensions of the prisms in the y and x directions

• nodeslist of lists = [x, y, z]

  Coordinates of the center of the top face of each prism.x, y, and z are lists with the x, y and z coordinates on a regular grid.

addprop(prop, values)

Add physical property values to the prisms.

警告

If the z value of any point in the relief is below the reference level, its corresponding prism will have the physical property value with oposite sign than was assigned to it.

Parameters:

• propstr

  Name of the physical property.

• valueslist

  List or array with the value of this physical property in each prism of the relief.

copy()

Return a deep copy of the current instance.

class geoist.inversion.mesh.SquareMesh(bounds, shape, props=None)

基类：object

A 2D regular mesh of squares.

For all purposes, SquareMesh can be used as a list of Square. The order of the squares in the list is: x directions varies first, then y.

Parameters:

• boundslist = [x1, x2, y1, y2]

  Boundaries of the mesh

• shapetuple = (ny, nx)

  Number of squares in the y and x dimension, respectively

• propsdict

  Physical properties of each square in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each square of the mesh.

实际案例

>>> mesh = SquareMesh((0, 4, 0, 6), (2, 2))
>>> for s in mesh:
...     print s
x1:0 | x2:2 | y1:0 | y2:3
x1:2 | x2:4 | y1:0 | y2:3
x1:0 | x2:2 | y1:3 | y2:6
x1:2 | x2:4 | y1:3 | y2:6
>>> print mesh[1]
x1:2 | x2:4 | y1:0 | y2:3
>>> print mesh[-1]
x1:2 | x2:4 | y1:3 | y2:6

With physical properties:

>>> mesh = SquareMesh((0, 4, 0, 6), (2, 1), {'slowness':[3.4, 8.6]})
>>> for s in mesh:
...     print s
x1:0 | x2:4 | y1:0 | y2:3 | slowness:3.4
x1:0 | x2:4 | y1:3 | y2:6 | slowness:8.6

Or:

>>> mesh = SquareMesh((0, 4, 0, 6), (2, 1))
>>> mesh.addprop('slowness', [3.4, 8.6])
>>> for s in mesh:
...     print s
x1:0 | x2:4 | y1:0 | y2:3 | slowness:3.4
x1:0 | x2:4 | y1:3 | y2:6 | slowness:8.6

addprop(prop, values)

Add physical property values to the cells in the mesh.

Different physical properties of the mesh are stored in a dictionary.

Parameters:

• propstr

  Name of the physical property

• valueslist or array

  The value of this physical property in each square of the mesh. For the ordering of squares in the mesh see SquareMesh

copy()

Return a deep copy of the current instance.

get_xs()

Get a list of the x coordinates of the corners of the cells in the mesh.

If the mesh has nx cells, get_xs() will return nx + 1 values.

get_ys()

Get a list of the y coordinates of the corners of the cells in the mesh.

If the mesh has ny cells, get_ys() will return ny + 1 values.

class geoist.inversion.mesh.TesseroidMesh(bounds, shape, props=None)

基类：geoist.inversion.mesh.PrismMesh

A 3D regular mesh of tesseroids.

Tesseroids are ordered as follows: first layers (height coordinate), then N-S rows and finaly E-W.

Ex: in a mesh with shape (3,3,3) the 15th element (index 14) has height index 1 (second layer), y index 1 (second row), and x index 2 ( third element in the column).

This class can used as list of tesseroids. It acts as an iteratior (so you can loop over tesseroids). It also has a __getitem__ method to access individual elements in the mesh. In practice, it should be able to be passed to any function that asks for a list of tesseroids, like geoist.pfm.tesseroid.gz.

To make the mesh incorporate a topography, use carvetopo

Parameters:

• boundslist = [w, e, s, n, top, bottom]

  Boundaries of the mesh. w, e, s, n in degrees, top and bottom are heights (positive upward) and in meters.

• shapetuple = (nr, nlat, nlon)

  Number of tesseroids in the radial, latitude, and longitude directions.

• propsdict

  Physical properties of each tesseroid in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each tesseroid of the mesh.

实际案例

>>> from geoist.inversion import TesseroidMesh
>>> mesh = TesseroidMesh((0, 1, 0, 2, 3, 0), (1, 2, 2))
>>> for p in mesh:
...     print p
w:0 | e:0.5 | s:0 | n:1 | top:3 | bottom:0
w:0.5 | e:1 | s:0 | n:1 | top:3 | bottom:0
w:0 | e:0.5 | s:1 | n:2 | top:3 | bottom:0
w:0.5 | e:1 | s:1 | n:2 | top:3 | bottom:0
>>> print mesh[0]
w:0 | e:0.5 | s:0 | n:1 | top:3 | bottom:0
>>> print mesh[-1]
w:0.5 | e:1 | s:1 | n:2 | top:3 | bottom:0

One with physical properties:

>>> props = {'density':[2670.0, 1000.0]}
>>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2), props=props)
>>> for p in mesh:
...     print p
w:0 | e:1 | s:0 | n:4 | top:3 | bottom:0 | density:2670
w:1 | e:2 | s:0 | n:4 | top:3 | bottom:0 | density:1000

or equivalently:

>>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2))
>>> mesh.addprop('density', [200, -1000.0])
>>> for p in mesh:
...     print p
w:0 | e:1 | s:0 | n:4 | top:3 | bottom:0 | density:200
w:1 | e:2 | s:0 | n:4 | top:3 | bottom:0 | density:-1000

You can use get_xs (and similar methods for y and z) to get the x coordinates of the tesseroidss in the mesh:

>>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2))
>>> print mesh.get_xs()
[ 0.  1.  2.]
>>> print mesh.get_ys()
[ 0.  4.]
>>> print mesh.get_zs()
[ 3.  0.]

You can iterate over the layers of the mesh:

>>> mesh = TesseroidMesh((0, 2, 0, 2, 2, 0), (2, 2, 2))
>>> for layer in mesh.layers():
...     for p in layer:
...         print p
w:0 | e:1 | s:0 | n:1 | top:2 | bottom:1
w:1 | e:2 | s:0 | n:1 | top:2 | bottom:1
w:0 | e:1 | s:1 | n:2 | top:2 | bottom:1
w:1 | e:2 | s:1 | n:2 | top:2 | bottom:1
w:0 | e:1 | s:0 | n:1 | top:1 | bottom:0
w:1 | e:2 | s:0 | n:1 | top:1 | bottom:0
w:0 | e:1 | s:1 | n:2 | top:1 | bottom:0
w:1 | e:2 | s:1 | n:2 | top:1 | bottom:0

The shape of the mesh must be integer!

>>> mesh = TesseroidMesh((0, 2, 0, 4, 0, 3), (1, 1, 2.5))
Traceback (most recent call last):
    ...
AttributeError: Invalid mesh shape (1, 1, 2.5). shape must be integers

celltype

geoist.inversion.geometry.Tesseroid 的别名

geoist.inversion.misfit module

Defines base classes to represent a data-misfit functions (l2-norm, etc)

These classes can be used to implement parameter estimation problems (inversions). They automate most of the boiler plate required and provide direct access to ready-made optimization routines and regularization.

For now, only implements an l2-norm data misfit:

• Misfit: an l2-norm data-misfit function

See the documentation for geoist.inversion for examples of using Misfit.

---

class geoist.inversion.misfit.Misfit(data, nparams, islinear, cache=True)

基类：geoist.inversion.base.OptimizerMixin, geoist.inversion.base.OperatorMixin

An l2-norm data-misfit function.

This is a kind of objective function that measures the misfit between observed data \(\bar{d}^o\) and data predicted by a set of model parameters \(\bar{d} = \bar{f}(\bar{p})\).

The l2-norm data-misfit is defined as:

\[\phi (\bar{p}) = \bar{r}^T \bar{r}\]

where \(\bar{r} = \bar{d}^o - \bar{d}\) is the residual vector and \(N\) is the number of data.

When subclassing this class, you must implement the method:

• predicted(self, p): calculates the predicted data \(\bar{d}\) for a given parameter vector p

If you want to use any gradient-based solver (you probably do), you'll need to implement the method:

• jacobian(self, p): calculates the Jacobian matrix of \(\bar{f}(\bar{p})\) evaluated at p

If \(\bar{f}\) is linear, then the Jacobian will be cached in memory so that it is only calculated once when using the class multiple times. So solving the same problem with different methods or using an iterative method doesn't have the penalty of recalculating the Jacobian.

警告

When subclassing, be careful not to set the following attributes: data, nparams, islinear, nparams, ndata, and (most importantly) regul_param and _regularizing_parameter. This could mess with internal behavior and break things in unexpected ways.

Parameters:

• data1d-array

  The observed data vector \(\bar{d}^o\)

• nparamsint

  The number of parameters in parameter vector \(\bar{p}\)

• islinearTrue or False

  Whether \(\bar{f}\) is linear or not.

• cacheTrue

  Whether or not to cache the output of some methods to avoid recomputing matrices and vectors when passed the same input parameter vector.

copy(deep=False)

Make a copy of me together with all the cached methods.

gradient(p)

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

• p1d-array

  The parameter vector where the gradient is evaluated

Returns:

• gradient1d-array

  The gradient vector.

hessian(p)

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

• p1d-array

  The parameter vector where the Hessian is evaluated

Returns:

• hessian2d-array

  The Hessian matrix

abstract predicted(p=None)

Calculate the predicted data for a given parameter vector.

Parameters:

• p1d-array or None

  The parameter vector used to calculate the predicted data. If None, will use the current estimate stored in estimate_.

Returns:

• predicted1d-array or list of 1d-arrays

  The predicted data. If this is the sum of 1 or more Misfit instances, will return the predicted data from each of the summed misfits in the order of the sum.

residuals(p=None)

Calculate the residuals vector (observed - predicted data).

Parameters:

• p1d-array or None

  The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

• residuals1d-array or list of 1d-arrays

  The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

• weights1d-array or 2d-array or None

  Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

• p1d-array or None

  The parameter vector.

Returns:

• valuefloat

  The value of the misfit function.

geoist.inversion.optimization module

Methods to optimize a given objective function.

All solvers are Python iterators. This means that should be used in a for loop, like so:

solver = newton(hess_func, grad_func, value_func, initial)
for i, p, stats in solver:
    ... do something or 'continue' to step through the iterations ...
    # 'p' is the current estimate for the parameter vector at the 'i'th
    # iteration.
    # 'stats' is a dictionary with some information about the optimization
    # process so far (number of attempted steps, value of objective
    # function per step, total number of iterations so far, etc).
# At the end, 'p' is the final estimate and 'stats' will contain the
# statistics for the whole iteration process.

Gradient descent

• linear: Solver for a linear problem

• newton: Newton's method

• levmarq: Levemberg-Marquardt algorithm

• steepest: Steepest Descent method

Heuristic methods

• acor: ACO-R: Ant Colony Optimization for Continuous Domains (Socha and Dorigo, 2008)

References

Socha, K., and M. Dorigo (2008), Ant colony optimization for continuous domains, European Journal of Operational Research, 185(3), 1155-1173, doi:10.1016/j.ejor.2006.06.046.

---

geoist.inversion.optimization.acor(value, bounds, nparams, nants=None, archive_size=None, maxit=1000, diverse=0.5, evap=0.85, seed=None)

Minimize the objective function using ACO-R.

ACO-R stands for Ant Colony Optimization for Continuous Domains (Socha and Dorigo, 2008).

Parameters:

• valuefunction

  Returns the value of the objective function at a given parameter vector

• boundslist

  The bounds of the search space. If only two values are given, will interpret as the minimum and maximum, respectively, for all parameters. Alternatively, you can given a minimum and maximum for each parameter, e.g., for a problem with 3 parameters you could give bounds = [min1, max1, min2, max2, min3, max3].

• nparamsint

  The number of parameters that the objective function takes.

• nantsint

  The number of ants to use in the search. Defaults to the number of parameters.

• archive_sizeint

  The number of solutions to keep in the solution archive. Defaults to 10 x nants

• maxitint

  The number of iterations to run.

• diversefloat

  Scalar from 0 to 1, non-inclusive, that controls how much better solutions are favored when constructing new ones.

• evapfloat

  The pheromone evaporation rate (evap > 0). Controls how spread out the search is.

• seedNone or int

  Seed for the random number generator.

Yields:

• i, estimate, stats:

  • iint

    The current iteration number

  • estimate1d-array

    The current best estimated parameter vector

  • statsdict

    Statistics about the optimization so far. Keys:

    • methodstf

      The name of the optimization algorithm

    • iterationsint

      The total number of iterations so far

    • objectivelist

      Value of the objective function corresponding to the best estimate per iteration.

geoist.inversion.optimization.levmarq(hessian, gradient, value, initial, maxit=30, maxsteps=20, lamb=10, dlamb=2, tol=1e-05, precondition=True)

Minimize an objective function using the Levemberg-Marquardt algorithm.

Parameters:

• hessianfunction

  A function that returns the Hessian matrix of the objective function when given a parameter vector.

• gradientfunction

  A function that returns the gradient vector of the objective function when given a parameter vector.

• valuefunction

  A function that returns the value of the objective function evaluated at a given parameter vector.

• initial1d-array

  The initial estimate for the gradient descent.

• maxitint

  The maximum number of iterations allowed.

• maxstepsint

  The maximum number of times to try to take a step before giving up.

• lambfloat

  Initial amount of step regularization. The larger this is, the more the algorithm will resemble Steepest Descent in the initial iterations.

• dlambfloat

  Factor by which lamb is divided or multiplied when taking steps.

• tolfloat

  The convergence criterion. The lower it is, the more steps are permitted.

• preconditionTrue or False

  If True, will use Jacobi preconditioning.

Yields:

• i, estimate, stats:

  • iint

    The current iteration number

  • estimate1d-array

    The current estimated parameter vector

  • statsdict

    Statistics about the optimization so far. Keys:

    • methodstr

      The name of the optimization method

    • iterationsint

      The total number of iterations so far

    • objectivelist

      Value of the objective function per iteration. First value corresponds to the inital estimate

    • step_attemptslist

      Number of attempts at taking a step per iteration. First number is zero, reflecting the initial estimate.

geoist.inversion.optimization.linear(hessian, gradient, precondition=True)

Find the parameter vector that minimizes a linear objective function.

The parameter vector \(\bar{p}\) that minimizes this objective function \(\phi\) is the one that solves the linear system

\[\bar{\bar{H}} \bar{p} = -\bar{g}\]

where \(\bar{\bar{H}}\) is the Hessian matrix of \(\phi\) and \(\bar{g}\) is the gradient vector of \(\phi\).

Parameters:

• hessian2d-array

  The Hessian matrix of the objective function.

• gradient1d-array

  The gradient vector of the objective function.

• preconditionTrue or False

  If True, will use Jacobi preconditioning.

Yields:

• i, estimate, stats:

  • iint

    The current iteration number

  • estimate1d-array

    The current estimated parameter vector

  • statsdict

    Statistics about the optimization so far

Linear solvers have only a single step, so i will be 0 and stats will only have the method name.

geoist.inversion.optimization.newton(hessian, gradient, value, initial, maxit=30, tol=1e-05, precondition=True)

Minimize an objective function using Newton's method.

Newton's method searches for the minimum of an objective function \(\phi(\bar{p})\) by successively incrementing the initial estimate. The increment is the solution of the linear system

\[\bar{\bar{H}}(\bar{p}^k) \bar{\Delta p}^k = -\bar{g}(\bar{p}^k)\]

where \(\bar{\bar{H}}\) is the Hessian matrix of \(\phi\) and \(\bar{g}\) is the gradient vector of \(\phi\). Both are evaluated at the previous estimate \(\bar{p}^k\).

Parameters:

• hessianfunction

  A function that returns the Hessian matrix of the objective function when given a parameter vector.

• gradientfunction

  A function that returns the gradient vector of the objective function when given a parameter vector.

• valuefunction

  A function that returns the value of the objective function evaluated at a given parameter vector.

• initial1d-array

  The initial estimate for the gradient descent.

• maxitint

  The maximum number of iterations allowed.

• tolfloat

  The convergence criterion. The lower it is, the more steps are permitted.

• preconditionTrue or False

  If True, will use Jacobi preconditioning.

Returns:

Yields:

• i, estimate, stats:

  • iint

    The current iteration number

  • estimate1d-array

    The current estimated parameter vector

  • statsdict

    Statistics about the optimization so far. Keys:

    • methodstr

      The name of the optimization method

    • iterationsint

      The total number of iterations so far

    • objectivelist

      Value of the objective function per iteration. First value corresponds to the inital estimate

geoist.inversion.optimization.steepest(gradient, value, initial, maxit=1000, linesearch=True, maxsteps=30, beta=0.1, tol=1e-05)

Minimize an objective function using the Steepest Descent method.

The increment to the initial estimate of the parameter vector \(\bar{p}\) is calculated by (Kelley, 1999)

\[\Delta\bar{p} = -\lambda\bar{g}\]

where \(\lambda\) is the step size and \(\bar{g}\) is the gradient vector.

The step size can be determined thought a line search algorithm using the Armijo rule (Kelley, 1999). In this case,

\[\lambda = \beta^m\]

where \(1 > \beta > 0\) and \(m \ge 0\) is an integer that controls the step size. The line search finds the smallest \(m\) that satisfies the Armijo rule

\[\phi(\bar{p} + \Delta\bar{p}) - \Gamma(\bar{p}) < \alpha\beta^m ||\bar{g}(\bar{p})||^2\]

where \(\phi(\bar{p})\) is the objective function evaluated at \(\bar{p}\) and \(\alpha = 10^{-4}\).

Parameters:

• gradientfunction

  A function that returns the gradient vector of the objective function when given a parameter vector.

• valuefunction

  A function that returns the value of the objective function evaluated at a given parameter vector.

• initial1d-array

  The initial estimate for the gradient descent.

• maxitint

  The maximum number of iterations allowed.

• linesearchTrue or False

  Whether or not to perform the line search to determine an optimal step size.

• maxstepsint

  The maximum number of times to try to take a step before giving up.

• betafloat

  The base factor used to determine the step size in line search algorithm. Must be 1 > beta > 0.

• tolfloat

  The convergence criterion. The lower it is, the more steps are permitted.

Yields:

• i, estimate, stats:

  • iint

    The current iteration number

  • estimate1d-array

    The current estimated parameter vector

  • statsdict

    Statistics about the optimization so far. Keys:

    • methodstf

      The name of the optimization algorithm

    • iterationsint

      The total number of iterations so far

    • objectivelist

      Value of the objective function per iteration. First value corresponds to the inital estimate

    • step_attemptslist

      Number of attempts at taking a step per iteration. First number is zero, reflecting the initial estimate. Will be empty if linesearch==False.

References:

Kelley, C. T., 1999, Iterative methods for optimization: Raleigh: SIAM.

geoist.inversion.pfmodel module

3D density structure inversion with gravity anomaly.

class geoist.inversion.pfmodel.AbicLSQOperator(toep, depth_constraint=None, smooth_components={}, refer_constraint=None, weights=None)

基类：geoist.inversion.toeplitz.LSQOperator

An operator doing matrix vector multiplication. The matrix is: $lpha_g G^TG + sum lpha_i W^TB_i^TB_iW$. Where $lpha$'s are weights, $G$ is kernel matrix, $W$ is depth constraint, $B_i$'s are other constrains.

matvec(v)

If do_precond is True, apply T^T*M^{-T}*M^{-1}*T + coef*diag**2 to vector v, else apply T^T*T + coef*diag**2 to vector v

class geoist.inversion.pfmodel.AbicLSQOperator2(op_in)

基类：object

matvec(v)

class geoist.inversion.pfmodel.InvModel(nzyx=[4, 4, 4], smooth_components=None, depth_constraint=None, model_density=None, refer_density=None, weights=None, source_volume=None, smooth_on='m', subtract_mean=True, data_dir='/data/gravity_inversion')

基类：object

Potential model define for inversion. eg. m1 = InvModel(nzyx, ...)

abic_optimize()

bound_optimize(x0=None)

calc_abic()

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_abic_quiet()

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_log_obs_det()

calc_log_obs_det_quiet()

calc_log_prior_total_det()

calc_log_prior_total_det_quiet()

calc_min_u(solved=False, x=None)

calc_res()

calc_u(solved=False, x=None)

calc_u_quiet(solved=False, x=None)

property depth_constraint

do_linear_solve()

do_linear_solve_quiet()

forward(model_density=None)

gen_kernel(process=1)

gen_mesh(height=- 1)

gen_model_name()

hessp_u(x, v)

jac_u(x=None)

lasso_hessp(x, v)

lasso_jac(x)

lasso_optimize(x0=None)

lasso_target(x)

property model_density

property nx

property ny

property nz

para_grad(x)

print_summary()

property refer_density

property smooth_components

property source_volume

u_bound()

walsh_transform(keys=None)

property weights

class geoist.inversion.pfmodel.SmoothOperator

基类：object

SmoothOperator for ABIC.

derivation(v, component='dx', array_order='zyx')

rderivation(v, component='dx', array_order='zyx')

geoist.inversion.pfmodel.check_gpu()

Check GPU setup.

geoist.inversion.pfmodel.free_gpu()

free up gpu memory consumption

geoist.inversion.pfmodel.timeit(f)

geoist.inversion.pfmodel_ts module

class geoist.inversion.pfmodel_ts.AbicLSQOperatorMS(kernel_matrices, ns, nzyx, smooth_components={}, weights=None)

基类：object

An operator doing matrix vector multiplication. The matrix is:

$lpha_g G^TG + sum lpha_i B_i^TB_i$. Where $lpha$'s are weights, $G$ is kernel matrix, $B_i$'s are smooth matrix.

参数

• kernel_matrix (list of ndarray) -- kernel matrix of each survey

• ns (int) -- number of surveys

• nzyx (tuple of int) -- number of cells along z-axis, y-axis and x-axis

• smmoth_components (list of str) -- which components should be smoothed, acceptable string could be 'dx','dxx','dy','dyy','dxy','dt',...

• weights (dict) -- name and values of weights.

matvec(v)

class geoist.inversion.pfmodel_ts.InvModelTS(nyx=20, 20, smooth_components=None, weights=None, optimize_weights=None, source_volume=None, margin=0, 0, 0, 0, cell_type='prism', data_dir='./data')

基类：object

Multiple survey gravity inversion use Abic :param nyx: number of cells along y-axis and x-axis :type nyx: tuple of int :param smooth_components: which components should be smoothed, acceptable

string could be 'dx','dxx','dy','dyy','dxy','dt',...

参数

• weights (dict) -- name and values of weights.

• optimize_weights (list of str) -- specify which weight should be optimized by Abic.

• source_volume (tuple of float) -- dimension of the underlying gravity source.

• margin (tuple of float) -- margin size around the source volume.

• cell_type (str) -- either 'prism' or 'Tesseroid'

• data_dir (str) -- a folder where observation data resides

ns

number of surveys

Type

int

nx,ny,nz

number of cells along x-axis, y-axis and z-axis

Type

int

source_volume

dimension of the underlying gravity source

Type

tuple of float

margin

margin size around the source volume.

Type

tuple of float

cell_type

either 'prism' or 'Tesseroid'

Type

str

smooth_components

which components should be smoothed

Type

list of str

weights

name and values of weights.

Type

dict

optimize_weights

specify which weight should be optimized by Abic.

Type

list of str

data_dir

a folder where observation data resides

Type

str

abic_val

abic value

Type

float

log_total_det_val

log determinate value of the total matrix

Type

float

log_obs_det_val

log determinate value of the observation matrix

Type

float

min_u_val

minimum of U

Type

float

min_density

minimum density

Type

float

max_density

maximum density

Type

float

kernel_matrices

kernels of each survey.

Type

list of ndarray

rhs

right hand side when solve min_u

Type

ndarray

solution

model densityies solve min_u

Type

ndarray

abic_log

save the weights and abic value during calculation

Type

dict

orig_data

original data loaded from files.

Type

DataFrame

mesh

mesh of the underlying density model.

Type

Mesh

R

radius of the Earth (m)

Type

float

abic_optimize()

calc_abic()

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_abic_quiet(precision=1e-06)

-log_prior_det_value+log_total_det-log_obs_det+min_u

calc_log_obs_det()

calc_log_obs_det_quiet()

calc_log_prior_total_det(precision=1e-06)

calc_log_prior_total_det_quiet(precision=1e-06)

calc_min_u(solved=False, x=None)

calc_res()

calc_u(solved=False, x=None)

calc_u_quiet(solved=False, x=None)

deg2xy()

do_linear_solve()

do_linear_solve_quiet()

forward(model_density=None)

gen_kernel()

gen_mesh(source_volume=None, margin=None)

gen_obs_grid(height=- 1)

generate obs grid

load_data(pattern='*.txt', names=['lon', 'lat', 'g'], **kwargs)

load all data inside data dir, result saved in self.orig_data. :param pattern: filename pattern. :type pattern: str :param names: name of the observation must be 'g', coordinates

must be 'lon','lat' or 'x','y'

参数

kwargs -- same as pd.read_csv

property nx

property ny

para_grad(x)

plot_density(density=None, surveys=None, fname=None)

plot_field(field=None, surveys=None, fname=None, plot_station=True)

print_summary()

property smooth_components

u_bound()

property weights

class geoist.inversion.pfmodel_ts.SmoothOperator

基类：object

derivation(v, component='dx', array_order='tzyx')

rderivation(v, component='dx', array_order='tzyx')

geoist.inversion.pfmodel_ts.timeit(f)

geoist.inversion.regularization module

Ready made classes for regularization.

Each class represents a regularizing function. They can be used by adding them to a Misfit derivative (all inversions in geoist are derived from Misfit). The regularization parameter is set by multiplying the regularization instance by a scalar, e.g., solver = misfit + 0.1*regularization.

See geoist.gravmag.eqlayer.EQLGravity for an example.

List of classes

• Damping: Damping regularization (or 0th order Tikhonov regularization)

• Smoothness: Generic smoothness regularization (or 1st order Tikhonov regularization). Requires a finite difference matrix to specify the parameter derivatives to minimize.

• Smoothness1D: Smoothness for 1D problems. Automatically builds a finite difference matrix based on the number of parameters

• Smoothness2D: Smoothness for 2D grid based problems. Automatically builds a finite difference matrix of derivatives in the two spacial dimensions based on the shape of the parameter grid

• TotalVariation: Generic total variation regularization (enforces sharpness of the solution). Requires a finite difference matrix to specify the parameter derivatives.

• TotalVariation1D: Total variation for 1D problems. Similar to Smoothness1D

• TotalVariation2D: Total variation for 2D grid based problems. Similar to Smoothness2D

---

class geoist.inversion.regularization.Damping(nparams)

基类：geoist.inversion.regularization.Regularization

Damping (0th order Tikhonov) regularization.

Imposes the minimum norm of the parameter vector.

The regularizing function if of the form

\[\theta^{NM}(\bar{p}) = \bar{p}^T\bar{p}\]

Its gradient and Hessian matrices are, respectively,

\[\bar{\nabla}\theta^{NM}(\bar{p}) = 2\bar{\bar{I}}\bar{p}\]

\[\bar{\bar{\nabla}}\theta^{NM}(\bar{p}) = 2\bar{\bar{I}}\]

Parameters:

• nparamsint

  The number of parameter

Examples:

>>> import numpy
>>> damp = Damping(3)
>>> p = numpy.array([0, 0, 0])
>>> damp.value(p)
0.0
>>> damp.hessian(p).todense()
matrix([[ 2.,  0.,  0.],
        [ 0.,  2.,  0.],
        [ 0.,  0.,  2.]])
>>> damp.gradient(p)
array([ 0.,  0.,  0.])
>>> p = numpy.array([1, 0, 0])
>>> damp.value(p)
1.0
>>> damp.hessian(p).todense()
matrix([[ 2.,  0.,  0.],
        [ 0.,  2.,  0.],
        [ 0.,  0.,  2.]])
>>> damp.gradient(p)
array([ 2.,  0.,  0.])

The Hessian matrix is cached so that it is only generated on the first call to damp.hessian (unlike the gradient, which is calculated every time).

>>> damp.hessian(p) is damp.hessian(p)
True
>>> damp.gradient(p) is damp.gradient(p)
False

gradient(p)

Calculate the gradient vector.

Parameters:

• p1d-array or None

  The parameter vector. If None, will return 0.

Returns:

• gradient1d-array

  The gradient

hessian(p)

Calculate the Hessian matrix.

Parameters:

• p1d-array

  The parameter vector

Returns:

• hessian2d-array

  The Hessian

value(p)

Calculate the value of this function.

Parameters:

• p1d-array

  The parameter vector

Returns:

• valuefloat

  The value of this function evaluated at p

class geoist.inversion.regularization.Regularization(nparams, islinear)

基类：geoist.inversion.base.OperatorMixin

Generic regularizing function.

When subclassing this, you must implemenent the value method.

copy(deep=False)

Make a copy of me together with all the cached methods.

class geoist.inversion.regularization.Smoothness(fdmat)

基类：geoist.inversion.regularization.Regularization

Smoothness (1st order Tikhonov) regularization.

Imposes that adjacent parameters have values close to each other.

The regularizing function if of the form

\[\theta^{SV}(\bar{p}) = \bar{p}^T \bar{\bar{R}}^T \bar{\bar{R}}\bar{p}\]

Its gradient and Hessian matrices are, respectively,

\[\bar{\nabla}\theta^{SV}(\bar{p}) = 2\bar{\bar{R}}^T \bar{\bar{R}}\bar{p}\]

\[\bar{\bar{\nabla}}\theta^{SV}(\bar{p}) = 2\bar{\bar{R}}^T \bar{\bar{R}}\]

in which matrix \(\bar{\bar{R}}\) is a finite difference matrix. It represents the differences between one parameter and another and is what indicates what adjacent means.

Parameters:

• fdmat2d-array or sparse matrix

  The finite difference matrix

Examples:

>>> import numpy as np
>>> fd = np.array([[1, -1, 0],
...                [0, 1, -1]])
>>> smooth = Smoothness(fd)
>>> p = np.array([0, 0, 0])
>>> smooth.value(p)
0.0
>>> smooth.gradient(p)
array([0, 0, 0])
>>> smooth.hessian(p)
array([[ 2, -2,  0],
       [-2,  4, -2],
       [ 0, -2,  2]])
>>> p = np.array([1, 0, 1])
>>> smooth.value(p)
2.0
>>> smooth.gradient(p)
array([ 2, -4,  2])
>>> smooth.hessian(p)
array([[ 2, -2,  0],
       [-2,  4, -2],
       [ 0, -2,  2]])

The Hessian matrix is cached so that it is only generated on the first call to hessian (unlike the gradient, which is calculated every time).

>>> smooth.hessian(p) is smooth.hessian(p)
True
>>> smooth.gradient(p) is smooth.gradient(p)
False

gradient(p)

Calculate the gradient vector.

Parameters:

• p1d-array or None

  The parameter vector. If None, will return 0.

Returns:

• gradient1d-array

  The gradient

hessian(p)

Calculate the Hessian matrix.

Parameters:

• p1d-array

  The parameter vector

Returns:

• hessian2d-array

  The Hessian

value(p)

Calculate the value of this function.

Parameters:

• p1d-array

  The parameter vector

Returns:

• valuefloat

  The value of this function evaluated at p

class geoist.inversion.regularization.Smoothness1D(npoints)

基类：geoist.inversion.regularization.Smoothness

Smoothness regularization for 1D problems.

Extends the generic Smoothness class by automatically building the finite difference matrix.

Parameters:

• npointsint

  The number of parameters

Examples:

>>> import numpy as np
>>> s = Smoothness1D(3)
>>> p = np.array([0, 0, 0])
>>> s.value(p)
0.0
>>> s.gradient(p)
array([0, 0, 0])
>>> s.hessian(p).todense()
matrix([[ 2, -2,  0],
        [-2,  4, -2],
        [ 0, -2,  2]])
>>> p = np.array([1, 0, 1])
>>> s.value(p)
2.0
>>> s.gradient(p)
array([ 2, -4,  2])
>>> s.hessian(p).todense()
matrix([[ 2, -2,  0],
        [-2,  4, -2],
        [ 0, -2,  2]])

class geoist.inversion.regularization.Smoothness2D(shape)

基类：geoist.inversion.regularization.Smoothness

Smoothness regularization for 2D problems.

Extends the generic Smoothness class by automatically building the finite difference matrix.

Parameters:

• shapetuple = (ny, nx)

  The shape of the parameter grid. Number of parameters in the y and x (or z and x, time and offset, etc) dimensions.

Examples:

>>> import numpy as np
>>> s = Smoothness2D((2, 2))
>>> p = np.array([[0, 0],
...               [0, 0]]).ravel()
>>> s.value(p)
0.0
>>> s.gradient(p)
array([0, 0, 0, 0])
>>> s.hessian(p).todense()
matrix([[ 4, -2, -2,  0],
        [-2,  4,  0, -2],
        [-2,  0,  4, -2],
        [ 0, -2, -2,  4]])
>>> p = np.array([[1, 0],
...               [2, 3]]).ravel()
>>> s.value(p)
12.0
>>> s.gradient(p)
array([ 0, -8,  0,  8])
>>> s.hessian(p).todense()
matrix([[ 4, -2, -2,  0],
        [-2,  4,  0, -2],
        [-2,  0,  4, -2],
        [ 0, -2, -2,  4]])

class geoist.inversion.regularization.TotalVariation(beta, fdmat)

基类：geoist.inversion.regularization.Regularization

Total variation regularization.

Imposes that adjacent parameters have a few sharp transitions.

The regularizing function if of the form

\[\theta^{VT}(\bar{p}) = \sum\limits_{k=1}^L |v_k|\]

where vector \(\bar{v} = \bar{\bar{R}}\bar{p}\). See Smoothness for the definition of the \(\bar{\bar{R}}\) matrix.

This functions is not differentiable at the null vector, so the following form is used to calculate the gradient and Hessian

\[\theta^{VT}(\bar{p}) \approx \theta^{VT}_\beta(\bar{p}) = \sum\limits_{k=1}^L \sqrt{v_k^2 + \beta}\]

Its gradient and Hessian matrices are, respectively,

\[\bar{\nabla}\theta^{VT}_\beta(\bar{p}) = \bar{\bar{R}}^T \bar{q}(\bar{p})\]

\[\bar{\bar{\nabla}}\theta^{VT}_\beta(\bar{p}) = \bar{\bar{R}}^T \bar{\bar{Q}}(\bar{p})\bar{\bar{R}}\]

and

\[\begin{split}\bar{q}(\bar{p}) = \begin{bmatrix} \dfrac{v_1}{\sqrt{v_1^2 + \beta}} \\ \dfrac{v_2}{\sqrt{v_2^2 + \beta}} \\ \vdots \\ \dfrac{v_L}{\sqrt{v_L^2 + \beta}} \end{bmatrix}\end{split}\]

\[\begin{split}\bar{\bar{Q}}(\bar{p}) = \begin{bmatrix} \dfrac{\beta}{(v_1^2 + \beta)^{\frac{3}{2}}} & 0 & \ldots & 0 \\ 0 & \dfrac{\beta}{(v_2^2 + \beta)^{\frac{3}{2}}} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \dfrac{\beta}{(v_L^2 + \beta)^{\frac{3}{2}}} \end{bmatrix}\end{split}\]

Parameters:

• betafloat

  The beta parameter for the differentiable approximation. The larger it is, the closer total variation is to Smoothness. Should be a small, positive value.

• fdmat2d-array or sparse matrix

  The finite difference matrix

gradient(p)

Calculate the gradient vector.

Parameters:

• p1d-array

  The parameter vector.

Returns:

• gradient1d-array

  The gradient

hessian(p)

Calculate the Hessian matrix.

Parameters:

• p1d-array

  The parameter vector

Returns:

• hessian2d-array

  The Hessian

value(p)

Calculate the value of this function.

Parameters:

• p1d-array

  The parameter vector

Returns:

• valuefloat

  The value of this function evaluated at p

class geoist.inversion.regularization.TotalVariation1D(beta, npoints)

基类：geoist.inversion.regularization.TotalVariation

Total variation regularization for 1D problems.

Extends the generic TotalVariation class by automatically building the finite difference matrix.

Parameters:

• betafloat

  The beta parameter for the differentiable approximation. The larger it is, the closer total variation is to Smoothness. Should be a small, positive value.

• npointsint

  The number of parameters

class geoist.inversion.regularization.TotalVariation2D(beta, shape)

基类：geoist.inversion.regularization.TotalVariation

Total variation regularization for 2D problems.

Extends the generic TotalVariation class by automatically building the finite difference matrix.

Parameters:

• betafloat

  The beta parameter for the differentiable approximation. The larger it is, the closer total variation is to Smoothness. Should be a small, positive value.

• shapetuple = (ny, nx)

  The shape of the parameter grid. Number of parameters in the y and x (or z and x, time and offset, etc) dimensions.

geoist.inversion.regularization.fd1d(size)

Produce a 1D finite difference matrix.

Parameters:

• sizeint

  The number of points

Returns:

• fdsparse CSR matrix

  The finite difference matrix

Examples:

>>> fd1d(2).todense()
matrix([[ 1, -1]])
>>> fd1d(3).todense()
matrix([[ 1, -1,  0],
        [ 0,  1, -1]])
>>> fd1d(4).todense()
matrix([[ 1, -1,  0,  0],
        [ 0,  1, -1,  0],
        [ 0,  0,  1, -1]])

geoist.inversion.regularization.fd2d(shape)

Produce a 2D finite difference matrix.

Parameters:

• shapetuple = (ny, nx)

  The shape of the parameter grid. Number of parameters in the y and x (or z and x, time and offset, etc) dimensions.

Returns:

• fdsparse CSR matrix

  The finite difference matrix

Examples:

>>> fd2d((2, 2)).todense()
matrix([[ 1, -1,  0,  0],
        [ 0,  0,  1, -1],
        [ 1,  0, -1,  0],
        [ 0,  1,  0, -1]])
>>> fd2d((2, 3)).todense()
matrix([[ 1, -1,  0,  0,  0,  0],
        [ 0,  1, -1,  0,  0,  0],
        [ 0,  0,  0,  1, -1,  0],
        [ 0,  0,  0,  0,  1, -1],
        [ 1,  0,  0, -1,  0,  0],
        [ 0,  1,  0,  0, -1,  0],
        [ 0,  0,  1,  0,  0, -1]])

geoist.inversion.toeplitz module

class geoist.inversion.toeplitz.GToepOperator(toep, *args, **kwargs)

基类：object

This class construct an linear operator of multilevel symmetric Toeplitz matrix.

cleanup()

get_m_eigs(toep)

M is designed to be an matrix approximate the underlying Toeplitz matrix. So that M could be used as preconditioner.

matvec(v)

Apply the underlying Toeplitz matrix T to a multivector v.

matvec_prec(v)

Apply M^{-1} to a vector v.

matvec_prec_sym(v)

Apply M^{-T}M^{-1} to a vector v.

rmatvec(v)

Apply the transpose of the underlying Toeplitz matrix (i.e. T^T) to a vector v.

class geoist.inversion.toeplitz.LSQOperator(toep, diag=None, coef=1.0, do_precond=False, gtoep=None, *args, **kwargs)

基类：object

This class used to define a Least Square problem by a multilevel symmetric Toepliz matrix.

cleanup()

matvec(v)

If do_precond is True, apply T^T*M^{-T}*M^{-1}*T + coef*diag**2 to vector v, else apply T^T*T + coef*diag**2 to vector v

geoist.inversion.toeplitz.block_circ(a)

generate full representation of 2-level circulant matrix

参数

a (ndarray) -- 1-st column of the circulant matrix in proper shape.

返回

Full filled circulant matrix

geoist.inversion.toeplitz.block_toep2_sym(a)

generate full representation of 2-level symmetric toeplitz matrix

参数

a (ndarray) -- 1-st column of the symmetrec toeplitz matrix in proper shape.

返回

Full filled toeplitz matrix.

geoist.inversion.toeplitz.cg(A, b, x=None, tol=1e-05, max_iter=None)

geoist.inversion.toeplitz.check_gpu()

geoist.inversion.toeplitz.circ_eigs(circ, dtype=<class 'numpy.complex64'>)

calculate eigenvalues of multilevel circulant matrix A

参数

circ (ndarray) -- representation of the multilevel circulant matrix A, i.e. the first column of A in proper shape.

返回

eigenvalues of A.

geoist.inversion.toeplitz.circ_mul_v(circ, v, eigs=None)

multiply a circulant matrix A by multi vector v.

参数

• circ (ndarray) -- representation of the multilevel circulant matrix A, i.e. the first column of A in proper shape.

• v (ndarray) -- vector to be multiplied. Should be reshaped to the same shape as circ. Should be the same reshape order (column first/row first) as circ.

返回

result of multiplication.

geoist.inversion.toeplitz.embed_toep2circ(toep, v=None)

embed toeplitz matrix to circulant matrix.

参数

• toep (ndarray) -- representation of multilevel toeplitz matrix, i.e. the first column of A in proper shape.

• v (ndarray) -- embed a vector together with toep.

返回

representation of embedded multilevel circulant matrix and embedded vector.

geoist.inversion.toeplitz.toeplitz_mul_v(toep, v)

multiply a toeplitz matrix A by vector v.

参数

• toep (ndarray) -- representation fo multilevel toeplitz matrix, i.e. the first column of A in proper shape.

• v (ndarray) -- vector to be multiplied. Should be reshaped to the same shape as toep. Should be the same reshape order (column first/row first) as circ.

返回

result of multiplication.

geoist.inversion.walsh module

geoist.inversion.walsh.fast_walsh_transform(a, normalized=True)

inplace FFT-like walsh transform. Let M be the walsh matrix, this routine will return M*a. This routine is designed to transform a vector. Although it still works on matrix, it misses the multiplier M^{-1} on the right side. If you want transform matrix A, first run fast_walsh_transform(A), then fast_walsh_transform(A.T). This routine use natural ordering walsh matrix.

参数

a (ndarray) -- ndarray to be transformed

返回

M*a. where M is the walsh matrix

geoist.inversion.walsh.linear_recover(M, IM, A=None, x=None, b=None)

Recover equation Ax=b from A'x'=b' where A'=M*A*M^{-1}

x'=M*x b'=M*b.

参数

• M (ndarray) -- trasform matrix.

• IM (ndarray) -- inverse of trasform matrix.

• A (ndarray) -- coefficient matrix A'.

• x (ndarray) -- unkowns x'.

• b (ndarray) -- righthand side b'.

返回

recovered A,x,b.

geoist.inversion.walsh.linear_transform(M, IM, A=None, x=None, b=None)

Transform equation Ax=b to A'x'=b' where A'=M*A*M^{-1}

x'=M*x b'=M*b.

参数

• M (ndarray) -- trasform matrix.

• IM (ndarray) -- inverse of trasform matrix.

• A (ndarray) -- coefficient matrix.

• x (ndarray) -- unkowns.

• b (ndarray) -- righthand side.

返回

transformed A',x',b'.

geoist.inversion.walsh.natural_walsh_matrix(n, normalized=False)

natural ordering walsh matrix

参数

• n (int) -- degree of walsh matrix.

• normalized (bool) -- whether or not normalize the generated matrix.

返回

walsh matrix in natural ordering.

geoist.inversion.walsh.natural_walsh_matrix_old(n, normalized=False)

natural ordering walsh matrix

参数

• n (int) -- degree of walsh matrix.

• normalized (bool) -- whether or not normalize the generated matrix.

返回

walsh matrix in natural ordering.

geoist.inversion.walsh.walsh_matrix(n, normalized=False, ordering='sequence', nxyz=None, None, None)

generate walsh matrix with given ordering.

参数

• n (int) -- degree of walsh matrix.

• normalized (bool) -- whether or not nomalize walsh matrix.

• ordering (string) -- ordering of walsh matrix. take values from 'natural', 'dyadic','sequence'

• nxyz (tuple of ints) -- (nx,ny,nz) specifies source points configuration.

返回

walsh matrix with given ordering.

geoist.inversion.walsh.walsh_matrix_reorder(m, origin_order='natural', target_order='sequence')

reorder walsh matrix

参数

• m (ndarray) -- walsh matrix

• from (string) -- orginal ordering. take values from 'natural','dyadic', 'sequence'

• to (string) -- target ordering. take values from 'natural','dyadic', 'sequence'

返回

reordered walsh matrix.

geoist.inversion.walsh.walsh_order(n)

generate 'natural','dyadic','sequence' ordering of walsh matrix.

参数

n (int) -- degree of walsh matrix.

geoist.inversion.walsh.walsh_recover(n, M=None, A=None, x=None, b=None, normalized=True, ordering='sequence2', nxyz=None, None, None)

Recover linear equation Ax=b from A'x'=b' where A'=M*A*M^{-1}

x'=M*x b'=M*b. M is the walsh matrix.

参数

• n (int) -- degree of walsh matrix.

• M (ndarray) -- trasform matrix.

• A (ndarray) -- coefficient matrix.

• x (ndarray) -- unkowns.

• b (ndarray) -- righthand side.

• ordering (string) -- ordering of walsh matrix. could be one of 'natural','dyadic','sequence','sequence2'.

• nxyz (tuple of ints) -- (nx,ny,nz) specifies source points configuration.

返回

recovered A,x,b.

geoist.inversion.walsh.walsh_reorder2(m, nxyz)

reorder walsh matrix from sequence order to sequence2 order. In sequence2 ordering, kernel.T*kernel will be a block diagonal matrix.

参数

• m (ndarray) -- walsh matrix in sequence ordering.

• nxyz (tuple of ints) -- (nx,ny,nz) specifies source points configuration.

返回

walsh matrix in sequence2.

geoist.inversion.walsh.walsh_transform(A=None, x=None, b=None, normalized=True, ordering='sequence2', nxyz=None, None, None)

Transform linear equation Ax=b to A'x'=b' where A'=M*A*M^{-1}

x'=M*x b'=M*b. M is the walsh matrix.

参数

• n (int) -- degree of walsh matrix.

• M (ndarray) -- trasform matrix.

• A (ndarray) -- coefficient matrix.

• x (ndarray) -- unkowns.

• b (ndarray) -- righthand side.

• ordering (string) -- ordering of walsh matrix. could be one of 'natural','dyadic','sequence','sequence2'.

• nxyz (tuple of ints) -- (nx,ny,nz) specifies source points configuration.

返回

transformed A',x',b'.

geoist.inversion.walshwavelet module

Created on Fri Apr 12 11:06:23 2019

@author: chens

geoist.inversion.walshwavelet.grav_matrix(number, normalized=False)

generate walsh matrix

Module contents

Everything you need to solve inverse problems!

This package provides the basic building blocks to implement inverse problem solvers. The main class for this is Misfit. It represents a data-misfit function and contains various tools to fit a model to some data. All you have to do is implement methods to calculate the predicted (or modeled) data and (optionally) the Jacobian (or sensitivity) matrix. With only that, you have access to a range of optimization methods, regularization, joint inversion, etc.

Modules

• misfit: Defines the data-misfit classes. Used to implement new inverse problems. The main class is Misfit. It offers a template for you to create standard least-squares inversion methods.

• regularization: Classes for common regularizing functions and base classes for building new ones.

• hyper_param: Classes hyper parameter optimization (estimating the regularization parameter), like L-curve analysis and (in the future) cross-validation.

• optimization: Functions for several optimization methods (used internally by Misfit). In most cases you won't need to touch this.

• base: Base classes used internally to define common operations and method. In most cases you won't need to touch this.

You can import the Misfit, regularization, and hyper parameter optimization classes directly from the geoist.inversion namespace:

>>> from geoist.inversion import Misfit, LCurve, Damping, Smoothness

The Misfit class is used internally in geoist to implement all of our inversion algorithms. Take a look at the modules below for some examples:

• geoist.seismic.srtomo

• geoist.gravmag.basin2d

• geoist.gravmag.eqlayer

• geoist.gravmag.euler

• geoist.gravmag.magdir

• geoist.seismic.profile

实际案例

The Misfit class works by subclassing it. Doing this gives you access to all optimization functions and possible combinations of regularization. When subclassing Misfit, you'll need to implement the predicted method that calculates a predicted data vector based on an input parameter vector. This is virtually all that is problem-specific in an inverse problem. If you want to use gradient-based optimization (or linear problems) you'll need to implement the jacobian method as well that calculates the Jacobian (or sensitivity) matrix.

Linear Regression

Here is an example of how to implement a simple linear regression using the Misfit class.

What we want to do is fit \(f(a, b, x) = y = ax + b\) to find a and b. Putting a and b in a parameter vector p = [a, b] we get:

\[\mathbf{d} = \mathbf{A} \mathbf{p}\]

where \(\mathbf{d}\) is the data vector containing all the values of y and \(\mathbf{A}\) is the Jacobian matrix with the values of x in the first column and 1 in the second column.

All we have to do to implement a solver for this problem is write the predicted (to calculate y from values of a and b) and jacobian (to calculate the Jacobian matrix):

First, I'll load numpy and the Misfit class.

>>> import numpy as np
>>> from geoist.inversion import Misfit

We'll also need some compatibility code to support Python 2 and 3 at the same time.

>>> from __future__ import division
>>> from future.builtins import super

Now, I'll make my regression class that inherits from Misfit. All of the least-squares solving code and much more we get for free just by using Misfit as template for our regression class. Note Misfit wants a 1D data vector, no matter what your data is (line, grid, volume).

>>> class Regression(Misfit):
...     "Perform a linear regression"
...     def __init__(self, x, y):
...         # Call the initialization of Misfit
...         super().__init__(data=y, nparams=2, islinear=True)
...         self.x = x  # Store this to use in the other methods
...     def predicted(self, p):
...         "Calculate the predicted data for a given parameter vector"
...         a, b = p
...         return a*self.x + b
...     def jacobian(self, p):
...         "Calculate the Jacobian (ignores p because the problem is linear)"
...         jac = np.empty((self.ndata, self.nparams))
...         jac[:, 0] = self.x
...         jac[:, 1] = 1
...         return jac

To test our regression, let's generate some data based on known values of a and b.

>>> x = np.linspace(0, 5, 6)
>>> x
array([ 0.,  1.,  2.,  3.,  4.,  5.])
>>> y = 2*x + 5
>>> y
array([  5.,   7.,   9.,  11.,  13.,  15.])

We must create a solver object to perform our regression. Fitting the data (running the optimization) is done by calling the fit method. fit returns the regression class itself we can chain calls to it like so:

>>> solver = Regression(x, y).fit()

The estimated parameter vector is stored in the p_ attribute. Misfit also provides a estimate_ attribute that can be a custom (user defined) formatted version of p_. It's better to use estimate_ if you're not interested in the parameter vector as it is. Since we didn't implement this custom formatting, both should give the same value.

>>> solver.p_
array([ 2.,  5.])
>>> solver.estimate_
array([ 2.,  5.])

The methods predicted and residuals will use the cached values based on p_ if the parameter vector is omitted as an argument.

>>> solver.predicted()
array([  5.,   7.,   9.,  11.,  13.,  15.])
>>> np.all(np.abs(solver.residuals()) < 10**10)
True

Caching

The Misfit class caches some of the matrices that it computes, like the Jacobian matrix. This is needed because separate methods call jacobian with the same input p, so recomputing the matrix would be a waste.

For linear problems, the Jacobian matrix is cached permanently. It is only calculated once, no matter what p is passed to it.

>>> A = solver.jacobian(solver.p_)
>>> A
array([[ 0.,  1.],
       [ 1.,  1.],
       [ 2.,  1.],
       [ 3.,  1.],
       [ 4.,  1.],
       [ 5.,  1.]])
>>> B = solver.jacobian(np.array([20, 30]))
>>> B
array([[ 0.,  1.],
       [ 1.,  1.],
       [ 2.,  1.],
       [ 3.,  1.],
       [ 4.,  1.],
       [ 5.,  1.]])
>>> A is B
True

For non-linear problems, the Jacobian is also cached but it will be recalculated if passed a different value of p (see the non-linear example below).

The Hessian matrix (\(2\mathbf{A}^T\mathbf{A}\)) is also cached permanently for linear problems.

>>> H = solver.hessian(solver.p_)
>>> H
array([[ 110.,   30.],
       [  30.,   12.]])
>>> H2 = solver.hessian(np.array([20, 30]))
>>> H2
array([[ 110.,   30.],
       [  30.,   12.]])
>>> H is H2
True

Non-linear optimization

You can configure the solver using the config method. This allows you to choose the optimization method you want to use and it's corresponding parameters. The config method also returns the solver class itself so it can be chained with fit:

>>> # Configure solver to use the Levemberg-Marquardt method
>>> solver.config('levmarq', initial=[1, 1]).fit().estimate_
array([ 2.,  5.])
>>> np.all(np.abs(solver.residuals()) < 10**10)
True
>>> # or the Steepest Descent method
>>> solver.config('steepest', initial=[1, 1]).fit().estimate_
array([ 2.,  5.])
>>> # or the Gauss-Newton method
>>> solver.config('newton', initial=[1, 1], maxit=5).fit().estimate_
array([ 2.,  5.])

The Misfit class keeps track of the optimization process and makes that data available through the stats_ attribute (a dictionary).

>>> list(sorted(solver.stats_))
['iterations', 'method', 'objective']
>>> solver.stats_['method']
"Newton's method"
>>> solver.stats_['iterations']
5

The 'objective' key holds a list of the objective function value per iteration of the optimization process.

Re-weighted least squares

Misfit allows you to set weights to the data in the form of a weight matrix or vector (the vector is assumed to be the diagonal of the weight matrix). We can use this to perform a re-weighted least-squares fit to remove outliers from our data.

>>> y_outlier = y.copy()
>>> y_outlier[3] += 20
>>> y_outlier
array([  5.,   7.,   9.,  31.,  13.,  15.])

First, we run the regression without any weights.

>>> solver = Regression(x, y_outlier).fit()
>>> print(np.array_repr(solver.estimate_, precision=3))
array([ 2.571,  6.905])

Now, we can use the inverse of the residuals to set the weights for our data. We repeat this for a few iterations and should have our robust estimate by the end of it.

>>> for i in range(20):
...     r = np.abs(solver.residuals())
...     # Avoid small residuals because of zero-division errors
...     r[r < 1e-10] = 1
...     _ = solver.set_weights(1/r).fit()
>>> solver.estimate_
array([ 2.,  5.])

Non-linear problems

In this example, I want to fit a Gaussian to my data:

\[f(x) = a\exp(-b(x + c)^{2})\]

Function f is non-linear with respect to inversion parameters a, b, c. Thus, we need to configure the solver and choose an optimization method before we can call fit().

First, lets create our solver class based on Misfit and implement the predicted and jacobian methods.

>>> class GaussianFit(Misfit):
...     def __init__(self, x, y):
...         super().__init__(
...             data=y, nparams=3, islinear=False)
...         self.x = x
...     def predicted(self, p):
...         a, b, c = p
...         return a*np.exp(-b*(self.x + c)**2)
...     def jacobian(self, p):
...         a, b, c = p
...         jac = np.empty((self.ndata, self.nparams))
...         var = self.x + c
...         exponential = np.exp(-b*var**2)
...         jac[:, 0] = exponential
...         jac[:, 1] = -a*exponential*(var**2)
...         jac[:, 2] = -a*exponential*2*b*var
...         return jac

Let's create some data to test this.

>>> x = np.linspace(0, 10, 1000)
>>> a, b, c = 100, 0.1, -2
>>> y = a*np.exp(-b*(x + c)**2)

For non-linear problems, we have to configure the optimization method. Lets use Levemberg-Marquardt because it generally offers good convergence.

>>> solver = GaussianFit(x, y).config('levmarq', initial=[1, 1, 1]).fit()
>>> # Print the estimated coefficients
>>> print(', '.join(['{:.1f}'.format(i) for i in solver.estimate_]))
100.0, 0.1, -2.0
>>> np.all(np.abs(solver.residuals()) < 10**-10)
True

We can use other optimization methods without having to re-implement our solution. For example, let's see how well the Ant Colony Optimization for Continuous Domains (ACO-R) does for this problem:

>>> # bounds are the min, max values of the search domain for each parameter
>>> _ = solver.config('acor', bounds=[50, 500, 0, 1, -20, 0], seed=0).fit()
>>> print(', '.join(['{:.1f}'.format(i) for i in solver.estimate_]))
100.0, 0.1, -2.0

For non-linear problems, the Jacobian and Hessian are cached but not permanently. Calling jacobian twice in a row with the same parameter vector will not trigger a computation and will return the cached value instead.

>>> A = solver.jacobian(np.array([1, 1, 1]))
>>> B = solver.jacobian(np.array([1, 1, 1]))
>>> A is B
True

But passing a different p will trigger a computation and the cache will be replaced by the new value.

>>> C = solver.jacobian(np.array([1, 1, 1.1]))
>>> A is C
False
>>> np.all(A == C)
False
>>> D = solver.jacobian(np.array([1, 1, 1.1]))
>>> D is C
True

---