diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 35d41cb8..304f3912 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -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 diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index e0b62866..5a71f6ec 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -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 diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index c2f95dd9..75209dc6 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -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 diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 72335661..3e429324 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -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 diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index d0512680..56df6250 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -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 diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index bc278c7f..2ba0f44d 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -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()