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

Merge branch 'master' of https://bitbucket.org/rcockett/simpeg into richards

Browse Source 
Conflicts:
	SimPEG/inverse/Optimize.py
	SimPEG/regularization/Regularization.py
This commit is contained in:
Rowan Cockett
2013-11-21 11:56:12 -08:00
parent 95fd98c23e f87899465d
commit 82edbf7a60
6 changed files with 223 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/forward/Problem.py
+9 -1
View File
@@ -1,5 +1,5 @@
import numpy as np
from SimPEG.utils import mkvc, sdiag
from SimPEG.utils import mkvc, sdiag, count, timeIt
import scipy.sparse as sp
norm = np.linalg.norm
@@ -37,6 +37,8 @@ class Problem(object):
to (locally) find how model parameters change the data, and optimize!
"""
counter = None
def __init__(self, mesh):
self.mesh = mesh
@@ -83,6 +85,7 @@ class Problem(object):
def dobs(self, value):
self._dobs = value
@count
def dpred(self, m, u=None):
"""
Predicted data.
@@ -94,6 +97,7 @@ class Problem(object):
u = self.field(m)
return self.P*u
@count
def dataResidual(self, m, u=None):
"""
:param numpy.array m: geophysical model
@@ -113,6 +117,7 @@ class Problem(object):
return self.dpred(m, u=u) - self.dobs
@timeIt
def J(self, m, v, u=None):
"""
:param numpy.array m: model
@@ -142,6 +147,7 @@ class Problem(object):
"""
raise NotImplementedError('J is not yet implemented.')
@timeIt
def Jt(self, m, v, u=None):
"""
:param numpy.array m: model
@@ -155,6 +161,7 @@ class Problem(object):
raise NotImplementedError('Jt is not yet implemented.')
@timeIt
def J_approx(self, m, v, u=None):
"""
@@ -169,6 +176,7 @@ class Problem(object):
"""
return self.J(m, v, u)
@timeIt
def Jt_approx(self, m, v, u=None):
"""
:param numpy.array m: model
SimPEG/inverse/Inversion.py
+8 -3
View File
@@ -1,7 +1,7 @@
import numpy as np
import scipy.sparse as sp
import SimPEG
from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers
from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt
from Optimize import Remember
from BetaSchedule import Cooling
@@ -13,6 +13,8 @@ class BaseInversion(object):
debug = False
beta0 = 1e4
counter = None
def __init__(self, prob, reg, opt, **kwargs):
setKwargs(self, **kwargs)
self.prob = prob
@@ -64,6 +66,7 @@ class BaseInversion(object):
def phi_d_target(self, value):
self._phi_d_target = value
@timeIt
def run(self, m0):
self.startup(m0)
while True:
@@ -145,7 +148,7 @@ class BaseInversion(object):
"""
printStoppers(self, self.stoppers)
@timeIt
def evalFunction(self, m, return_g=True, return_H=True):
u = self.prob.field(m)
@@ -176,7 +179,7 @@ class BaseInversion(object):
out += (operator,)
return out if len(out) > 1 else out[0]
@timeIt
def dataObj(self, m, u=None):
"""
:param numpy.array m: geophysical model
@@ -198,6 +201,7 @@ class BaseInversion(object):
R = mkvc(R)
return 0.5*np.vdot(R, R)
@timeIt
def dataObjDeriv(self, m, u=None):
"""
:param numpy.array m: geophysical model
@@ -238,6 +242,7 @@ class BaseInversion(object):
return dmisfit
@timeIt
def dataObj2Deriv(self, m, v, u=None):
"""
:param numpy.array m: geophysical model
SimPEG/inverse/Optimize.py
+22 -3
View File
@@ -1,6 +1,6 @@
import numpy as np
import matplotlib.pyplot as plt
from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers
from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt
norm = np.linalg.norm
import scipy.sparse as sp
from SimPEG import Solver
@@ -107,6 +107,7 @@ class Minimize(object):
debugLS = False
comment = ''
counter = None
def __init__(self, **kwargs):
self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub
@@ -118,6 +119,7 @@ class Minimize(object):
setKwargs(self, **kwargs)
@timeIt
def minimize(self, evalFunction, x0):
"""
Minimizes the function (evalFunction) starting at the location x0.
@@ -265,6 +267,7 @@ class Minimize(object):
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.
@@ -296,14 +299,14 @@ class Minimize(object):
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
self.g0 = self.g
return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
@timeIt
def projection(self, p):
"""
projects the search direction.
@@ -316,6 +319,7 @@ class Minimize(object):
"""
return p
@timeIt
def findSearchDirection(self):
"""
**findSearchDirection** should return an approximation of:
@@ -345,6 +349,7 @@ class Minimize(object):
"""
return -self.g
@count
def scaleSearchDirection(self, p):
"""
**scaleSearchDirection** should scale the search direction if appropriate.
@@ -362,6 +367,7 @@ class Minimize(object):
nameLS = "Armijo linesearch"
@timeIt
def modifySearchDirection(self, p):
"""
**modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
@@ -398,6 +404,7 @@ class Minimize(object):
return self._LS_xt, self._iterLS < self.maxIterLS
@count
def modifySearchDirectionBreak(self, p):
"""
Code is called if modifySearchDirection fails
@@ -418,6 +425,7 @@ class Minimize(object):
print 'The linesearch got broken. Boo.'
return p, False
@count
def doEndIteration(self, xt):
"""
**doEndIteration** is called at the end of each minimize iteration.
@@ -536,18 +544,22 @@ class ProjectedGradient(Minimize, Remember):
self.aSet_prev = self.activeSet(x0)
@count
def projection(self, 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"""
return np.logical_or(x == self.lower, x == self.upper)
@count
def inactiveSet(self, x):
"""The free variables."""
return np.logical_not(self.activeSet(x))
@count
def bindingSet(self, x):
"""
If we are on a bound and the negative gradient points away from the feasible set.
@@ -559,6 +571,7 @@ class ProjectedGradient(Minimize, Remember):
bind_low = np.logical_and(x == self.upper, self.g <= 0)
return np.logical_or(bind_up, bind_low)
@timeIt
def findSearchDirection(self):
self.aSet_prev = self.activeSet(self.xc)
allBoundsAreActive = sum(self.aSet_prev) == self.xc.size
@@ -599,6 +612,7 @@ class ProjectedGradient(Minimize, Remember):
# aSet_after = self.activeSet(self.xc+p)
return p
@timeIt
def _doEndIteration_ProjectedGradient(self, xt):
aSet = self.activeSet(xt)
bSet = self.bindingSet(xt)
@@ -697,6 +711,8 @@ class BFGS(Minimize, Remember):
class GaussNewton(Minimize, Remember):
name = 'Gauss Newton'
@timeIt
def findSearchDirection(self):
return Solver(self.H).solve(-self.g)
@@ -743,6 +759,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
def approxHinv(self, value):
self._approxHinv = value
@timeIt
def findSearchDirection(self):
Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG})
p = Hinv.solve(-self.g)
@@ -751,6 +768,8 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
class SteepestDescent(Minimize, Remember):
name = 'Steepest Descent'
@timeIt
def findSearchDirection(self):
return -self.g
SimPEG/regularization/Regularization.py
+9 -8
View File
@@ -1,4 +1,4 @@
from SimPEG.utils import sdiag
from SimPEG.utils import sdiag, count, timeIt
import numpy as np
class Regularization(object):
@@ -40,21 +40,20 @@ class Regularization(object):
self._Wz = sdiag(a)*self.mesh.cellGradz
return self._Wz
alpha_s = 1e-6
alpha_x = 1
alpha_y = 1
alpha_z = 1
counter = None
def __init__(self, mesh):
self.mesh = mesh
self._Wx = None
self._Wy = None
self._Wz = None
self.alpha_s = 1e-6
self.alpha_x = 1
self.alpha_y = 1
self.alpha_z = 1
def pnorm(self, r):
return 0.5*r.dot(r)
@timeIt
def modelObj(self, m):
mresid = m - self.mref
@@ -69,6 +68,7 @@ class Regularization(object):
return mobj
@timeIt
def modelObjDeriv(self, m):
"""
@@ -104,6 +104,7 @@ class Regularization(object):
return mobjDeriv
@timeIt
def modelObj2Deriv(self):
mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
SimPEG/utils/ModelBuilder.py
+64
View File
@@ -1,4 +1,6 @@
import numpy as np
import scipy.ndimage as ndi
import scipy.sparse as sp
def getIndecesBlock(p0,p1,ccMesh):
@@ -130,6 +132,68 @@ def scalarConductivity(ccMesh,pFunction):
return sigma
def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
"""
Create a random model by convolving a kernal with a
uniformly distributed model.
:param int,tuple shape: shape of the model.
:param int seed: pick which model to produce, prints the seed if you don't choose.
:param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
:param int its: number of smoothing iterations
:param list bounds: bounds on the model, len(list) == 2
:rtype: numpy.ndarray
:return: M, the model
.. plot::
import matplotlib.pyplot as plt
import SimPEG.utils.ModelBuilder as MB
plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0])))
plt.title('A very cool, yet completely random model.')
plt.show()
"""
if seed is None:
seed = np.random.randint(1e3)
print 'Using a seed of: ', seed
if type(shape) in [int, long, float]:
shape = (shape,) # make it a tuple for consistency
np.random.seed(seed)
mr = np.random.rand(*shape)
if anisotropy is None:
if len(shape) is 1:
smth = np.array([1,10.,1],dtype=float)
elif len(shape) is 2:
smth = np.array([[1,2,1],[7,10,7],[1,2,1]],dtype=float)
elif len(shape) is 3:
kernal = np.array([1,4,1], dtype=float).reshape((1,3))
smth = np.array(sp.kron(sp.kron(kernal,kernal.T).todense()[:],kernal).todense()).reshape((3,3,3))
else:
assert len(anisotropy.shape) is len(shape), 'Anisotropy must be the same shape.'
smth = np.array(anisotropy,dtype=float)
smth = smth/smth.sum() # normalize
mi = mr
for i in range(its):
mi = ndi.convolve(mi, smth)
# scale the model to live between the bounds.
mi = (mi - mi.min())/(mi.max()-mi.min()) # scaled between 0 and 1
mi = mi*(bounds[1]-bounds[0])+bounds[0]
return mi
if __name__ == '__main__':
from SimPEG.mesh import TensorMesh
SimPEG/utils/__init__.py
+111
View File
@@ -60,3 +60,114 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
r = stopper['right'](obj)
print pad + stopper['str'] % (l<=r,l,r)
print pad + "%s%s%s" % ('-'*25,done,'-'*25)
import time
import numpy as np
class Counter(object):
"""
Counter allows anything that calls it to record iterations and
timings in a simple way.
Also has plotting functions that allow quick recalls of data.
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):
self.counter = Counter()
@count
def MyMethod(self):
pass
@timeIt
def MySecondMethod(self):
pass
c = MyClass('blah')
for i in range(100): c.MyMethod()
for i in range(300): c.MySecondMethod()
c.counter.summary()
"""
def __init__(self):
self._countList = {}
self._timeList = {}
def count(self, prop):
"""
Increases the count of the property.
"""
assert type(prop) is str, 'The property must be a string.'
if prop not in self._countList:
self._countList[prop] = 0
self._countList[prop] += 1
def countTic(self, prop):
"""
Times a property call, this is the init call.
"""
assert type(prop) is str, 'The property must be a string.'
if prop not in self._timeList:
self._timeList[prop] = []
self._timeList[prop].append(-time.time())
def countToc(self, prop):
"""
Times a property call, this is the end call.
"""
assert type(prop) is str, 'The property must be a string.'
assert prop in self._timeList, 'The property must already be in the dictionary.'
self._timeList[prop][-1] += time.time()
def summary(self):
"""
Provides a text summary of the current counters and timers.
"""
print 'Counters:'
for prop in sorted(self._countList):
print " {0:<40}: {1:8d}".format(prop,self._countList[prop])
print '\nTimes:'+' '*40+'mean sum'
for prop in sorted(self._timeList):
l = len(self._timeList[prop])
a = np.array(self._timeList[prop])
print " {0:<40}: {1:4.2e}, {2:4.2e}, {3:4d}x".format(prop,a.mean(),a.sum(),l)
def count(f):
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
out = f(self,*args,**kwargs)
return out
return wrapper
def timeIt(f):
def wrapper(self,*args,**kwargs):
counter = getattr(self,'counter',None)
if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
out = f(self,*args,**kwargs)
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):
self.counter = Counter()
@count
def MyMethod(self):
pass
@timeIt
def MySecondMethod(self):
pass
c = MyClass('blah')
for i in range(100): c.MyMethod()
for i in range(300): c.MySecondMethod()
c.counter.summary()