diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 0b47de66..74a92efd 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): @@ -22,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 @@ -32,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. @@ -98,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.""" @@ -119,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 @@ -155,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): @@ -220,6 +226,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/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/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) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py new file mode 100644 index 00000000..3bc11d3c --- /dev/null +++ b/SimPEG/PropMaps.py @@ -0,0 +1,239 @@ +import Utils, Maps, numpy as np + +class Property(object): + + name = '' + doc = '' + + defaultVal = None + defaultInvProp = False + + def __init__(self, doc, **kwargs): + # Set the default after all other params are set + 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) + + 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 and prop.propertyLink is None: + return prop.defaultVal + + 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): + mapping = getattr(self, '%sMap'%prop.name) + 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) + + 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, None) + return property(fget=fget) + + + +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 + + + +_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 + + 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] + + 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, mappings): + """ + PropMap takes a multi parameter model and maps it to the equivalent PropModel + """ + 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']) + elif 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(slices[s]) in [slice, list] or isinstance(slices[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) + diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py new file mode 100644 index 00000000..bb5788bc --- /dev/null +++ b/SimPEG/Tests/test_PropMaps.py @@ -0,0 +1,124 @@ +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("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): + + 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 == 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 m.sigmaMap is expMap + 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]}}) + 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])) + + 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]) + + 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() + + 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): 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: