From 9e26543931b6a0e35f5698fd586370d51b0a556a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 14:19:37 -0800 Subject: [PATCH] fixes to the meta classes including soft linking inside the database. These are reflected in the loaded outputs which return the same object. --- SimPEG/inverse/Inversion.py | 2 +- SimPEG/mesh/BaseMesh.py | 1 - SimPEG/mesh/LogicallyOrthogonalMesh.py | 96 +++++++++++++------------ SimPEG/mesh/TensorMesh.py | 60 ++++++++-------- SimPEG/regularization/Regularization.py | 41 +++++------ SimPEG/utils/Save.py | 74 ++++++++++++++----- 6 files changed, 157 insertions(+), 117 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 86a5dd14..19280e0a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -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): diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index ced9bd59..6151d20b 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -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): diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index 5c4a73db..b3dcd095 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -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') diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 65f94b3e..d4ab86b0 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -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 diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 41938152..5c50ff59 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -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 diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index 318723ce..d9b3dd58 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -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)