mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-08 22:55:26 +08:00
Merged in updateFramework (pull request #28)
BFGS, Iterative Solvers, and Updates to framework.
This commit is contained in:
rowanc1
@@ -1,5 +1,5 @@
|
||||
from SimPEG.mesh import TensorMesh
|
||||
from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms
|
||||
from SimPEG.forward import Problem, ModelTransforms
|
||||
from SimPEG.tests import checkDerivative
|
||||
from SimPEG.utils import ModelBuilder, sdiag, mkvc
|
||||
from SimPEG import Solver
|
||||
@@ -201,23 +201,17 @@ if __name__ == '__main__':
|
||||
P = Q.T
|
||||
|
||||
# Create some data
|
||||
class syntheticDCProblem(DCProblem, SyntheticProblem):
|
||||
pass
|
||||
problem = DCProblem(mesh)
|
||||
problem.P = P
|
||||
problem.RHS = q
|
||||
dobs, Wd = problem.createSyntheticData(mSynth, std=0.05)
|
||||
|
||||
synthetic = syntheticDCProblem(mesh);
|
||||
synthetic.P = P
|
||||
synthetic.RHS = q
|
||||
dobs, Wd = synthetic.createData(mSynth, std=0.05)
|
||||
|
||||
u = synthetic.field(mSynth)
|
||||
u = synthetic.reshapeFields(u)
|
||||
u = problem.field(mSynth)
|
||||
u = problem.reshapeFields(u)
|
||||
mesh.plotImage(u[:,10])
|
||||
# plt.show()
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
problem = DCProblem(mesh)
|
||||
problem.P = P
|
||||
problem.RHS = q
|
||||
problem.dobs = dobs
|
||||
problem.std = dobs*0 + 0.05
|
||||
m0 = mesh.gridCC[:,0]*0+sig2
|
||||
|
||||
@@ -226,27 +226,15 @@ class Problem(object):
|
||||
"""
|
||||
return sp.eye(m.size)
|
||||
|
||||
|
||||
|
||||
|
||||
class SyntheticProblem(object):
|
||||
"""
|
||||
Has helpful functions when dealing with synthetic problems
|
||||
|
||||
To use this class, inherit to your problem::
|
||||
|
||||
class mySyntheticExample(Problem, SyntheticProblem):
|
||||
pass
|
||||
"""
|
||||
def createData(self, m, std=0.05):
|
||||
def createSyntheticData(self, m, std=0.05):
|
||||
"""
|
||||
Create synthetic data given a model, and a standard deviation.
|
||||
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array std: standard deviation
|
||||
:rtype: numpy.array, numpy.array
|
||||
:return: dobs, Wd
|
||||
|
||||
Create synthetic data given a model, and a standard deviation.
|
||||
|
||||
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.
|
||||
"""
|
||||
|
||||
@@ -31,6 +31,11 @@ class BaseInversion(object):
|
||||
self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m)
|
||||
self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize)
|
||||
|
||||
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}
|
||||
|
||||
|
||||
@property
|
||||
def Wd(self):
|
||||
"""
|
||||
@@ -40,6 +45,9 @@ class BaseInversion(object):
|
||||
eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5
|
||||
self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
|
||||
return self._Wd
|
||||
@Wd.setter
|
||||
def Wd(self, value):
|
||||
self._Wd = value
|
||||
|
||||
@property
|
||||
def phi_d_target(self):
|
||||
@@ -90,9 +98,15 @@ class BaseInversion(object):
|
||||
if self.debug: print 'startup is calling self.'+method
|
||||
getattr(self,method)(m0)
|
||||
|
||||
if not hasattr(self.reg, '_mref'):
|
||||
print 'Regularization has not set mref. SimPEG will set it to m0.'
|
||||
self.reg.mref = m0
|
||||
|
||||
self.m = m0
|
||||
self._iter = 0
|
||||
self._beta = None
|
||||
self.phi_d_last = np.nan
|
||||
self.phi_m_last = np.nan
|
||||
|
||||
def doEndIteration(self):
|
||||
"""
|
||||
@@ -157,7 +171,7 @@ class BaseInversion(object):
|
||||
if return_H:
|
||||
def H_fun(v):
|
||||
phi_d2Deriv = self.dataObj2Deriv(m, v, u=u)
|
||||
phi_m2Deriv = self.reg.modelObj2Deriv(m)*v
|
||||
phi_m2Deriv = self.reg.modelObj2Deriv()*v
|
||||
|
||||
return phi_d2Deriv + self._beta * phi_m2Deriv
|
||||
|
||||
|
||||
+204
-9
@@ -76,7 +76,7 @@ class IterationPrinters(object):
|
||||
itType = {"title": "itType", "value": lambda M: M._itType, "width": 8, "format": "%s"}
|
||||
aSet = {"title": "aSet", "value": lambda M: np.sum(M.activeSet(M.xc)), "width": 8, "format": "%d"}
|
||||
bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"}
|
||||
comment = {"title": "Comment", "value": lambda M: M.projComment, "width": 7, "format": "%s"}
|
||||
comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"}
|
||||
|
||||
beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format": "%1.2e"}
|
||||
phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"}
|
||||
@@ -106,6 +106,7 @@ class Minimize(object):
|
||||
debug = False
|
||||
debugLS = False
|
||||
|
||||
comment = ''
|
||||
counter = None
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
@@ -275,6 +276,11 @@ class Minimize(object):
|
||||
parent.printIter function and call that.
|
||||
|
||||
"""
|
||||
|
||||
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)
|
||||
@@ -534,7 +540,7 @@ class ProjectedGradient(Minimize, Remember):
|
||||
self.stopDoingPG = False
|
||||
|
||||
self._itType = 'SD'
|
||||
self.projComment = ''
|
||||
self.comment = ''
|
||||
|
||||
self.aSet_prev = self.activeSet(x0)
|
||||
|
||||
@@ -600,7 +606,7 @@ class ProjectedGradient(Minimize, Remember):
|
||||
def reduceHess(v):
|
||||
# Z is tall and skinny
|
||||
return Z.T*(self.H*(Z*v))
|
||||
operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=float )
|
||||
operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=self.xc.dtype )
|
||||
p, info = sp.linalg.cg(operator, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG)
|
||||
p = Z*p # bring up to full size
|
||||
# aSet_after = self.activeSet(self.xc+p)
|
||||
@@ -615,7 +621,7 @@ class ProjectedGradient(Minimize, Remember):
|
||||
self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient
|
||||
|
||||
f_current_decrease = self.f_last - self.f
|
||||
self.projComment = ''
|
||||
self.comment = ''
|
||||
if self._iter < 1:
|
||||
# Note that this is reset on every CG iteration.
|
||||
self.f_decrease_max = -np.inf
|
||||
@@ -623,7 +629,7 @@ class ProjectedGradient(Minimize, Remember):
|
||||
self.f_decrease_max = max(self.f_decrease_max, f_current_decrease)
|
||||
self.stopDoingPG = f_current_decrease < 0.25 * self.f_decrease_max
|
||||
if self.stopDoingPG:
|
||||
self.projComment = 'Stop SD'
|
||||
self.comment = 'Stop SD'
|
||||
self.explorePG = False
|
||||
self.exploreCG = True
|
||||
# implement 3.8, MoreToraldo91
|
||||
@@ -635,6 +641,74 @@ class ProjectedGradient(Minimize, Remember):
|
||||
if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max
|
||||
if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
|
||||
|
||||
|
||||
|
||||
class BFGS(Minimize, Remember):
|
||||
name = 'BFGS'
|
||||
nbfgs = 10
|
||||
|
||||
@property
|
||||
def bfgsH0(self):
|
||||
"""
|
||||
Approximate Hessian used in preconditioning the problem.
|
||||
|
||||
Must be a SimPEG.Solver
|
||||
"""
|
||||
_bfgsH0 = getattr(self,'_bfgsH0',None)
|
||||
if _bfgsH0 is None:
|
||||
return Solver(sp.identity(self.xc.size).tocsc(), flag='D')
|
||||
return _bfgsH0
|
||||
@bfgsH0.setter
|
||||
def bfgsH0(self, value):
|
||||
assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver'
|
||||
self._bfgsH0 = value
|
||||
|
||||
def _startup_BFGS(self,x0):
|
||||
self._bfgscnt = -1
|
||||
self._bfgsY = np.zeros((x0.size, self.nbfgs))
|
||||
self._bfgsS = np.zeros((x0.size, self.nbfgs))
|
||||
if not np.any([p is IterationPrinters.comment for p in self.printers]):
|
||||
self.printers.append(IterationPrinters.comment)
|
||||
|
||||
def bfgs(self, d):
|
||||
n = self._bfgscnt
|
||||
nn = ktop = min(self._bfgsS.shape[1],n)
|
||||
return self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d)
|
||||
|
||||
def bfgsrec(self,k,n,nn,S,Y,d):
|
||||
"""BFGS recursion"""
|
||||
if k < 0:
|
||||
d = self.bfgsH0.solve(d)
|
||||
else:
|
||||
khat = 0 if nn is 0 else np.mod(n-nn+k,nn)
|
||||
gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])
|
||||
d = d - gamma*Y[:,khat]
|
||||
d = self.bfgsrec(k-1,n,nn,S,Y,d)
|
||||
d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]
|
||||
return d
|
||||
|
||||
def findSearchDirection(self):
|
||||
return self.bfgs(-self.g)
|
||||
|
||||
def _doEndIteration_BFGS(self, xt):
|
||||
if self._iter is 0:
|
||||
self.g_last = self.g
|
||||
return
|
||||
|
||||
yy = self.g - self.g_last;
|
||||
ss = self.xc - xt;
|
||||
self.g_last = self.g
|
||||
|
||||
if yy.dot(ss) > 0:
|
||||
self._bfgscnt += 1
|
||||
ktop = np.mod(self._bfgscnt,self.nbfgs)
|
||||
self._bfgsY[:,ktop] = yy
|
||||
self._bfgsS[:,ktop] = ss
|
||||
self.comment = ''
|
||||
else:
|
||||
self.comment = 'Skip BFGS'
|
||||
|
||||
|
||||
class GaussNewton(Minimize, Remember):
|
||||
name = 'Gauss Newton'
|
||||
|
||||
@@ -643,16 +717,52 @@ class GaussNewton(Minimize, Remember):
|
||||
return Solver(self.H).solve(-self.g)
|
||||
|
||||
|
||||
class InexactGaussNewton(Minimize, Remember):
|
||||
class InexactGaussNewton(BFGS, Minimize, Remember):
|
||||
"""
|
||||
Minimizes using CG as the inexact solver of
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{H p = -g}
|
||||
|
||||
By default BFGS is used as the preconditioner.
|
||||
|
||||
Use *nbfgs* to set the memory limitation of BFGS.
|
||||
|
||||
To set the initial H0 to be used in BFGS, set *bfgsH0* to be a SimPEG.Solver
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
Minimize.__init__(self, **kwargs)
|
||||
|
||||
name = 'Inexact Gauss Newton'
|
||||
|
||||
maxIterCG = 10
|
||||
tolCG = 1e-5
|
||||
tolCG = 1e-3
|
||||
|
||||
@property
|
||||
def approxHinv(self):
|
||||
"""
|
||||
The approximate Hessian inverse is used to precondition CG.
|
||||
|
||||
Default uses BFGS, with an initial H0 of *bfgsH0*.
|
||||
|
||||
Must be a scipy.sparse.linalg.LinearOperator
|
||||
"""
|
||||
_approxHinv = getattr(self,'_approxHinv',None)
|
||||
if _approxHinv is None:
|
||||
M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype )
|
||||
return M
|
||||
return _approxHinv
|
||||
@approxHinv.setter
|
||||
def approxHinv(self, value):
|
||||
self._approxHinv = value
|
||||
|
||||
@timeIt
|
||||
def findSearchDirection(self):
|
||||
# TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly
|
||||
p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG)
|
||||
Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG})
|
||||
p = Hinv.solve(-self.g)
|
||||
return p
|
||||
|
||||
|
||||
@@ -663,6 +773,84 @@ class SteepestDescent(Minimize, Remember):
|
||||
def findSearchDirection(self):
|
||||
return -self.g
|
||||
|
||||
|
||||
class NewtonRoot(object):
|
||||
"""
|
||||
Newton Method - Root Finding
|
||||
|
||||
root = newtonRoot(fun,x);
|
||||
|
||||
Where fun is the function that returns the function value as well as the
|
||||
gradient.
|
||||
|
||||
For iterative solving of dh = -J\\r, use O.solveTol = TOL. For direct
|
||||
solves, use SOLVETOL = 0 (default)
|
||||
|
||||
Rowan Cockett
|
||||
16-May-2013 16:29:51
|
||||
University of British Columbia
|
||||
rcockett@eos.ubc.ca
|
||||
|
||||
"""
|
||||
|
||||
tol = 1.000e-06
|
||||
solveTol = 0 # Default direct solve.
|
||||
maxIter = 20
|
||||
stepDcr = 0.5
|
||||
maxLS = 30
|
||||
comments = False
|
||||
doLS = True
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
setKwargs(self, **kwargs)
|
||||
|
||||
def root(self, fun, x):
|
||||
if self.comments: print 'Newton Method:\n'
|
||||
|
||||
self._iter = 0
|
||||
while True:
|
||||
|
||||
[r,J] = fun(x);
|
||||
if self.solveTol == 0:
|
||||
Jinv = Solver(J)
|
||||
dh = - Jinv.solve(r)
|
||||
else:
|
||||
raise NotImplementedError('Iterative solve on NewtonRoot is not yet implemented.')
|
||||
# M = @(x) tril(J)\(diag(J).*(triu(J)\x));
|
||||
# [dh, ~] = bicgstab(J,-r,O.solveTol,500,M);
|
||||
|
||||
muLS = 1.
|
||||
LScnt = 1
|
||||
xt = x + dh
|
||||
rt, Jt = fun(xt) # TODO: get rid of Jt
|
||||
|
||||
if self.comments: print '\tLinesearch:\n'
|
||||
# Enter Linesearch
|
||||
while True and self.doLS:
|
||||
if self.comments:
|
||||
print '\t\tResid: %e\n'%norm(rt)
|
||||
if norm(rt) <= norm(r) or norm(rt) < self.tol:
|
||||
break
|
||||
|
||||
muLS = muLS*self.stepDcr
|
||||
LScnt = LScnt + 1
|
||||
print '.'
|
||||
if LScnt > self.maxLS:
|
||||
print 'Newton Method: Line search break.'
|
||||
root = NaN
|
||||
return
|
||||
xt = x + muLS*dh
|
||||
rt, Jt = fun(xt) # TODO: get rid of Jt
|
||||
|
||||
x = xt
|
||||
self._iter += 1
|
||||
if norm(rt) < self.tol or self._iter > self.maxIter:
|
||||
break
|
||||
|
||||
return x
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG.tests import Rosenbrock, checkDerivative
|
||||
import matplotlib.pyplot as plt
|
||||
@@ -677,3 +865,10 @@ if __name__ == '__main__':
|
||||
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)
|
||||
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
|
||||
|
||||
|
||||
print 'test the newtonRoot finding.'
|
||||
fun = lambda x: (np.sin(x), sdiag(np.cos(x)))
|
||||
x = np.array([np.pi-0.3, np.pi+0.1, 0])
|
||||
pnt = NewtonRoot(comments=False).root(fun,x)
|
||||
print pnt
|
||||
|
||||
@@ -7,7 +7,7 @@ class Regularization(object):
|
||||
@property
|
||||
def mref(self):
|
||||
if getattr(self, '_mref', None) is None:
|
||||
self._mref = np.zeros(self.mesh.nC);
|
||||
return np.zeros(self.mesh.nC);
|
||||
return self._mref
|
||||
@mref.setter
|
||||
def mref(self, value):
|
||||
@@ -105,8 +105,7 @@ class Regularization(object):
|
||||
|
||||
|
||||
@timeIt
|
||||
def modelObj2Deriv(self, m):
|
||||
mresid = m - self.mref
|
||||
def modelObj2Deriv(self):
|
||||
|
||||
mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
|
||||
|
||||
|
||||
@@ -241,7 +241,8 @@ function showClassDetail(cid, count) {
|
||||
for (var i = 0; i < count; i++) {
|
||||
tid = id_list[i];
|
||||
if (toHide) {
|
||||
document.getElementById('div_'+tid).style.display = 'none'
|
||||
var divTid = document.getElementById('div_'+tid);
|
||||
if(divTid !== null){divTid.style.display = 'none';}
|
||||
document.getElementById(tid).className = 'hiddenRow';
|
||||
}
|
||||
else {
|
||||
|
||||
@@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True):
|
||||
out += (H,)
|
||||
return out if len(out) > 1 else out[0]
|
||||
|
||||
def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10):
|
||||
"""
|
||||
Basic derivative check
|
||||
|
||||
@@ -201,6 +201,9 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
:param int num: number of times to reduce step length, h
|
||||
:param bool plotIt: if you would like to plot
|
||||
:param numpy.array dx: step direction
|
||||
:param int expectedOrder: The order that you expect the derivative to yield.
|
||||
:param float tolerance: The tolerance on the expected order.
|
||||
:param float eps: What is zero?
|
||||
:rtype: bool
|
||||
:return: did you pass the test?!
|
||||
|
||||
@@ -243,9 +246,6 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
order1 = np.log10(E1[:-1]/E1[1:])
|
||||
print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
|
||||
|
||||
tolerance = 0.9
|
||||
expectedOrder = 2
|
||||
eps = 1e-10
|
||||
order0 = order0[E0[1:] > eps]
|
||||
order1 = order1[E1[1:] > eps]
|
||||
belowTol = order1.size == 0 and order0.size > 0
|
||||
@@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
|
||||
|
||||
|
||||
|
||||
def getQuadratic(A, b):
|
||||
def getQuadratic(A, b, c=0):
|
||||
"""
|
||||
Given A and b, this returns a quadratic, Q
|
||||
Given A, b and c, this returns a quadratic, Q
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{Q( x ) = 0.5 x A x + b x}
|
||||
\mathbf{Q( x ) = 0.5 x A x + b x} + c
|
||||
"""
|
||||
def Quadratic(x, return_g=True, return_H=True):
|
||||
f = 0.5 * x.dot( A.dot(x)) + b.dot( x )
|
||||
f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c
|
||||
out = (f,)
|
||||
if return_g:
|
||||
g = A.dot(x) + b
|
||||
|
||||
@@ -2,7 +2,7 @@ import numpy as np
|
||||
import unittest
|
||||
from SimPEG.mesh import TensorMesh
|
||||
from SimPEG.utils import ModelBuilder, sdiag
|
||||
from SimPEG.forward import Problem, SyntheticProblem
|
||||
from SimPEG.forward import Problem
|
||||
from SimPEG.forward.DCProblem import *
|
||||
from TestUtils import checkDerivative
|
||||
from scipy.sparse.linalg import dsolve
|
||||
@@ -40,18 +40,13 @@ class DCProblemTests(unittest.TestCase):
|
||||
P = Q.T
|
||||
|
||||
# Create some data
|
||||
class syntheticDCProblem(DCProblem, SyntheticProblem):
|
||||
pass
|
||||
|
||||
synthetic = syntheticDCProblem(mesh);
|
||||
synthetic.P = P
|
||||
synthetic.RHS = q
|
||||
dobs, Wd = synthetic.createData(mSynth, std=0.05)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
problem = DCProblem(mesh)
|
||||
problem.P = P
|
||||
problem.RHS = q
|
||||
dobs, Wd = problem.createSyntheticData(mSynth, std=0.05)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
problem.W = Wd
|
||||
problem.dobs = dobs
|
||||
problem.std = dobs*0 + 0.05
|
||||
|
||||
+69
-12
@@ -1,9 +1,10 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sparse
|
||||
import scipy.sparse as sp
|
||||
import scipy.sparse.linalg as linalg
|
||||
from SimPEG.utils import mkvc
|
||||
from SimPEG.utils import mkvc, sdiag
|
||||
import warnings
|
||||
|
||||
DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
|
||||
DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
|
||||
|
||||
try:
|
||||
import TriSolve
|
||||
@@ -45,13 +46,44 @@ class Solver(object):
|
||||
def __init__(self, A, doDirect=True, flag=None, options={}):
|
||||
assert type(doDirect) is bool, 'doDirect must be a boolean'
|
||||
assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
|
||||
|
||||
assert type(options) is dict, 'options must be a dictionary object'
|
||||
self.A = A
|
||||
|
||||
self.dsolve = None
|
||||
self.doDirect = doDirect
|
||||
self.flag = flag
|
||||
self.options = options
|
||||
if doDirect: return
|
||||
|
||||
# Now deal with iterative stuff only
|
||||
if 'M' not in options:
|
||||
warnings.warn("You should provide a preconditioner, M.", UserWarning)
|
||||
return
|
||||
M = options['M']
|
||||
if type(M) is sp.linalg.LinearOperator:
|
||||
return
|
||||
PreconditionerList = ['J','GS']
|
||||
if type(M) is str:
|
||||
assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']"
|
||||
M = (M,A) # use A as the base for the preconditioner.
|
||||
if type(M) is tuple:
|
||||
assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix."
|
||||
if M[0] is 'J':
|
||||
Jacobi = sdiag(1.0/M[1].diagonal())
|
||||
options['M'] = Jacobi
|
||||
elif M[0] is 'GS':
|
||||
LL = sp.tril(M[1])
|
||||
UU = sp.triu(M[1])
|
||||
DD = sdiag(M[1].diagonal())
|
||||
Uinv = Solver(UU, flag='U')
|
||||
Linv = Solver(LL, flag='L')
|
||||
def GS(f):
|
||||
return Uinv.solve(DD*Linv.solve(f))
|
||||
options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype )
|
||||
|
||||
else:
|
||||
raise Exception('M must be a LinearOperator or a tuple')
|
||||
|
||||
|
||||
def solve(self, b):
|
||||
"""
|
||||
@@ -118,8 +150,20 @@ class Solver(object):
|
||||
|
||||
return X
|
||||
|
||||
def solveIter(self, b, M=None, iterSolver='CG'):
|
||||
pass
|
||||
def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
|
||||
if backend is None: backend = DEFAULTS['iter']
|
||||
|
||||
algorithms = {'CG':sp.linalg.cg}
|
||||
assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
|
||||
alg = algorithms[iterSolver]
|
||||
|
||||
if len(b.shape) == 1 or b.shape[1] == 1:
|
||||
x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter)
|
||||
else:
|
||||
x = np.empty_like(b)
|
||||
for i in range(b.shape[1]):
|
||||
x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter)
|
||||
return x
|
||||
|
||||
def solveBackward(self, b, backend=None):
|
||||
"""
|
||||
@@ -132,9 +176,8 @@ class Solver(object):
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['backward']
|
||||
if type(self.A) is not sparse.csr.csr_matrix:
|
||||
from scipy.sparse import csr_matrix
|
||||
self.A = csr_matrix(self.A)
|
||||
if type(self.A) is not sp.csr.csr_matrix:
|
||||
self.A = sp.csr_matrix(self.A)
|
||||
vals = self.A.data
|
||||
rowptr = self.A.indptr
|
||||
colind = self.A.indices
|
||||
@@ -164,7 +207,7 @@ class Solver(object):
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['forward']
|
||||
if type(self.A) is not sparse.csr.csr_matrix:
|
||||
if type(self.A) is not sp.csr.csr_matrix:
|
||||
from scipy.sparse import csr_matrix
|
||||
self.A = csr_matrix(self.A)
|
||||
vals = self.A.data
|
||||
@@ -240,13 +283,13 @@ if __name__ == '__main__':
|
||||
print np.linalg.norm(e-x,np.inf)
|
||||
|
||||
|
||||
n = 6000
|
||||
n = 600
|
||||
A_dense = np.random.random((n,n))
|
||||
L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned.
|
||||
e = np.ones(n)
|
||||
b = np.dot(L, e)
|
||||
|
||||
A = sparse.csr_matrix(L)
|
||||
A = sp.csr_matrix(L)
|
||||
pSolve = Solver(A,flag='L',options={'backend':'python'});
|
||||
fSolve = Solver(A,flag='L',options={'backend':'fortran'})
|
||||
tic = time()
|
||||
@@ -257,3 +300,17 @@ if __name__ == '__main__':
|
||||
x = fSolve.solve(b)
|
||||
toc = time() - tic
|
||||
print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc
|
||||
|
||||
|
||||
|
||||
A = -D*D.T
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
e = np.ones(M.nC)
|
||||
b = A.dot(e)
|
||||
|
||||
iSolve = Solver(A, doDirect=False,options={'M':('GS',A)})
|
||||
tic = time()
|
||||
x = iSolve.solve(b)
|
||||
toc = time() - tic
|
||||
print x
|
||||
print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info
|
||||
|
||||
@@ -26,3 +26,4 @@ Linear Problem
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
|
||||
@@ -1,9 +0,0 @@
|
||||
.. _api_Solver:
|
||||
|
||||
Solver
|
||||
******
|
||||
|
||||
.. automodule:: SimPEG.utils.Solver
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
+1024
-657
File diff suppressed because it is too large
Load Diff
@@ -1,24 +1,52 @@
|
||||
.. _api_Utils:
|
||||
|
||||
|
||||
Solver
|
||||
******
|
||||
|
||||
.. automodule:: SimPEG.utils.Solver
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
|
||||
Utilities
|
||||
*********
|
||||
|
||||
.. automodule:: SimPEG.utils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Matrix Utilities
|
||||
****************
|
||||
|
||||
.. automodule:: SimPEG.utils.matutils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Sparse Utilities
|
||||
****************
|
||||
|
||||
.. automodule:: SimPEG.utils.sputils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
LOM Utilities
|
||||
*************
|
||||
|
||||
.. automodule:: SimPEG.utils.lomutils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Model Builder Utilities
|
||||
***********************
|
||||
|
||||
.. automodule:: SimPEG.utils.ModelBuilder
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Interpolation Utilities
|
||||
***********************
|
||||
|
||||
.. automodule:: SimPEG.utils.interputils
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
@@ -56,7 +56,6 @@ Utility Codes
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
|
||||
api_Solver
|
||||
api_Utils
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,158 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": ""
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import SimPEG\n",
|
||||
"from SimPEG import Solver\n",
|
||||
"from SimPEG.mesh import TensorMesh\n",
|
||||
"from SimPEG.regularization import Regularization\n",
|
||||
"import SimPEG.inverse as inverse\n",
|
||||
"from SimPEG.inverse import Minimize, Remember, IterationPrinters\n",
|
||||
"import numpy as np\n",
|
||||
"import scipy.sparse as sp"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 2
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"FUN = SimPEG.tests.Rosenbrock\n",
|
||||
"FUN = SimPEG.tests.getQuadratic(sp.csr_matrix(([100,1],([0,1],[0,1])),shape=(2,2)),np.array([-5,-5]),100)\n",
|
||||
"\n",
|
||||
"x0 = np.array([1,0])\n",
|
||||
"opt = inverse.BFGS()\n",
|
||||
"xopt = opt.minimize(FUN,x0)\n",
|
||||
"print xopt\n",
|
||||
"opt = inverse.GaussNewton()\n",
|
||||
"xopt = opt.minimize(FUN,x0)\n",
|
||||
"print xopt\n",
|
||||
"opt = inverse.SteepestDescent()\n",
|
||||
"xopt = opt.minimize(FUN,x0)\n",
|
||||
"print xopt"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"===================== BFGS =====================\n",
|
||||
" # f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------\n",
|
||||
" 0 1.45e+02 9.51e+01 0 \n",
|
||||
" 1 1.14e+02 5.37e+01 6 \n",
|
||||
" 2 1.04e+02 3.04e+01 6 \n",
|
||||
" 3 8.83e+01 1.37e+01 0 \n",
|
||||
" 4 8.76e+01 5.97e+00 0 Skip BFGS \n",
|
||||
" 5 8.74e+01 2.61e+00 0 Skip BFGS \n",
|
||||
" 6 8.74e+01 1.14e+00 0 Skip BFGS \n",
|
||||
" 7 8.74e+01 5.01e-01 0 Skip BFGS \n",
|
||||
" 8 8.74e+01 2.19e-01 0 Skip BFGS \n",
|
||||
" 9 8.74e+01 9.60e-02 0 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 1.9437e-04 <= tolF*(1+|f0|) = 1.4600e+01\n",
|
||||
"1 : |xc-x_last| = 1.2663e-03 <= tolX*(1+|x0|) = 2.0000e-01\n",
|
||||
"1 : |proj(x-g)-x| = 9.5952e-02 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 9.5952e-02 <= 1e3*eps = 1.0000e-02\n",
|
||||
"0 : maxIter = 20 <= iter = 9\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"[ 0.05095952 4.99977449]\n",
|
||||
"=========== Gauss Newton ===========\n",
|
||||
" # f |proj(x-g)-x| LS \n",
|
||||
"-----------------------------------\n",
|
||||
" 0 1.45e+02 9.51e+01 0 \n",
|
||||
" 1 8.74e+01 4.44e-15 0 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"0 : |fc-fOld| = 5.7625e+01 <= tolF*(1+|f0|) = 1.4600e+01\n",
|
||||
"0 : |xc-x_last| = 5.0894e+00 <= tolX*(1+|x0|) = 2.0000e-01\n",
|
||||
"1 : |proj(x-g)-x| = 4.4409e-15 <= tolG = 1.0000e-01\n",
|
||||
"1 : |proj(x-g)-x| = 4.4409e-15 <= 1e3*eps = 1.0000e-02\n",
|
||||
"0 : maxIter = 20 <= iter = 1\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"[ 0.05 5. ]\n",
|
||||
"========= Steepest Descent =========\n",
|
||||
" # f |proj(x-g)-x| LS \n",
|
||||
"-----------------------------------\n",
|
||||
" 0 1.45e+02 9.51e+01 0 \n",
|
||||
" 1 1.14e+02 5.37e+01 6 \n",
|
||||
" 2 1.04e+02 3.04e+01 6 \n",
|
||||
" 3 1.00e+02 1.76e+01 6 \n",
|
||||
" 4 9.88e+01 1.06e+01 6 \n",
|
||||
" 5 9.82e+01 7.07e+00 6 \n",
|
||||
" 6 9.80e+01 1.22e+01 5 \n",
|
||||
" 7 9.73e+01 7.77e+00 6 \n",
|
||||
" 8 9.68e+01 5.64e+00 6 \n",
|
||||
" 9 9.65e+01 8.72e+00 5 \n",
|
||||
" 10 9.60e+01 5.97e+00 6 \n",
|
||||
" 11 9.58e+01 9.98e+00 5 \n",
|
||||
" 12 9.53e+01 6.48e+00 6 \n",
|
||||
" 13 9.53e+01 1.16e+01 5 \n",
|
||||
" 14 9.46e+01 7.20e+00 6 \n",
|
||||
" 15 9.43e+01 5.07e+00 6 \n",
|
||||
" 16 9.41e+01 8.17e+00 5 \n",
|
||||
" 17 9.37e+01 5.43e+00 6 \n",
|
||||
" 18 9.36e+01 9.42e+00 5 \n",
|
||||
" 19 9.32e+01 5.98e+00 6 \n",
|
||||
" 20 9.29e+01 4.32e+00 6 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 2.5913e-01 <= tolF*(1+|f0|) = 1.4600e+01\n",
|
||||
"1 : |xc-x_last| = 9.3379e-02 <= tolX*(1+|x0|) = 2.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 4.3246e+00 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 4.3246e+00 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 20 <= iter = 20\n",
|
||||
"------------------------- DONE! -------------------------\n",
|
||||
"[ 0.07777107 1.6849632 ]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stderr",
|
||||
"text": [
|
||||
"/Users/rowan/git/simpeg/SimPEG/inverse/Optimize.py:664: RuntimeWarning: divide by zero encountered in remainder\n",
|
||||
" khat = np.mod(n-nn+k,nn)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"A = sp.identity(2)\n",
|
||||
"S = Solver(A)\n",
|
||||
"\n",
|
||||
"assert type(S) is Solver"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": []
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
Reference in New Issue
Block a user