From a953a52ccc528726ea03ae33c514ce0fec5d1c69 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 31 May 2015 10:43:23 -0700 Subject: [PATCH 1/9] initial commit of PropMap --- SimPEG/Models.py | 4 +- SimPEG/Utils/PropMaps.py | 177 ++++++++++++++++++++++++++++++++++++++ SimPEG/Utils/codeutils.py | 4 +- 3 files changed, 182 insertions(+), 3 deletions(-) create mode 100644 SimPEG/Utils/PropMaps.py diff --git a/SimPEG/Models.py b/SimPEG/Models.py index f69b749c..65cfdeee 100644 --- a/SimPEG/Models.py +++ b/SimPEG/Models.py @@ -31,5 +31,5 @@ class Model(np.ndarray): @property def transformDeriv(self): if getattr(self, '_transformDeriv', None) is None: - self.deriv = self.mapping.deriv(self.view(np.ndarray)) - return self.deriv + self._transformDeriv = self.mapping.deriv(self.view(np.ndarray)) + return self._transformDeriv diff --git a/SimPEG/Utils/PropMaps.py b/SimPEG/Utils/PropMaps.py new file mode 100644 index 00000000..e4c49407 --- /dev/null +++ b/SimPEG/Utils/PropMaps.py @@ -0,0 +1,177 @@ +from SimPEG import Utils, Maps + +class Property(object): + + name = '' + doc = '' + # mappingPair = None + + defaultVal = None + defaultMap = None + defaultInvProp = False + + def __init__(self, doc, **kwargs): + # Set the default after all other params are set + self.doc = doc + Utils.setKwargs(self, **kwargs) + + def _getMapProperty(self): + prop = self + def fget(self): + return getattr(self, '_%sMap'%prop.name, prop.defaultMap) + def fset(self, val): + # TODO: Check if the mapping can be correct + setattr(self, '_%sMap'%prop.name, val) + return property(fget=fget, fset=fset, doc=prop.doc) + + def _getIndexProperty(self): + prop = self + def fget(self): + return getattr(self, '_%sIndex'%prop.name, slice(None)) + def fset(self, val): + setattr(self, '_%sIndex'%prop.name, val) + return property(fget=fget, fset=fset, doc=prop.doc) + + def _getProperty(self): + prop = self + def fget(self): + mapping = getattr(self, '%sMap'%prop.name) + if mapping is None: + return prop.defaultVal + m = getattr(self, '%sModel'%prop.name) + return mapping * m + return property(fget=fget) + + def _getModelDerivProperty(self): + prop = self + def fget(self): + m = getattr(self, '%sModel'%prop.name) + mapping = getattr(self, '%sMap'%prop.name) + if mapping is None: + return None + return mapping.deriv( m ) + return property(fget=fget) + + def _getModelProperty(self): + prop = self + def fget(self): + mapping = getattr(self, '%sMap'%prop.name) + if mapping is None: + return None + index = getattr(self.propMap, '_%sIndex'%prop.name, slice(None)) + return self.vector[index] + return property(fget=fget) + + def _getModelMapProperty(self): + prop = self + def fget(self): + return getattr(self.propMap, '_%sMap'%prop.name, prop.defaultMap) + return property(fget=fget) + + + +class PropModel(object): + def __init__(self, propMap, vector): + self.propMap = propMap + self.vector = vector + + # TODO: nP + + +_PROPMAPCLASSREGISTRY = {} + +class _PropMapMetaClass(type): + def __new__(cls, name, bases, attrs): + assert name.endswith('PropMap'), 'Please use convention: ___PropMap, e.g. ElectromagneticPropMap' + _properties = {} + for base in bases: + for baseProp in getattr(base, '_properties', {}): + _properties[baseProp] = base._properties[baseProp] + keys = [key for key in attrs] + for attr in keys: + if isinstance(attrs[attr], Property): + attrs[attr].name = attr + attrs[attr + 'Map' ] = attrs[attr]._getMapProperty() + attrs[attr + 'Index'] = attrs[attr]._getIndexProperty() + _properties[attr] = attrs[attr] + attrs.pop(attr) + + attrs['_properties'] = _properties + + newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs) + + newClass.PropModel = cls.createPropModelClass(newClass, name, _properties) + + _PROPMAPCLASSREGISTRY[name] = newClass + return newClass + + def createPropModelClass(self, name, _properties): + + attrs = dict() + + for attr in _properties: + prop = _properties[attr] + + attrs[attr ] = prop._getProperty() + attrs[attr + 'Map' ] = prop._getModelMapProperty() + attrs[attr + 'Model'] = prop._getModelProperty() + attrs[attr + 'Deriv'] = prop._getModelDerivProperty() + + return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs) + + + + + +class PropMap(object): + __metaclass__ = _PropMapMetaClass + + def __init__(self, ): + """ + PropMap takes a multi parameter model and maps it to the equivalent PropModel + """ + pass + + def __call__(self, vec): + return self.PropModel(self, vec) + + +class MyPropMap(PropMap): + sigma = Property("Electrical Conductivity", defaultInvProp=True) + # rho = InveseProperty(sigma) + +class My2PropMap(MyPropMap): + mu = Property("Electrical Conductivity", defaultVal=4e-10) + + +if __name__ == '__main__': + + + from SimPEG import Mesh + import numpy as np + + expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) + print expMap.nP + + # propMap = MyPropMap([('sigma', expMap), ('mu', IMap)], indices={'sigma':[1,2,3,4,7,8]}) + propMap = MyPropMap() + print [n for n in dir(propMap) if n[0] is not '_'] + propMap = My2PropMap() + print [n for n in dir(propMap) if n[0] is not '_'] + + propMap.sigmaMap = expMap + + print propMap.sigmaMap + print propMap.sigmaIndex + + mod = propMap(np.r_[1,2,3]) + print [n for n in dir(mod) if n[0] is not '_'] + print mod.sigmaModel + print mod.sigma + print mod.sigmaMap + print mod.sigmaDeriv + + print mod.mu + print mod.muMap + print mod.muModel + print mod.muDeriv diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 1399c1b5..0ba57b2e 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -58,9 +58,11 @@ def hook(obj, method, name=None, overwrite=False, silent=False): print 'Method '+name+' was not overwritten.' -def setKwargs(obj, **kwargs): +def setKwargs(obj, ignore=[], **kwargs): """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: + if attr in ignore: + continue if hasattr(obj, attr): setattr(obj, attr, kwargs[attr]) else: From fbda6ab53b3d5c47c4f0e507f2c64ea0745484db Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 09:34:16 -0700 Subject: [PATCH 2/9] updates to init of propMap --- SimPEG/Utils/PropMaps.py | 76 +++++++++++++++++++++++++++++++++++----- 1 file changed, 68 insertions(+), 8 deletions(-) diff --git a/SimPEG/Utils/PropMaps.py b/SimPEG/Utils/PropMaps.py index e4c49407..7b40a0ec 100644 --- a/SimPEG/Utils/PropMaps.py +++ b/SimPEG/Utils/PropMaps.py @@ -98,6 +98,14 @@ class _PropMapMetaClass(type): attrs['_properties'] = _properties + # You are only allowed one default inversion property. + defaultInvProps = [] + for p in _properties: + if _properties[p].defaultInvProp: + defaultInvProps += [p] + if len(defaultInvProps) > 1: + raise Exception('You have more than one default inversion property: %s' % defaultInvProps) + newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs) newClass.PropModel = cls.createPropModelClass(newClass, name, _properties) @@ -126,11 +134,65 @@ class _PropMapMetaClass(type): class PropMap(object): __metaclass__ = _PropMapMetaClass - def __init__(self, ): + def __init__(self, mappings): """ PropMap takes a multi parameter model and maps it to the equivalent PropModel """ - pass + if type(mappings) is dict: + assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"' + self.setup(mappings['maps'], slices=mappings['slices']) + if type(mappings) is list: + self.setup(mappings) + elif isinstance(mappings, Maps.IdentityMap): + self.setup([(self.defaultInvProp, mappings)]) + else: + raise Exception('mappings must be a dict, a mapping, or a list of tuples.') + + + def setup(self, maps, slices=None): + """ + Sets up the maps and slices for the PropertyMap + + + :param list maps: [('sigma', sigmaMap), ('mu', muMap), ...] + :param list slices: [('sigma', slice(0,nP)), ('mu', [1,2,5,6]), ...] + + """ + assert np.all([ + type(m) is tuple and + len(m)==2 and + type(m[0]) is str and + m[0] in self._properties and + isinstance(m[1], Maps.IdentityMap) + for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties])) + if slices is None: + slices = dict() + else: + assert np.all([ + s in self._properties and + (type(s) in [slice, list] or isinstance(s, np.ndarray)) + for s in slices]), 'Slices must be for each property' + + self.clearMaps() + + nP = 0 + for name, mapping in maps: + setattr(self, '%sMap'%name, mapping) + setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP))) + nP += mapping.nP + + + @property + def defaultInvProp(self): + for name in self._properties: + p = self._properties[name] + if p.defaultInvProp: + return p.name + + def clearMaps(self): + for name in self._properties: + setattr(self, '%sMap'%name, None) + setattr(self, '%sIndex'%name, None) def __call__(self, vec): return self.PropModel(self, vec) @@ -138,9 +200,6 @@ class PropMap(object): class MyPropMap(PropMap): sigma = Property("Electrical Conductivity", defaultInvProp=True) - # rho = InveseProperty(sigma) - -class My2PropMap(MyPropMap): mu = Property("Electrical Conductivity", defaultVal=4e-10) @@ -154,13 +213,14 @@ if __name__ == '__main__': print expMap.nP # propMap = MyPropMap([('sigma', expMap), ('mu', IMap)], indices={'sigma':[1,2,3,4,7,8]}) - propMap = MyPropMap() - print [n for n in dir(propMap) if n[0] is not '_'] - propMap = My2PropMap() + propMap = MyPropMap([('sigma',expMap)]) print [n for n in dir(propMap) if n[0] is not '_'] + # propMap = My2PropMap() + # print [n for n in dir(propMap) if n[0] is not '_'] propMap.sigmaMap = expMap + print propMap.defaultInvProp print propMap.sigmaMap print propMap.sigmaIndex From 8475eadccebba8d95e18fbee1e8b688b86136619 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 09:55:14 -0700 Subject: [PATCH 3/9] updates to propMap location and testing --- SimPEG/Maps.py | 1 + SimPEG/{Utils => }/PropMaps.py | 53 +++---------------------------- SimPEG/Tests/test_PropMaps.py | 57 ++++++++++++++++++++++++++++++++++ 3 files changed, 63 insertions(+), 48 deletions(-) rename SimPEG/{Utils => }/PropMaps.py (82%) create mode 100644 SimPEG/Tests/test_PropMaps.py diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 0b47de66..4ab79e64 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -1,5 +1,6 @@ import Utils, numpy as np, scipy.sparse as sp from Tests import checkDerivative +from PropMaps import PropMap, Property class IdentityMap(object): diff --git a/SimPEG/Utils/PropMaps.py b/SimPEG/PropMaps.py similarity index 82% rename from SimPEG/Utils/PropMaps.py rename to SimPEG/PropMaps.py index 7b40a0ec..7c22e149 100644 --- a/SimPEG/Utils/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -1,4 +1,4 @@ -from SimPEG import Utils, Maps +import Utils, Maps, numpy as np class Property(object): @@ -7,7 +7,6 @@ class Property(object): # mappingPair = None defaultVal = None - defaultMap = None defaultInvProp = False def __init__(self, doc, **kwargs): @@ -18,7 +17,7 @@ class Property(object): def _getMapProperty(self): prop = self def fget(self): - return getattr(self, '_%sMap'%prop.name, prop.defaultMap) + return getattr(self, '_%sMap'%prop.name, None) def fset(self, val): # TODO: Check if the mapping can be correct setattr(self, '_%sMap'%prop.name, val) @@ -65,7 +64,7 @@ class Property(object): def _getModelMapProperty(self): prop = self def fget(self): - return getattr(self.propMap, '_%sMap'%prop.name, prop.defaultMap) + return getattr(self.propMap, '_%sMap'%prop.name, None) return property(fget=fget) @@ -128,9 +127,6 @@ class _PropMapMetaClass(type): return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs) - - - class PropMap(object): __metaclass__ = _PropMapMetaClass @@ -141,7 +137,7 @@ class PropMap(object): if type(mappings) is dict: assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"' self.setup(mappings['maps'], slices=mappings['slices']) - if type(mappings) is list: + elif type(mappings) is list: self.setup(mappings) elif isinstance(mappings, Maps.IdentityMap): self.setup([(self.defaultInvProp, mappings)]) @@ -170,7 +166,7 @@ class PropMap(object): else: assert np.all([ s in self._properties and - (type(s) in [slice, list] or isinstance(s, np.ndarray)) + (type(slices[s]) in [slice, list] or isinstance(slices[s], np.ndarray)) for s in slices]), 'Slices must be for each property' self.clearMaps() @@ -181,7 +177,6 @@ class PropMap(object): setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP))) nP += mapping.nP - @property def defaultInvProp(self): for name in self._properties: @@ -197,41 +192,3 @@ class PropMap(object): def __call__(self, vec): return self.PropModel(self, vec) - -class MyPropMap(PropMap): - sigma = Property("Electrical Conductivity", defaultInvProp=True) - mu = Property("Electrical Conductivity", defaultVal=4e-10) - - -if __name__ == '__main__': - - - from SimPEG import Mesh - import numpy as np - - expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) - print expMap.nP - - # propMap = MyPropMap([('sigma', expMap), ('mu', IMap)], indices={'sigma':[1,2,3,4,7,8]}) - propMap = MyPropMap([('sigma',expMap)]) - print [n for n in dir(propMap) if n[0] is not '_'] - # propMap = My2PropMap() - # print [n for n in dir(propMap) if n[0] is not '_'] - - propMap.sigmaMap = expMap - - print propMap.defaultInvProp - print propMap.sigmaMap - print propMap.sigmaIndex - - mod = propMap(np.r_[1,2,3]) - print [n for n in dir(mod) if n[0] is not '_'] - print mod.sigmaModel - print mod.sigma - print mod.sigmaMap - print mod.sigmaDeriv - - print mod.mu - print mod.muMap - print mod.muModel - print mod.muDeriv diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py new file mode 100644 index 00000000..4e7c8bb7 --- /dev/null +++ b/SimPEG/Tests/test_PropMaps.py @@ -0,0 +1,57 @@ +import unittest +from SimPEG import * + + + +class MyPropMap(Maps.PropMap): + sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True) + mu = Maps.Property("Electrical Conductivity", defaultVal=4e-10) + + +class TestPropMaps(unittest.TestCase): + + def setUp(self): + pass + + def test_setup(self): + expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) + assert expMap.nP == 3 + + PM1 = MyPropMap(expMap) + PM2 = MyPropMap([('sigma', expMap)]) + PM3 = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':slice(0,3)}}) + + for PM in [PM1,PM2,PM3]: + assert PM.defaultInvProp == 'sigma' + assert PM.sigmaMap is not None + assert PM.sigmaMap is expMap + assert PM.sigmaIndex == slice(0,3) + assert getattr(PM, 'sigma', None) is None + assert PM.muMap is None + assert PM.muIndex is None + + m = PM(np.r_[1,2,3]) + assert m.mu == 4e-10 + assert m.muModel is None + assert m.muMap is None + assert m.muDeriv is None + + assert np.all(m.sigmaModel == np.r_[1,2,3]) + assert m.sigmaMap is expMap + assert np.all(m.sigma == np.exp(np.r_[1,2,3])) + assert m.sigmaDeriv is not None + + def test_slices(self): + expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) + PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}}) + assert PM.sigmaIndex == [2,1,0] + m = PM(np.r_[1,2,3]) + assert np.all(m.sigmaModel == np.r_[3,2,1]) + assert np.all(m.sigma == np.exp(np.r_[3,2,1])) + + + +if __name__ == '__main__': + unittest.main() + + From 2c48a69fb2283b4a7e6b249f3a23538fc0889267 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 10:15:20 -0700 Subject: [PATCH 4/9] testing of compressed maps. --- SimPEG/PropMaps.py | 16 ++++++++++++++- SimPEG/Tests/test_PropMaps.py | 37 ++++++++++++++++++++++++++++++++--- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 7c22e149..2a8e7ffc 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -73,8 +73,22 @@ class PropModel(object): def __init__(self, propMap, vector): self.propMap = propMap self.vector = vector + assert len(self.vector) == self.nP + + @property + def nP(self): + inds = [] + if getattr(self, '_nP', None) is None: + for name in self.propMap._properties: + index = getattr(self.propMap, '%sIndex'%name, None) + if index is not None: + if type(index) is slice: + inds += range(*index.indices(len(self.vector))) + else: + inds += list(index) + self._nP = len(set(inds)) + return self._nP - # TODO: nP _PROPMAPCLASSREGISTRY = {} diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index 4e7c8bb7..b869502b 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -1,11 +1,11 @@ import unittest from SimPEG import * - +from scipy.constants import mu_0 class MyPropMap(Maps.PropMap): sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True) - mu = Maps.Property("Electrical Conductivity", defaultVal=4e-10) + mu = Maps.Property("Electrical Conductivity", defaultVal=mu_0) class TestPropMaps(unittest.TestCase): @@ -31,7 +31,7 @@ class TestPropMaps(unittest.TestCase): assert PM.muIndex is None m = PM(np.r_[1,2,3]) - assert m.mu == 4e-10 + assert m.mu == mu_0 assert m.muModel is None assert m.muMap is None assert m.muDeriv is None @@ -41,6 +41,8 @@ class TestPropMaps(unittest.TestCase): assert np.all(m.sigma == np.exp(np.r_[1,2,3])) assert m.sigmaDeriv is not None + assert m.nP == 3 + def test_slices(self): expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}}) @@ -49,7 +51,36 @@ class TestPropMaps(unittest.TestCase): assert np.all(m.sigmaModel == np.r_[3,2,1]) assert np.all(m.sigma == np.exp(np.r_[3,2,1])) + def test_multiMap(self): + m = Mesh.TensorMesh((3,)) + expMap = Maps.ExpMap(m) + iMap = Maps.IdentityMap(m) + PM = MyPropMap([('sigma', expMap), ('mu', iMap)]) + pm = PM(np.r_[1,2,3,4,5,6]) + + assert pm.nP == 6 + + assert np.all(pm.sigmaModel == [1,2,3]) + assert np.all(pm.sigma == np.exp([1,2,3])) + assert np.all(pm.muModel == [4,5,6]) + assert np.all(pm.mu == [4,5,6]) + + + def test_multiMapCompressed(self): + m = Mesh.TensorMesh((3,)) + expMap = Maps.ExpMap(m) + iMap = Maps.IdentityMap(m) + PM = MyPropMap({'maps':[('sigma', expMap), ('mu', iMap)],'slices':{'mu':[0,1,2]}}) + + pm = PM(np.r_[1,2,3]) + + assert pm.nP == 3 + + assert np.all(pm.sigmaModel == [1,2,3]) + assert np.all(pm.sigma == np.exp([1,2,3])) + assert np.all(pm.muModel == [1,2,3]) + assert np.all(pm.mu == [1,2,3]) if __name__ == '__main__': unittest.main() From 592d169f9dc36380370afe7887b328eaa8830be7 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 13:28:57 -0700 Subject: [PATCH 5/9] Changes in base problem for propmap --- SimPEG/Problem.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 75e9d4bf..607cff5b 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -18,15 +18,27 @@ class BaseProblem(object): Solver = Solver #: A SimPEG Solver class. solverOpts = {} #: Sovler options as a kwarg dict - mapping = None #: A SimPEG.Map instance. mesh = None #: A SimPEG.Mesh instance. + PropMap = None #: A SimPEG PropertyMap class. + + @property + def mapping(self): + "A SimPEG.Map instance or a property map is PropMap is not None" + return getattr(self, '_mapping', None) + @mapping.setter + def mapping(self, val): + if self.PropMap is None: + val._assertMatchesPair(self.mapPair) + self._mapping = val + else: + self._mapping = self.PropMap(val) + def __init__(self, mesh, mapping=None, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." self.mesh = mesh self.mapping = mapping or Maps.IdentityMap(mesh) - self.mapping._assertMatchesPair(self.mapPair) @property def survey(self): @@ -62,7 +74,10 @@ class BaseProblem(object): def curModel(self, value): if value is self.curModel: return # it is the same! - self._curModel = Models.Model(value, self.mapping) + if self.PropMap is not None: + self._curModel = self.mapping(value) + else: + self._curModel = Models.Model(value, self.mapping) for prop in self.deleteTheseOnModelUpdate: if hasattr(self, prop): delattr(self, prop) From 3f4f71bf3ca3f504795cfaf1f6be31b1a4df46e9 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 14:01:51 -0700 Subject: [PATCH 6/9] reciprocal map --- SimPEG/Maps.py | 20 ++++++++++++++++++++ SimPEG/Tests/test_maps.py | 4 ++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 4ab79e64..4ceb35bd 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -221,6 +221,26 @@ class ExpMap(IdentityMap): """ return Utils.sdiag(np.exp(Utils.mkvc(m))) +class ReciprocalMap(IdentityMap): + """ + Reciprocal mapping. For example, electrical resistivity and conductivity. + + .. math:: + + \\rho = \\frac{1}{\sigma} + + """ + def _transform(self, m): + return 1.0 / Utils.mkvc(m) + + def inverse(self, D): + return 1.0 / Utils.mkvc(m) + + def deriv(self, m): + # TODO: if this is a tensor, you might have a problem. + return Utils.sdiag( - Utils.mkvc(m)**(-2) ) + + class LogMap(IdentityMap): """ diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py index 9f8fef3c..437e51fa 100644 --- a/SimPEG/Tests/test_maps.py +++ b/SimPEG/Tests/test_maps.py @@ -6,8 +6,8 @@ from scipy.sparse.linalg import dsolve TOL = 1e-14 -MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"] -MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"] +MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "ReciprocalMap"] +MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "ReciprocalMap"] class MapTests(unittest.TestCase): From 760f24ea33b296816a24ece5cc0680df73454026 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 14:29:58 -0700 Subject: [PATCH 7/9] mesh independent maps --- SimPEG/Maps.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 4ceb35bd..74a92efd 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -23,6 +23,8 @@ class IdentityMap(object): :rtype: int :return: number of parameters in the model """ + if self.mesh is None: + return '*' return self.mesh.nC @property @@ -33,8 +35,11 @@ class IdentityMap(object): :rtype: (int,int) :return: shape of the operator as a tuple """ + if self.mesh is None: + return ('*', self.nP) return (self.mesh.nC, self.nP) + def _transform(self, m): """ Changes the model into the physical property. @@ -99,17 +104,17 @@ class IdentityMap(object): def __mul__(self, val): if isinstance(val, IdentityMap): - if not self.shape[1] == val.shape[0]: + if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]: raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val))) return ComboMap([self, val]) elif isinstance(val, np.ndarray): - if not self.shape[1] == val.shape[0]: + if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]: raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape))) return self._transform(val) raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!') def __str__(self): - return "%s(%d,%d)" % (self.__class__.__name__, self.shape[0], self.shape[1]) + return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) class ComboMap(IdentityMap): """Combination of various maps.""" @@ -120,10 +125,10 @@ class ComboMap(IdentityMap): self.maps = [] for ii, m in enumerate(maps): assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!' - if ii > 0 and not self.shape[1] == m.shape[0]: + if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]: prev = self.maps[-1] - errArgs = (prev.__name__, prev.shape[0], prev.shape[1], m.__name__, m.shape[0], m.shape[1]) - raise ValueError('Dimension mismatch in map[%s] (%i, %i) and map[%s] (%i, %i).' % errArgs) + errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1]) + raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs) if isinstance(m, ComboMap): self.maps += m.maps @@ -156,7 +161,7 @@ class ComboMap(IdentityMap): return deriv def __str__(self): - return 'ComboMap[%s]%s' % (' * '.join([m.__str__() for m in self.maps]), str(self.shape)) + return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1]) class ExpMap(IdentityMap): From 4df00148a367f54869f010126f72a19237061fc2 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 15:06:58 -0700 Subject: [PATCH 8/9] property links --- SimPEG/PropMaps.py | 35 ++++++++++++++++++--- SimPEG/Tests/test_PropMaps.py | 58 ++++++++++++++++++++++++++++------- 2 files changed, 77 insertions(+), 16 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 2a8e7ffc..d4fcb85d 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -4,7 +4,6 @@ class Property(object): name = '' doc = '' - # mappingPair = None defaultVal = None defaultInvProp = False @@ -14,11 +13,23 @@ class Property(object): self.doc = doc Utils.setKwargs(self, **kwargs) + @property + def propertyLink(self): + "Can be something like: ('sigma', Maps.ReciprocalMap)" + return getattr(self, '_propertyLink', None) + @propertyLink.setter + def propertyLink(self, value): + assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name + self._propertyLink = value + def _getMapProperty(self): prop = self def fget(self): return getattr(self, '_%sMap'%prop.name, None) def fset(self, val): + if prop.propertyLink is not None: + linkName, linkMap = prop.propertyLink + assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.' # TODO: Check if the mapping can be correct setattr(self, '_%sMap'%prop.name, val) return property(fget=fget, fset=fset, doc=prop.doc) @@ -35,19 +46,33 @@ class Property(object): prop = self def fget(self): mapping = getattr(self, '%sMap'%prop.name) - if mapping is None: + if mapping is None and prop.propertyLink is None: return prop.defaultVal - m = getattr(self, '%sModel'%prop.name) + + if mapping is None and prop.propertyLink is not None: + linkName, linkMapClass = prop.propertyLink + linkMap = linkMapClass(None) + m = getattr(self, '%s'%linkName) + return linkMap * m + + m = getattr(self, '%sModel'%prop.name) return mapping * m return property(fget=fget) def _getModelDerivProperty(self): prop = self def fget(self): - m = getattr(self, '%sModel'%prop.name) mapping = getattr(self, '%sMap'%prop.name) - if mapping is None: + if mapping is None and prop.propertyLink is None: return None + + if mapping is None and prop.propertyLink is not None: + linkName, linkMapClass = prop.propertyLink + linkMap = linkMapClass(None) * getattr(self, '%sMap'%linkName) + m = getattr(self, '%s'%linkName) + return linkMap.deriv( m ) + + m = getattr(self, '%sModel'%prop.name) return mapping.deriv( m ) return property(fget=fget) diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index b869502b..bb5788bc 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -5,7 +5,12 @@ from scipy.constants import mu_0 class MyPropMap(Maps.PropMap): sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True) - mu = Maps.Property("Electrical Conductivity", defaultVal=mu_0) + mu = Maps.Property("Mu", defaultVal=mu_0) + +class MyReciprocalPropMap(Maps.PropMap): + sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True, propertyLink=('rho', Maps.ReciprocalMap)) + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + mu = Maps.Property("Mu", defaultVal=mu_0) class TestPropMaps(unittest.TestCase): @@ -30,15 +35,15 @@ class TestPropMaps(unittest.TestCase): assert PM.muMap is None assert PM.muIndex is None - m = PM(np.r_[1,2,3]) + m = PM(np.r_[1.,2,3]) assert m.mu == mu_0 assert m.muModel is None assert m.muMap is None assert m.muDeriv is None - assert np.all(m.sigmaModel == np.r_[1,2,3]) + assert np.all(m.sigmaModel == np.r_[1.,2,3]) assert m.sigmaMap is expMap - assert np.all(m.sigma == np.exp(np.r_[1,2,3])) + assert np.all(m.sigma == np.exp(np.r_[1.,2,3])) assert m.sigmaDeriv is not None assert m.nP == 3 @@ -47,7 +52,7 @@ class TestPropMaps(unittest.TestCase): expMap = Maps.ExpMap(Mesh.TensorMesh((3,))) PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}}) assert PM.sigmaIndex == [2,1,0] - m = PM(np.r_[1,2,3]) + m = PM(np.r_[1.,2,3]) assert np.all(m.sigmaModel == np.r_[3,2,1]) assert np.all(m.sigma == np.exp(np.r_[3,2,1])) @@ -57,14 +62,14 @@ class TestPropMaps(unittest.TestCase): iMap = Maps.IdentityMap(m) PM = MyPropMap([('sigma', expMap), ('mu', iMap)]) - pm = PM(np.r_[1,2,3,4,5,6]) + pm = PM(np.r_[1.,2,3,4,5,6]) assert pm.nP == 6 - assert np.all(pm.sigmaModel == [1,2,3]) - assert np.all(pm.sigma == np.exp([1,2,3])) - assert np.all(pm.muModel == [4,5,6]) - assert np.all(pm.mu == [4,5,6]) + assert np.all(pm.sigmaModel == [1.,2,3]) + assert np.all(pm.sigma == np.exp([1.,2,3])) + assert np.all(pm.muModel == [4.,5,6]) + assert np.all(pm.mu == [4.,5,6]) def test_multiMapCompressed(self): @@ -73,7 +78,7 @@ class TestPropMaps(unittest.TestCase): iMap = Maps.IdentityMap(m) PM = MyPropMap({'maps':[('sigma', expMap), ('mu', iMap)],'slices':{'mu':[0,1,2]}}) - pm = PM(np.r_[1,2,3]) + pm = PM(np.r_[1,2.,3]) assert pm.nP == 3 @@ -82,6 +87,37 @@ class TestPropMaps(unittest.TestCase): assert np.all(pm.muModel == [1,2,3]) assert np.all(pm.mu == [1,2,3]) + def test_Links(self): + m = Mesh.TensorMesh((3,)) + expMap = Maps.ExpMap(m) + iMap = Maps.IdentityMap(m) + PM = MyReciprocalPropMap([('sigma', iMap)]) + pm = PM(np.r_[1,2.,3]) + # print pm.sigma + # print pm.sigmaMap + assert np.all(pm.sigma == [1,2,3]) + assert np.all(pm.rho == 1./np.r_[1,2,3]) + assert pm.sigmaMap is iMap + assert pm.rhoMap is None + assert pm.sigmaDeriv is not None + assert pm.rhoDeriv is not None + + PM = MyReciprocalPropMap([('rho', iMap)]) + pm = PM(np.r_[1,2.,3]) + # print pm.sigma + # print pm.sigmaMap + assert np.all(pm.sigma == 1./np.r_[1,2,3]) + assert np.all(pm.rho == [1,2,3]) + assert pm.sigmaMap is None + assert pm.rhoMap is iMap + assert pm.sigmaDeriv is not None + assert pm.rhoDeriv is not None + + self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)]) + self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)]) + + MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine + if __name__ == '__main__': unittest.main() From 5b45fc628edc83f6e921063c31ee15a722f58f05 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 15:17:35 -0700 Subject: [PATCH 9/9] updates to testing class creation --- SimPEG/PropMaps.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index d4fcb85d..3bc11d3c 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -136,13 +136,19 @@ class _PropMapMetaClass(type): attrs['_properties'] = _properties - # You are only allowed one default inversion property. defaultInvProps = [] + for p in _properties: + prop = _properties[p] + if prop.defaultInvProp: + defaultInvProps += [p] + if prop.propertyLink is not None: + assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0]) + if len(defaultInvProps) > 1: + raise Exception('You have more than one default inversion property: %s' % defaultInvProps) + for p in _properties: if _properties[p].defaultInvProp: defaultInvProps += [p] - if len(defaultInvProps) > 1: - raise Exception('You have more than one default inversion property: %s' % defaultInvProps) newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)