wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 14:24:21 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #38 from simpeg/develop

Browse Source 
Develop
This commit is contained in:
Rowan Cockett
2013-11-26 17:51:06 -08:00
parent 545b1637d8 83bdf73a63
commit 86e8952327
22 changed files with 3581 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yaml
+7
View File
@@ -0,0 +1,7 @@
language: python
python:
- "2.7"
# command to install dependencies
install: "pip install -r requirements.txt"
# command to run tests
script: nosetests
SimPEG/GCEtools/gceStartup.txt
+11
View File
@@ -0,0 +1,11 @@
# Check project status
gcutil getproject --project=<ProjectName> --cache_flag_values
# Start an instance
gcutil addinstance <instanceName>
# Log in
gcutil ssh <instanceName>
# Shut down
gcutil deleteinstance <instanceName>
SimPEG/GCEtools/startup.sh
+22
View File
@@ -0,0 +1,22 @@
#! /bin/bash
sudo aptitude -y update
sudo aptitude -y upgrade
sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev
sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip python-matplotlib
sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev
sudo aptitude -y install libblas-dev liblapack-dev
sudo pip install mpi4py
sudo pip install pymumps
sudo pip install scipy --upgrade
sudo pip install numpy --upgrade
sudo pip install ipython --upgrade
git clone https://bitbucket.org/rcockett/simpeg.git
cd simpeg/SimPEG/
python setup.py
cd ~
echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc
source .bashrc
SimPEG/forward/Problem.py
+2 -3
View File
@@ -226,7 +226,7 @@ class Problem(object):
"""
return sp.eye(m.size)
def createSyntheticData(self, m, std=0.05):
def createSyntheticData(self, m, std=0.05, u=None):
"""
Create synthetic data given a model, and a standard deviation.
@@ -238,8 +238,7 @@ class Problem(object):
Returns the observed data with random Gaussian noise
and Wd which is the same size as data, and can be used to weight the inversion.
"""
dobs = self.dpred(m)
dobs = dobs
dobs = self.dpred(m,u=u)
noise = std*abs(dobs)*np.random.randn(*dobs.shape)
dobs = dobs+noise
eps = np.linalg.norm(mkvc(dobs),2)*1e-5
SimPEG/inverse/BetaSchedule.py
+2 -2
View File
@@ -3,8 +3,8 @@
class Cooling(object):
"""Simple Beta Schedule"""
beta0 = 1.e6
beta_coolingFactor = 5.
beta0 = None #: The initial beta value, set to none means that it will be approximated in the first iteration.
beta_coolingFactor = 2.
def getBeta(self):
if self._beta is None:
SimPEG/inverse/Inversion.py
+136 -28
View File
@@ -1,19 +1,24 @@
import numpy as np
import scipy.sparse as sp
import SimPEG
from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt
from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks
from Optimize import Remember
from BetaSchedule import Cooling
from SimPEG.inverse import IterationPrinters, StoppingCriteria
class BaseInversion(object):
"""docstring for BaseInversion"""
maxIter = 1
maxIter = 1 #: Maximum number of iterations
name = 'BaseInversion'
debug = False
beta0 = 1e4
counter = None
debug = False #: Print debugging information
comment = '' #: Used by some functions to indicate what is going on in the algorithm
counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things
beta0 = None #: The initial Beta (regularization parameter)
def __init__(self, prob, reg, opt, **kwargs):
setKwargs(self, **kwargs)
@@ -22,18 +27,17 @@ class BaseInversion(object):
self.opt = opt
self.opt.parent = self
self.stoppers = [SimPEG.inverse.StoppingCriteria.iteration, SimPEG.inverse.StoppingCriteria.phi_d_target_Inversion]
self.stoppers = [StoppingCriteria.iteration]
# Check if we have inserted printers into the optimization
if not np.any([p is SimPEG.inverse.IterationPrinters.phi_d for p in self.opt.printers]):
self.opt.printers.insert(1,SimPEG.inverse.IterationPrinters.beta)
self.opt.printers.insert(2,SimPEG.inverse.IterationPrinters.phi_d)
self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m)
self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize)
if IterationPrinters.phi_d not in self.opt.printers:
self.opt.printers.insert(1,IterationPrinters.beta)
self.opt.printers.insert(2,IterationPrinters.phi_d)
self.opt.printers.insert(3,IterationPrinters.phi_m)
if not hasattr(opt, '_bfgsH0'): # Check if it has been set by the user and the default is not being used.
print 'Setting bfgsH0 to the inverse of the modelObj2Deriv.'
opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv(),doDirect=True,options={'factorize':True}) # False, options={'M':'GS','maxIter':15}
if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used.
print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv())
@property
@@ -68,9 +72,14 @@ class BaseInversion(object):
@timeIt
def run(self, m0):
"""run(m0)
Runs the inversion!
"""
self.startup(m0)
while True:
self._beta = self.getBeta()
self.doStartIteration()
self.m = self.opt.minimize(self.evalFunction, self.m)
self.doEndIteration()
if self.stoppingCriteria(): break
@@ -94,9 +103,7 @@ class BaseInversion(object):
:rtype: None
:return: None
"""
for method in [posible for posible in dir(self) if '_startup' in posible]:
if self.debug: print 'startup is calling self.'+method
getattr(self,method)(m0)
callHooks(self,'startup',m0)
if not hasattr(self.reg, '_mref'):
print 'Regularization has not set mref. SimPEG will set it to m0.'
@@ -108,25 +115,42 @@ class BaseInversion(object):
self.phi_d_last = np.nan
self.phi_m_last = np.nan
def doStartIteration(self):
"""
**doStartIteration** is called at the end of each run iteration.
If you have things that also need to run at the end of every iteration, you can create a method::
def _doStartIteration*(self):
pass
Where the * can be any string. If present, _doStartIteration* will be called at the start of the default doStartIteration call.
You may also completely overwrite this function.
:rtype: None
:return: None
"""
callHooks(self,'doStartIteration')
self._beta = self.getBeta()
def doEndIteration(self):
"""
**doEndIteration** is called at the end of each run iteration.
If you have things that also need to run at the end of every iteration, you can create a method::
def _doEndIteration*(self, xt):
def _doEndIteration*(self):
pass
Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call.
You may also completely overwrite this function.
:param numpy.ndarray xt: tested new iterate that ensures a descent direction.
:rtype: None
:return: None
"""
for method in [posible for posible in dir(self) if '_doEndIteration' in posible]:
if self.debug: print 'doEndIteration is calling self.'+method
getattr(self,method)()
callHooks(self,'doEndIteration')
# store old values
self.phi_d_last = self.phi_d
@@ -136,6 +160,47 @@ class BaseInversion(object):
def getBeta(self):
return self.beta0
def estimateBeta0(self, u=None, ratio=0.1):
"""estimateBeta0(u=None, ratio=0.1)
The initial beta is calculated by comparing the estimated
eigenvalues of JtJ and WtW.
To estimate the eigenvector of **A**, we will use one iteration
of the *Power Method*:
.. math::
\mathbf{x_1 = A x_0}
Given this (very course) approximation of the eigenvector,
we can use the *Rayleigh quotient* to approximate the largest eigenvalue.
.. math::
\lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}}
We will approximate the largest eigenvalue for both JtJ and WtW, and
use some ratio of the quotient to estimate beta0.
.. math::
\\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}}
:param numpy.array u: fields
:param float ratio: desired ratio of the eigenvalues, default is 0.1
:rtype: float
:return: beta0
"""
if u is None:
u = self.prob.field(self.m)
x0 = np.random.rand(*self.m.shape)
t = x0.dot(self.dataObj2Deriv(self.m,x0,u=u))
b = x0.dot(self.reg.modelObj2Deriv()*x0)
return ratio*(t/b)
def stoppingCriteria(self):
if self.debug: print 'checking stoppingCriteria'
return checkStoppers(self, self.stoppers)
@@ -150,11 +215,20 @@ class BaseInversion(object):
@timeIt
def evalFunction(self, m, return_g=True, return_H=True):
"""evalFunction(m, return_g=True, return_H=True)
"""
u = self.prob.field(m)
if self._iter is 0 and self._beta is None:
self._beta = self.beta0 = self.estimateBeta0(u=u)
phi_d = self.dataObj(m, u)
phi_m = self.reg.modelObj(m)
self.dpred = self.prob.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.phi_d = phi_d
self.phi_m = phi_m
@@ -175,13 +249,14 @@ class BaseInversion(object):
return phi_d2Deriv + self._beta * phi_m2Deriv
operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype )
out += (operator,)
return out if len(out) > 1 else out[0]
@timeIt
def dataObj(self, m, u=None):
"""
"""dataObj(m, u=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: float
@@ -203,7 +278,8 @@ class BaseInversion(object):
@timeIt
def dataObjDeriv(self, m, u=None):
"""
"""dataObjDeriv(m, u=None)
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: numpy.array
@@ -244,8 +320,10 @@ class BaseInversion(object):
@timeIt
def dataObj2Deriv(self, m, v, u=None):
"""
"""dataObj2Deriv(m, v, u=None)
:param numpy.array m: geophysical model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: data misfit derivative
@@ -282,7 +360,7 @@ class BaseInversion(object):
R = self.Wd*self.prob.dataResidual(m, u=u)
# TODO: abstract to different norms a little cleaner.
# \/ it goes here. in l2 it is the identity.
# \/ it goes here. in l2 it is the identity.
dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u)
return dmisfit
@@ -294,3 +372,33 @@ class Inversion(Cooling, Remember, BaseInversion):
def __init__(self, prob, reg, opt, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, **kwargs)
self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
class TimeSteppingInversion(Remember, BaseInversion):
"""
A slightly different view on regularization parameters,
let Beta be viewed as 1/dt, and timestep by updating the
reference model every optimization iteration.
"""
maxIter = 1
name = "Time-Stepping SimPEG Inversion"
def __init__(self, prob, reg, opt, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, **kwargs)
self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
def _startup_TimeSteppingInversion(self, m0):
def _doEndIteration_updateMref(self, xt):
if self.debug: 'Updating the reference model.'
self.parent.reg.mref = self.xc
self.opt.hook(_doEndIteration_updateMref, overwrite=True)
SimPEG/inverse/Optimize.py
+77 -91
View File
@@ -1,16 +1,13 @@
import numpy as np
import matplotlib.pyplot as plt
from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt
from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt, callHooks
norm = np.linalg.norm
import scipy.sparse as sp
from SimPEG import Solver
try:
from pubsub import pub
doPub = True
except Exception, e:
print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.'
doPub = False
__all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters']
class StoppingCriteria(object):
"""docstring for StoppingCriteria"""
@@ -85,32 +82,28 @@ class IterationPrinters(object):
class Minimize(object):
"""
Minimize is a general class for derivative based optimization.
"""
name = "General Optimization Algorithm"
name = "General Optimization Algorithm" #: The name of the optimization algorithm
maxIter = 20
maxIterLS = 10
maxStep = np.inf
LSreduction = 1e-4
LSshorten = 0.5
tolF = 1e-1
tolX = 1e-1
tolG = 1e-1
eps = 1e-5
maxIter = 20 #: Maximum number of iterations
maxIterLS = 10 #: Maximum number of iterations for the line-search
maxStep = np.inf #: Maximum step possible, used in scaling before the line-search.
LSreduction = 1e-4 #: Expected decrease in the line-search
LSshorten = 0.5 #: Line-search step is shortened by this amount each time.
tolF = 1e-1 #: Tolerance on function value decrease
tolX = 1e-1 #: Tolerance on norm(x) movement
tolG = 1e-1 #: Tolerance on gradient norm
eps = 1e-5 #: Small value
debug = False
debugLS = False
debug = False #: Print debugging information
debugLS = False #: Print debugging information for the line-search
comment = ''
counter = None
comment = '' #: Used by some functions to indicate what is going on in the algorithm
counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things
def __init__(self, **kwargs):
self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub
self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration]
self.stoppersLS = [StoppingCriteria.armijoGoldstein, StoppingCriteria.iterationLS]
@@ -121,7 +114,8 @@ class Minimize(object):
@timeIt
def minimize(self, evalFunction, x0):
"""
"""minimize(evalFunction, x0)
Minimizes the function (evalFunction) starting at the location x0.
:param def evalFunction: function handle that evaluates: f, g, H = F(x)
@@ -142,33 +136,13 @@ class Minimize(object):
return out if len(out) > 1 else out[0]
Events are fired with the following inputs via pypubsub::
Minimize.printInit (minimize)
Minimize.evalFunction (minimize, f, g, H)
Minimize.printIter (minimize)
Minimize.searchDirection (minimize, p)
Minimize.scaleSearchDirection (minimize, p)
Minimize.modifySearchDirection (minimize, xt, passLS)
Minimize.endIteration (minimize, xt)
Minimize.printDone (minimize)
To hook into one of these events (must have pypubsub installed)::
from pubsub import pub
def listener(minimize,p):
print 'The search direction is: ', p
pub.subscribe(listener, 'Minimize.searchDirection')
You can use pubsub communication to debug your code, it is not used internally.
The algorithm for general minimization is as follows::
startup(x0)
printInit()
while True:
doStartIteration()
f, g, H = evalFunction(xc)
printIter()
if stoppingCriteria(): break
@@ -188,21 +162,17 @@ class Minimize(object):
self.printInit()
while True:
self.doStartIteration()
self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True)
if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H)
self.printIter()
if self.stoppingCriteria(): break
p = self.findSearchDirection()
if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p)
p = self.scaleSearchDirection(p)
if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p)
self.searchDirection = self.findSearchDirection()
p = self.scaleSearchDirection(self.searchDirection)
xt, passLS = self.modifySearchDirection(p)
if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt, passLS=passLS)
if not passLS:
xt, caught = self.modifySearchDirectionBreak(p)
if not caught: return self.xc
self.doEndIteration(xt)
if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt)
self.printDone()
@@ -240,9 +210,7 @@ class Minimize(object):
:rtype: None
:return: None
"""
for method in [posible for posible in dir(self) if '_startup' in posible]:
if self.debug: print 'startup is calling self.'+method
getattr(self,method)(x0)
callHooks(self,'startup',x0)
self._iter = 0
self._iterLS = 0
@@ -253,6 +221,16 @@ class Minimize(object):
self.f_last = np.nan
self.x_last = x0
@count
def doStartIteration(self):
"""doStartIteration()
**doStartIteration** is called at the start of each minimize iteration.
:rtype: None
:return: None
"""
callHooks(self,'doStartIteration')
def printInit(self, inLS=False):
"""
@@ -262,12 +240,10 @@ class Minimize(object):
parent.printInit function and call that.
"""
if doPub and not inLS: pub.sendMessage('Minimize.printInit', minimize=self)
pad = ' '*10 if inLS else ''
name = self.name if not inLS else self.nameLS
printTitles(self, self.printers if not inLS else self.printersLS, name, pad)
@count
def printIter(self, inLS=False):
"""
**printIter** is called directly after function evaluations.
@@ -276,12 +252,8 @@ class Minimize(object):
parent.printIter function and call that.
"""
callHooks(self,'printIter',inLS)
for method in [posible for posible in dir(self) if '_printIter' in posible]:
if self.debug: print 'printIter is calling self.'+method
getattr(self,method)(inLS)
if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self)
pad = ' '*10 if inLS else ''
printLine(self, self.printers if not inLS else self.printersLS, pad=pad)
@@ -293,13 +265,11 @@ class Minimize(object):
parent.printDone function and call that.
"""
if doPub and not inLS: pub.sendMessage('Minimize.printDone', minimize=self)
pad = ' '*10 if inLS else ''
stop, done = (' STOP! ', ' DONE! ') if not inLS else ('----------------', ' End Linesearch ')
stoppers = self.stoppers if not inLS else self.stoppersLS
printStoppers(self, stoppers, pad='', stop=stop, done=done)
@timeIt
def stoppingCriteria(self, inLS=False):
if self._iter == 0:
self.f0 = self.f
@@ -308,7 +278,8 @@ class Minimize(object):
@timeIt
def projection(self, p):
"""
"""projection(p)
projects the search direction.
by default, no projection is applied.
@@ -317,11 +288,13 @@ class Minimize(object):
:rtype: numpy.ndarray
:return: p, projected search direction
"""
callHooks(self,'projection',p)
return p
@timeIt
def findSearchDirection(self):
"""
"""findSearchDirection()
**findSearchDirection** should return an approximation of:
.. math::
@@ -351,7 +324,8 @@ class Minimize(object):
@count
def scaleSearchDirection(self, p):
"""
"""scaleSearchDirection(p)
**scaleSearchDirection** should scale the search direction if appropriate.
Set the parameter **maxStep** in the minimize object, to scale back the gradient to a maximum size.
@@ -365,11 +339,12 @@ class Minimize(object):
p = self.maxStep*p/np.abs(p.max())
return p
nameLS = "Armijo linesearch"
nameLS = "Armijo linesearch" #: The line-search name
@timeIt
def modifySearchDirection(self, p):
"""
"""modifySearchDirection(p)
**modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
By default, an Armijo backtracking linesearch is preformed with the following parameters:
@@ -406,7 +381,8 @@ class Minimize(object):
@count
def modifySearchDirectionBreak(self, p):
"""
"""modifySearchDirectionBreak(p)
Code is called if modifySearchDirection fails
to find a descent direction.
@@ -427,7 +403,8 @@ class Minimize(object):
@count
def doEndIteration(self, xt):
"""
"""doEndIteration(xt)
**doEndIteration** is called at the end of each minimize iteration.
By default, function values and x locations are shuffled to store 1 past iteration in memory.
@@ -447,9 +424,7 @@ class Minimize(object):
:rtype: None
:return: None
"""
for method in [posible for posible in dir(self) if '_doEndIteration' in posible]:
if self.debug: print 'doEndIteration is calling self.'+method
getattr(self,method)(xt)
callHooks(self,'doEndIteration',xt)
# store old values
self.f_last = self.f
@@ -458,7 +433,6 @@ class Minimize(object):
if self.debug: self.printDone()
class Remember(object):
"""
This mixin remembers all the things you tend to forget.
@@ -509,12 +483,11 @@ class Remember(object):
self._rememberList[param[0]].append( param[1](self) )
class ProjectedGradient(Minimize, Remember):
name = 'Projected Gradient'
maxIterCG = 10
tolCG = 1e-3
maxIterCG = 5
tolCG = 1e-1
lower = -np.inf
upper = np.inf
@@ -546,22 +519,35 @@ class ProjectedGradient(Minimize, Remember):
@count
def projection(self, x):
"""Make sure we are feasible."""
"""projection(x)
Make sure we are feasible.
"""
return np.median(np.c_[self.lower,x,self.upper],axis=1)
@count
def activeSet(self, x):
"""If we are on a bound"""
"""activeSet(x)
If we are on a bound
"""
return np.logical_or(x == self.lower, x == self.upper)
@count
def inactiveSet(self, x):
"""The free variables."""
"""inactiveSet(x)
The free variables.
"""
return np.logical_not(self.activeSet(x))
@count
def bindingSet(self, x):
"""
"""bindingSet(x)
If we are on a bound and the negative gradient points away from the feasible set.
Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91
@@ -573,6 +559,10 @@ class ProjectedGradient(Minimize, Remember):
@timeIt
def findSearchDirection(self):
"""findSearchDirection()
Finds the search direction based on either CG or steepest descent.
"""
self.aSet_prev = self.activeSet(self.xc)
allBoundsAreActive = sum(self.aSet_prev) == self.xc.size
@@ -614,6 +604,7 @@ class ProjectedGradient(Minimize, Remember):
@timeIt
def _doEndIteration_ProjectedGradient(self, xt):
"""_doEndIteration_ProjectedGradient(xt)"""
aSet = self.activeSet(xt)
bSet = self.bindingSet(xt)
@@ -642,7 +633,6 @@ class ProjectedGradient(Minimize, Remember):
if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
class BFGS(Minimize, Remember):
name = 'BFGS'
nbfgs = 10
@@ -738,8 +728,8 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
name = 'Inexact Gauss Newton'
maxIterCG = 10
tolCG = 1e-3
maxIterCG = 5
tolCG = 1e-1
@property
def approxHinv(self):
@@ -857,10 +847,6 @@ if __name__ == '__main__':
x0 = np.array([2.6, 3.7])
checkDerivative(Rosenbrock, x0, plotIt=False)
# def listener1(minimize,p):
# print 'hi: ', p
# if doPub: pub.subscribe(listener1, 'Minimize.searchDirection')
xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0)
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0)
SimPEG/mesh/DiffOperators.py
+71 -6
View File
@@ -152,6 +152,73 @@ class DiffOperators(object):
_faceDiv = None
faceDiv = property(**faceDiv())
def faceDivx():
doc = "Construct divergence operator in the x component (face-stg to cell-centres)."
def fget(self):
if(self._faceDivx is None):
# The number of cell centers in each direction
n = self.n
# Compute faceDivergence operator on faces
if(self.dim == 1):
D1 = ddx(n[0])
elif(self.dim == 2):
D1 = sp.kron(speye(n[1]), ddx(n[0]))
elif(self.dim == 3):
D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fx', 'V')
V = self.vol
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
return self._faceDivx
return locals()
_faceDivx = None
faceDivx = property(**faceDivx())
def faceDivy():
doc = "Construct divergence operator in the y component (face-stg to cell-centres)."
def fget(self):
if(self.dim < 2): return None
if(self._faceDivy is None):
# The number of cell centers in each direction
n = self.n
# Compute faceDivergence operator on faces
if(self.dim == 2):
D2 = sp.kron(ddx(n[1]), speye(n[0]))
elif(self.dim == 3):
D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fy', 'V')
V = self.vol
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
return self._faceDivy
return locals()
_faceDivy = None
faceDivy = property(**faceDivy())
def faceDivz():
doc = "Construct divergence operator in the z component (face-stg to cell-centres)."
def fget(self):
if(self.dim < 3): return None
if(self._faceDivz is None):
# The number of cell centers in each direction
n = self.n
# Compute faceDivergence operator on faces
D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fz', 'V')
V = self.vol
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
return self._faceDivz
return locals()
_faceDivz = None
faceDivz = property(**faceDivz())
def nodalGrad():
doc = "Construct gradient operator (nodes to edges)."
@@ -279,12 +346,12 @@ class DiffOperators(object):
elif(self.dim == 2):
G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
G = sp.vstack((G1, G2), format="csr")
G = sp.block_diag((G1, G2), format="csr")
elif(self.dim == 3):
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
G = sp.vstack((G1, G2, G3), format="csr")
G = sp.block_diag((G1, G2, G3), format="csr")
# Compute areas of cell faces & volumes
S = self.area
V = self.aveCC2F*self.vol # Average volume between adjacent cells
@@ -318,8 +385,7 @@ class DiffOperators(object):
def cellGrady():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
if self.dim < 2:
return None
if self.dim < 2: return None
if getattr(self, '_cellGrady', None) is None:
BC = ['neumann', 'neumann']
n = self.n
@@ -338,8 +404,7 @@ class DiffOperators(object):
def cellGradz():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
if self.dim < 3:
return None
if self.dim < 3: return None
if getattr(self, '_cellGradz', None) is None:
BC = ['neumann', 'neumann']
n = self.n
SimPEG/mesh/TensorMesh.py
+5 -3
View File
@@ -33,14 +33,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
"""
_meshType = 'TENSOR'
def __init__(self, h, x0=None):
for i, h_i in enumerate(h):
def __init__(self, h_in, x0=None):
assert type(h_in) is list, 'h_in must be a list'
h = range(len(h_in))
for i, h_i in enumerate(h_in):
if type(h_i) in [int, long, float]:
# This gives you something over the unit cube.
h_i = np.ones(int(h_i))/int(h_i)
h[i] = h_i
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
h[i] = h_i[:] # make a copy.
BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
SimPEG/mesh/TensorView.py
+67 -4
View File
@@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from SimPEG.utils import mkvc
from SimPEG.utils import mkvc, animate
class TensorView(object):
@@ -14,7 +14,7 @@ class TensorView(object):
def __init__(self):
pass
def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False):
def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False,clim=None):
"""
Mesh.plotImage(I)
@@ -141,7 +141,9 @@ class TensorView(object):
C = I[:].reshape(self.nEy, order='F')
C = 0.5*(C[:-1,:] + C[1:,:] )
ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T)
if clim is None:
clim = [C.min(),C.max()]
ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T, vmin=clim[0], vmax=clim[1])
ax.axis('tight')
ax.set_xlabel("x")
ax.set_ylabel("y")
@@ -196,7 +198,10 @@ class TensorView(object):
xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())]
yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())]
# Plot the mesh
ph = ax.pcolormesh(xx, yy, C.T)
if clim is None:
clim = [C.min(),C.max()]
ph = ax.pcolormesh(xx, yy, C.T, vmin=clim[0], vmax=clim[1])
# Plot the lines
gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0])
gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1])
@@ -336,3 +341,61 @@ class TensorView(object):
ax.set_ylabel('x2')
ax.set_zlabel('x3')
if showIt: plt.show()
def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):
assert normal in 'xyz', 'normal must be x, y, or z'
if ax is None: ax = plt.subplot(111)
I = mesh.r(var,'CC','CC','M')
axes = [p for p in 'xyz' if p not in normal.lower()]
if normal is 'x': I = I[index,:,:]
if normal is 'y': I = I[:,index,:]
if normal is 'z': I = I[:,:,index]
if clim is None: clim = [I.min(),I.max()]
p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1])
ax.axis('tight')
ax.set_xlabel(axes[0])
ax.set_ylabel(axes[1])
return p
def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)):
assert mesh.dim > 2, 'This is for 3D meshes only.'
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=figsize)
ax = plt.axes()
clim = [var.min(),var.max()]
plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim))
tlt = plt.title(normal)
def animateFrame(i):
mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)])
def video(mesh,var,function,figsize=(10,8),colorbar=True):
"""
Call a function for a list of models to create a video.
::
def function(var, ax, clim, tlt, i):
tlt.set_text('%%d'%%i)
return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
mesh.video([model1, model2, ..., modeln],function)
"""
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=figsize)
ax = plt.axes()
VAR = np.concatenate(var)
clim = [VAR.min(),VAR.max()]
tlt = plt.title('')
if colorbar:
plt.colorbar(function(var[0],ax,clim,tlt,0))
def animateFrame(i):
function(var[i],ax,clim,tlt,i)
return animate(fig, animateFrame, frames=len(var))
SimPEG/regularization/Regularization.py
+10 -11
View File
@@ -1,4 +1,4 @@
from SimPEG.utils import sdiag, count, timeIt
from SimPEG.utils import sdiag, count, timeIt, setKwargs
import numpy as np
class Regularization(object):
@@ -22,34 +22,33 @@ class Regularization(object):
@property
def Wx(self):
if getattr(self, '_Wx', None) is None:
a = self.mesh.r(self.mesh.area,'F','Fx','V')
self._Wx = sdiag(a)*self.mesh.cellGradx
self._Wx = self.mesh.cellGradx*sdiag(self.mesh.vol)
return self._Wx
@property
def Wy(self):
if getattr(self, '_Wy', None) is None:
a = self.mesh.r(self.mesh.area,'F','Fy','V')
self._Wy = sdiag(a)*self.mesh.cellGrady
self._Wy = self.mesh.cellGrady*sdiag(self.mesh.vol)
return self._Wy
@property
def Wz(self):
if getattr(self, '_Wz', None) is None:
a = self.mesh.r(self.mesh.area,'F','Fz','V')
self._Wz = sdiag(a)*self.mesh.cellGradz
self._Wz = self.mesh.cellGradz*sdiag(self.mesh.vol)
return self._Wz
alpha_s = 1e-6
alpha_x = 1
alpha_y = 1
alpha_z = 1
alpha_x = 1.0
alpha_y = 1.0
alpha_z = 1.0
counter = None
def __init__(self, mesh):
def __init__(self, mesh, **kwargs):
setKwargs(self, **kwargs)
self.mesh = mesh
def pnorm(self, r):
return 0.5*r.dot(r)
SimPEG/utils/Solver.py
+67
View File
@@ -19,6 +19,11 @@ except Exception, e:
DEFAULTS['forward'] = 'python'
DEFAULTS['backward'] = 'python'
try:
import mumps
except Exception, e:
print 'Warning: mumps solver not available.'
class Solver(object):
"""
Solver is a light wrapper on the various types of
@@ -113,6 +118,9 @@ class Solver(object):
def clean(self):
"""Cleans up the memory"""
if self.options.has_key('backend'):
if self.options['backend'] == 'mumps':
self.mctx.destroy()
del self.dsolve
self.dsolve = None
@@ -120,6 +128,7 @@ class Solver(object):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:param str backend: which backend to use. Default is scipy
:rtype: numpy.ndarray
@@ -129,6 +138,22 @@ class Solver(object):
assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
if backend == 'scipy':
X = self.solveDirect_scipy(b, factorize)
elif backend == 'mumps':
X = self.solveDirect_mumps(b, factorize)
return X
def solveDirect_scipy(self, b, factorize):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:rtype: numpy.ndarray
:return: x
"""
if factorize and self.dsolve is None:
self.A = self.A.tocsc() # for efficiency
self.dsolve = linalg.factorized(self.A)
@@ -150,6 +175,48 @@ class Solver(object):
return X
def solveDirect_mumps(self, b, factorize):
"""
Use solve instead of this interface.
:param numpy.ndarray b: the right hand side
:param bool factorize: if you want to factorize and store factors
:rtype: numpy.ndarray
:return: x
"""
if factorize and self.dsolve is None:
self.mctx = mumps.DMumpsContext()
self.mctx.set_icntl(14, 60)
# self.mctx.set_silent()
self.mctx.set_centralized_sparse(self.A)
self.mctx.run(job=4)
def mdsolve(rhs):
x = rhs.copy()
self.mctx.set_rhs(x)
self.mctx.run(job=3)
return x
self.dsolve = mdsolve
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
if factorize:
X = self.dsolve(b)
else:
X = mumps.spsolve(self.A, b)
else:
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
if factorize:
X[:,i] = self.dsolve(b[:,i])
else:
X[:,i] = mumps.spsolve(self.A,b[:,i])
return X
def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
if backend is None: backend = DEFAULTS['iter']
SimPEG/utils/__init__.py
+35 -3
View File
@@ -12,6 +12,29 @@ import Solver
from Solver import Solver
import Geophysics
import types
import time
import numpy as np
from functools import wraps
def hook(obj, method, name=None, overwrite=False, silent=False):
"""
This dynamically binds a method to the instance of the class.
If name is None, the name of the method is used.
"""
if name is None:
name = method.__name__
if name == '<lambda>':
raise Exception('Must provide name to hook lambda functions.')
if not hasattr(obj,name) or overwrite:
setattr(obj, name, types.MethodType( method, obj ))
if getattr(obj,'debug',False):
print 'Method '+name+' was added to class.'
elif not silent or getattr(obj,'debug',False):
print 'Method '+name+' was not overwritten.'
def setKwargs(obj, **kwargs):
"""Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
for attr in kwargs:
@@ -19,6 +42,9 @@ def setKwargs(obj, **kwargs):
setattr(obj, attr, kwargs[attr])
else:
raise Exception('%s attr is not recognized' % attr)
hook(obj,callHooks, silent=True)
hook(obj,hook, silent=True)
hook(obj,setKwargs, silent=True)
def printTitles(obj, printers, name='Print Titles', pad=''):
titles = ''
@@ -61,9 +87,11 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
print pad + stopper['str'] % (l<=r,l,r)
print pad + "%s%s%s" % ('-'*25,done,'-'*25)
def callHooks(obj, match, *args, **kwargs):
for method in [posible for posible in dir(obj) if ('_'+match) in posible]:
if getattr(obj,'debug',False): print (match+' is calling self.'+method)
getattr(obj,method)(*args, **kwargs)
import time
import numpy as np
class Counter(object):
@@ -75,7 +103,7 @@ class Counter(object):
If you want to use this, import *count* or *timeIt* and use them as decorators on class methods.
.. ::
::
class MyClass(object):
def __init__(self, url):
@@ -139,6 +167,7 @@ class Counter(object):
print " {0:<40}: {1:4.2e}, {2:4.2e}, {3:4d}x".format(prop,a.mean(),a.sum(),l)
def count(f):
@wraps(f)
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
@@ -147,6 +176,7 @@ def count(f):
return wrapper
def timeIt(f):
@wraps(f)
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
@@ -154,6 +184,8 @@ def timeIt(f):
if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
return out
return wrapper
if __name__ == '__main__':
class MyClass(object):
def __init__(self, url):
docs/api_LogicallyOrthogonalMesh.rst
+2
View File
@@ -4,8 +4,10 @@ Logically Orthogonal Mesh
*************************
.. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
LOM View
docs/api_Optimize.rst
+3
View File
@@ -4,6 +4,8 @@ Optimize
********
.. automodule:: SimPEG.inverse.Optimize
:show-inheritance:
:private-members:
:members:
:undoc-members:
@@ -12,6 +14,7 @@ Inversion
*********
.. automodule:: SimPEG.inverse.Inversion
:show-inheritance:
:members:
:undoc-members:
docs/api_Problem.rst
+6 -3
View File
@@ -1,21 +1,24 @@
.. _api_Problem:
Problem
*******
.. automodule:: SimPEG.forward.Problem
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
DCProblem
*********
.. automodule:: SimPEG.forward.DCProblem
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
@@ -23,7 +26,7 @@ Linear Problem
**************
.. automodule:: SimPEG.forward.LinearProblem
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
docs/api_TensorMesh.rst
+2
View File
@@ -4,8 +4,10 @@ Tensor Mesh
***********
.. automodule:: SimPEG.mesh.TensorMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
Tensor View
***********
docs/index.rst
+4 -2
View File
@@ -1,5 +1,7 @@
SimPEG
======
.. image:: simpeg-logo.png
:width: 300 px
:alt: SimPEG
:align: center
SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the context of geophysical applications.
docs/requirements.txt
+1 -1
View File
@@ -1,3 +1,3 @@
numpy
pypubsub
scipy
ipython
docs/simpeg-logo.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 23 KiB

notebooks/SliceThroughModel.ipynb
+3048
View File
File diff suppressed because it is too large Load Diff
requirements.txt
+3
View File
@@ -0,0 +1,3 @@
numpy
pypubsub
ipython