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

fixes to the meta classes including soft linking inside the database. These are reflected in the loaded outputs which return the same object.

Browse Source
This commit is contained in:
rowanc1
2013-12-06 14:19:37 -08:00
parent 5e0fb8642d
commit 9e26543931
6 changed files with 157 additions and 117 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/inverse/Inversion.py
+1 -1
View File
@@ -187,7 +187,7 @@ class BaseInversion(object):
**printDone** is called at the end of the inversion routine.
"""
printStoppers(self, self.stoppers)
utils.printStoppers(self, self.stoppers)
@utils.callHooks('finish')
def finish(self):
SimPEG/mesh/BaseMesh.py
-1
View File
@@ -11,7 +11,6 @@ class BaseMesh(object):
:param numpy.array,list x0: Origin of the mesh (dim, )
"""
__metaclass__ = utils.Save.Savable
def __init__(self, n, x0=None):
SimPEG/mesh/LogicallyOrthogonalMesh.py
+49 -47
View File
@@ -1,15 +1,14 @@
import numpy as np
from SimPEG import utils, np
from BaseMesh import BaseMesh
from DiffOperators import DiffOperators
from InnerProducts import InnerProducts
from LomView import LomView
from SimPEG.utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
# Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), mkvc(length3D(x), 2))
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2))
class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
@@ -21,6 +20,9 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
.. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py
"""
__metaclass__ = utils.Save.Savable
_meshType = 'LOM'
def __init__(self, nodes):
@@ -38,7 +40,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
# Save nodes to private variable _gridN as vectors
self._gridN = np.ones((nodes[0].size, self.dim))
for i, node_i in enumerate(nodes):
self._gridN[:, i] = mkvc(node_i.astype(float))
self._gridN[:, i] = utils.mkvc(node_i.astype(float))
def gridCC():
doc = "Cell-centered grid."
@@ -69,10 +71,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridFx is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
self._gridFx = np.c_[XY[0], XY[1]]
elif self.dim == 3:
XYZ = [mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
XYZ = [utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFx
return locals()
@@ -86,10 +88,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridFy is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
self._gridFy = np.c_[XY[0], XY[1]]
elif self.dim == 3:
XYZ = [mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
XYZ = [utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFy
return locals()
@@ -102,7 +104,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if self._gridFz is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M')
XYZ = [mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
XYZ = [utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFz
return locals()
@@ -116,10 +118,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridEx is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
self._gridEx = np.c_[XY[0], XY[1]]
elif self.dim == 3:
XYZ = [mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
XYZ = [utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEx
return locals()
@@ -133,10 +135,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridEy is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
self._gridEy = np.c_[XY[0], XY[1]]
elif self.dim == 3:
XYZ = [mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
XYZ = [utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEy
return locals()
@@ -149,7 +151,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if self._gridEz is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M')
XYZ = [mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
XYZ = [utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEz
return locals()
@@ -192,25 +194,25 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if(self._vol is None):
if self.dim == 2:
A, B, C, D = indexCube('ABCD', self.n+1)
normal, area = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
A, B, C, D = utils.indexCube('ABCD', self.n+1)
normal, area = utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
self._vol = area
elif self.dim == 3:
# Each polyhedron can be decomposed into 5 tetrahedrons
# However, this presents a choice so we may as well divide in two ways and average.
A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', self.n+1)
A, B, C, D, E, F, G, H = utils.indexCube('ABCDEFGH', self.n+1)
vol1 = (volTetra(self.gridN, A, B, D, E) + # cutted edge top
volTetra(self.gridN, B, E, F, G) + # cutted edge top
volTetra(self.gridN, B, D, E, G) + # middle
volTetra(self.gridN, B, C, D, G) + # cutted edge bottom
volTetra(self.gridN, D, E, G, H)) # cutted edge bottom
vol1 = (utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
utils.volTetra(self.gridN, B, D, E, G) + # middle
utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom
utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom
vol2 = (volTetra(self.gridN, A, F, B, C) + # cutted edge top
volTetra(self.gridN, A, E, F, H) + # cutted edge top
volTetra(self.gridN, A, H, F, C) + # middle
volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
vol2 = (utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top
utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top
utils.volTetra(self.gridN, A, H, F, C) + # middle
utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
self._vol = (vol1 + vol2)/2
return self._vol
@@ -226,30 +228,30 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
# Compute areas of cell faces
if(self.dim == 2):
xy = self.gridN
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
edge1 = xy[B, :] - xy[A, :]
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
area1 = length2D(edge1)
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
# Note that we are doing A-D to make sure the normal points the right way.
# Think about it. Look at the picture. Normal points towards C iff you do this.
edge2 = xy[A, :] - xy[D, :]
normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
area2 = length2D(edge2)
self._area = np.r_[mkvc(area1), mkvc(area2)]
self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)]
self._normals = [normalize2D(normal1), normalize2D(normal2)]
elif(self.dim == 3):
A, E, F, B = indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
normal1, area1 = faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
A, E, F, B = utils.indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
normal1, area1 = utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
A, D, H, E = indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
normal2, area2 = faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
A, D, H, E = utils.indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
normal2, area2 = utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
A, B, C, D = indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
normal3, area3 = faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
A, B, C, D = utils.indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
normal3, area3 = utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)]
self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)]
self._normals = [normal1, normal2, normal3]
return self._area
return locals()
@@ -289,21 +291,21 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if(self._edge is None or self._tangents is None):
if(self.dim == 2):
xy = self.gridN
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
edge1 = xy[D, :] - xy[A, :]
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
edge2 = xy[B, :] - xy[A, :]
self._edge = np.r_[mkvc(length2D(edge1)), mkvc(length2D(edge2))]
self._edge = np.r_[utils.mkvc(length2D(edge1)), utils.mkvc(length2D(edge2))]
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
elif(self.dim == 3):
xyz = self.gridN
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
edge1 = xyz[D, :] - xyz[A, :]
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
edge2 = xyz[B, :] - xyz[A, :]
A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
A, E = utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
edge3 = xyz[E, :] - xyz[A, :]
self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), mkvc(length3D(edge3))]
self._edge = np.r_[utils.mkvc(length3D(edge1)), utils.mkvc(length3D(edge2)), utils.mkvc(length3D(edge3))]
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
return self._edge
return locals()
@@ -329,10 +331,10 @@ if __name__ == '__main__':
h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
dee3 = True
if dee3:
X, Y, Z = ndgrid(h1, h2, h3, vector=False)
X, Y, Z = utils.ndgrid(h1, h2, h3, vector=False)
M = LogicallyOrthogonalMesh([X, Y, Z])
else:
X, Y = ndgrid(h1, h2, vector=False)
X, Y = utils.ndgrid(h1, h2, vector=False)
M = LogicallyOrthogonalMesh([X, Y])
print M.r(M.normals, 'F', 'Fx', 'V')
SimPEG/mesh/TensorMesh.py
+30 -30
View File
@@ -1,11 +1,8 @@
import numpy as np
import scipy.sparse as sp
from SimPEG import utils, np, sp
from BaseMesh import BaseMesh
from TensorView import TensorView
from DiffOperators import DiffOperators
from InnerProducts import InnerProducts
from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat
class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
"""
@@ -35,6 +32,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
mesh = TensorMesh([10, 12, 15])
"""
__metaclass__ = utils.Save.Savable
_meshType = 'TENSOR'
def __init__(self, h_in, x0=None):
@@ -52,7 +52,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
# Ensure h contains 1D vectors
self._h = [mkvc(x.astype(float)) for x in h]
self._h = [utils.mkvc(x.astype(float)) for x in h]
def __str__(self):
outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim)
@@ -170,7 +170,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridCC is None:
self._gridCC = ndgrid(self.getTensor('CC'))
self._gridCC = utils.ndgrid(self.getTensor('CC'))
return self._gridCC
return locals()
_gridCC = None # Store grid by default
@@ -181,7 +181,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridN is None:
self._gridN = ndgrid(self.getTensor('N'))
self._gridN = utils.ndgrid(self.getTensor('N'))
return self._gridN
return locals()
_gridN = None # Store grid by default
@@ -192,7 +192,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFx is None:
self._gridFx = ndgrid(self.getTensor('Fx'))
self._gridFx = utils.ndgrid(self.getTensor('Fx'))
return self._gridFx
return locals()
_gridFx = None # Store grid by default
@@ -203,7 +203,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFy is None and self.dim > 1:
self._gridFy = ndgrid(self.getTensor('Fy'))
self._gridFy = utils.ndgrid(self.getTensor('Fy'))
return self._gridFy
return locals()
_gridFy = None # Store grid by default
@@ -214,7 +214,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFz is None and self.dim > 2:
self._gridFz = ndgrid(self.getTensor('Fz'))
self._gridFz = utils.ndgrid(self.getTensor('Fz'))
return self._gridFz
return locals()
_gridFz = None # Store grid by default
@@ -225,7 +225,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEx is None:
self._gridEx = ndgrid(self.getTensor('Ex'))
self._gridEx = utils.ndgrid(self.getTensor('Ex'))
return self._gridEx
return locals()
_gridEx = None # Store grid by default
@@ -236,7 +236,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEy is None and self.dim > 1:
self._gridEy = ndgrid(self.getTensor('Ey'))
self._gridEy = utils.ndgrid(self.getTensor('Ey'))
return self._gridEy
return locals()
_gridEy = None # Store grid by default
@@ -247,7 +247,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEz is None and self.dim > 2:
self._gridEz = ndgrid(self.getTensor('Ez'))
self._gridEz = utils.ndgrid(self.getTensor('Ez'))
return self._gridEz
return locals()
_gridEz = None # Store grid by default
@@ -262,13 +262,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
vh = self.h
# Compute cell volumes
if(self.dim == 1):
self._vol = mkvc(vh[0])
self._vol = utils.mkvc(vh[0])
elif(self.dim == 2):
# Cell sizes in each direction
self._vol = mkvc(np.outer(vh[0], vh[1]))
self._vol = utils.mkvc(np.outer(vh[0], vh[1]))
elif(self.dim == 3):
# Cell sizes in each direction
self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2]))
self._vol = utils.mkvc(np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2]))
return self._vol
return locals()
_vol = None
@@ -289,12 +289,12 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
elif(self.dim == 2):
area1 = np.outer(np.ones(n[0]+1), vh[1])
area2 = np.outer(vh[0], np.ones(n[1]+1))
self._area = np.r_[mkvc(area1), mkvc(area2)]
self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)]
elif(self.dim == 3):
area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2])))
area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2])))
area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1))))
self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)]
area1 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2])))
area2 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
area3 = np.outer(vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)]
return self._area
return locals()
_area = None
@@ -311,16 +311,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
n = self.n
# Compute edge lengths
if(self.dim == 1):
self._edge = mkvc(vh[0])
self._edge = utils.mkvc(vh[0])
elif(self.dim == 2):
l1 = np.outer(vh[0], np.ones(n[1]+1))
l2 = np.outer(np.ones(n[0]+1), vh[1])
self._edge = np.r_[mkvc(l1), mkvc(l2)]
self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)]
elif(self.dim == 3):
l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1))))
l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2])))
self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)]
l1 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
l2 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
l3 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3)]
return self._edge
return locals()
_edge = None
@@ -410,11 +410,11 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
nF_nE = self.nFv if 'F' in locType else self.nEv
components = [spzeros(loc.shape[0], n) for n in nF_nE]
components[ind] = interpmat(loc, *self.getTensor(locType))
components = [utils.spzeros(loc.shape[0], n) for n in nF_nE]
components[ind] = utils.interpmat(loc, *self.getTensor(locType))
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
Q = interpmat(loc, *self.getTensor(locType))
Q = utils.interpmat(loc, *self.getTensor(locType))
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
return Q
SimPEG/regularization/Regularization.py
+21 -20
View File
@@ -1,9 +1,21 @@
from SimPEG.utils import sdiag, count, timeIt, setKwargs
import numpy as np
from SimPEG import utils, np
class Regularization(object):
"""docstring for Regularization"""
__metaclass__ = utils.Save.Savable
alpha_s = 1e-6
alpha_x = 1.0
alpha_y = 1.0
alpha_z = 1.0
counter = None
def __init__(self, mesh, **kwargs):
utils.setKwargs(self, **kwargs)
self.mesh = mesh
@property
def mref(self):
if getattr(self, '_mref', None) is None:
@@ -16,43 +28,32 @@ class Regularization(object):
@property
def Ws(self):
if getattr(self,'_Ws', None) is None:
self._Ws = sdiag(self.mesh.vol)
self._Ws = utils.sdiag(self.mesh.vol)
return self._Ws
@property
def Wx(self):
if getattr(self, '_Wx', None) is None:
self._Wx = self.mesh.cellGradx*sdiag(self.mesh.vol)
self._Wx = self.mesh.cellGradx*utils.sdiag(self.mesh.vol)
return self._Wx
@property
def Wy(self):
if getattr(self, '_Wy', None) is None:
self._Wy = self.mesh.cellGrady*sdiag(self.mesh.vol)
self._Wy = self.mesh.cellGrady*utils.sdiag(self.mesh.vol)
return self._Wy
@property
def Wz(self):
if getattr(self, '_Wz', None) is None:
self._Wz = self.mesh.cellGradz*sdiag(self.mesh.vol)
self._Wz = self.mesh.cellGradz*utils.sdiag(self.mesh.vol)
return self._Wz
alpha_s = 1e-6
alpha_x = 1.0
alpha_y = 1.0
alpha_z = 1.0
counter = None
def __init__(self, mesh, **kwargs):
setKwargs(self, **kwargs)
self.mesh = mesh
def pnorm(self, r):
return 0.5*r.dot(r)
@timeIt
@utils.timeIt
def modelObj(self, m):
mresid = m - self.mref
@@ -67,7 +68,7 @@ class Regularization(object):
return mobj
@timeIt
@utils.timeIt
def modelObjDeriv(self, m):
"""
@@ -103,7 +104,7 @@ class Regularization(object):
return mobjDeriv
@timeIt
@utils.timeIt
def modelObj2Deriv(self):
mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
SimPEG/utils/Save.py
+56 -18
View File
@@ -5,7 +5,7 @@ import re
try:
import h5py
except Exception, e:
print 'Warning: SimPEG table needs h5py to be installed.'
print 'Warning: SimPEG.utils.Save needs h5py to be installed.'
SAVEABLES = {}
@@ -47,7 +47,10 @@ class SimPEGTable:
# Create a new inversion anytime this is run.
def _startup_hdf5_inv(invObj, m0):
invObj._invNode = self.inversions.addGroup('%d'%self.inversions.numChildren)
node = self.inversions.addGroup('%d'%self.inversions.numChildren)
saveSavable(invObj,node.addGroup('rebuild'))
results = node.addGroup('results')
invObj._invNode = results
invObj.hook(_startup_hdf5_inv, overwrite=True)
# At the start of every iteration we will create a inversion iteration node.
@@ -196,18 +199,25 @@ class hdf5InversionGroup(hdf5Group):
hdf5Group.__init__(self, T, groupNode)
self.childClass = hdf5Inversion
class hdf5Inversion(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5InversionGroup
self.childClass = hdf5InversionIteration
self.childClass = hdf5InversionResults
def rebuild(self):
return loadSavable(self['rebuild'])
class hdf5InversionResults(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5Inversion
self.childClass = hdf5InversionIteration
class hdf5InversionIteration(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5Inversion
self.parentClass = hdf5InversionResults
@@ -225,38 +235,61 @@ class Savable(type):
return newClass
def saveSavable(obj, group):
def saveSavable(obj, group, debug=False):
"""
This creates softlinks if _savable exists in children object.
The first object is always created.
"""
assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.'
def doSave(grp, name, val):
if debug: print name, val
if type(val.__class__) is Savable:
subgrp = grp.addGroup(name)
saveInitArgs(val, subgrp)
elif type(val) is np.ndarray:
grp.setArray(name, val)
link = getattr(val,'_savable',None)
if link is not None:
group.node[name] = h5py.SoftLink(link.path)
if debug: 'Created a softlink path to %s' % link.path
else:
subgrp = grp.addGroup(name)
saveSavable(val, subgrp, debug=debug)
elif type(val) in [list, tuple]:
# Split up, and save each element
for i, v in enumerate(val):
doSave(grp, name + '[%d]'%i, v)
elif type(val) is np.ndarray:
grp.setArray(name, val)
elif val is None:
grp.attrs[name] = 'None'
else:
# just try saving it as an attr
grp.attrs[name] = val
try:
grp.attrs[name] = val
except Exception, e:
print 'Warning: Could not save %s, problems may arise is loading.' % name
group.attrs['__class__'] = obj.__class__.__name__
for arg in obj._kwargs_init:
doSave(group, '_kwarg_'+arg, obj._kwargs_init[arg])
for i, arg in enumerate(obj._args_init):
doSave(group, '_arg%d'%i, arg)
obj._savable = group
def loadSavable(node):
def loadSavable(node, pointers=None):
"""
pointers allow things that point to the same node in the h5py file to
be returned as the same object, if they have already been created.
"""
if pointers is None: pointers = []
for pointer in pointers:
if pointer._savable.node == node.node: return pointer
args = ([a for a in node.attrs if '_arg' in a] + [a for a in node.children if '_arg' in a])
kwargs = ([a for a in node.attrs if '_kwarg' in a] + [a for a in node.children if '_kwarg' in a])
args.sort(key=utils.Save.natural_keys)
kwargs.sort(key=utils.Save.natural_keys)
args.sort(key=natural_keys)
kwargs.sort(key=natural_keys)
def get(node,key):
if key in node.children: return node[key]
@@ -266,6 +299,7 @@ def loadSavable(node):
for name in args:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
ind = int(name[4:name.index('[')])
if len(ARGS) is ind: # Create the list
@@ -273,7 +307,7 @@ def loadSavable(node):
else:
ARGS[ind].append(val)
elif issubclass(val.__class__,hdf5Group):
ARGS.append(load(val))
ARGS.append(loadSavable(val,pointers=pointers))
else:
ind = int(name[4:])
ARGS.append(val)
@@ -282,6 +316,7 @@ def loadSavable(node):
for name in kwargs:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
key = name[7:name.index('[')]
if key not in KWARGS: # Create the list
@@ -290,15 +325,18 @@ def loadSavable(node):
KWARGS[key].append(val)
elif issubclass(val.__class__,hdf5Group):
key = name[7:]
KWARGS[key] = load(val)
KWARGS[key] = loadSavable(val,pointers=pointers)
else:
key = name[7:]
KWARGS[key] = val
cls = get(node, '__class__')
if cls in SAVEABLES:
return SAVEABLES[cls](*ARGS,**KWARGS)
out = SAVEABLES[cls](*ARGS, **KWARGS)
out._savable = node
pointers.append(out) # Because this is recursive.
return out
else:
print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls
return (cls, ARGS, KWARGS)
return (cls, ARGS, KWARGS, node)