Merge pull request #101 from simpeg/PropMap

Property Maps
This commit is contained in:
Rowan Cockett
2015-06-01 15:31:51 -07:00
7 changed files with 421 additions and 15 deletions
+33 -7
View File
@@ -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):
"""
+2 -2
View File
@@ -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
+18 -3
View File
@@ -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)
+239
View File
@@ -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)
+124
View File
@@ -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()
+2 -2
View File
@@ -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):
+3 -1
View File
@@ -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: