From 428e34f631384a28f5ca9fc7c709dca2c887132e Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 6 May 2015 13:16:56 -0700 Subject: [PATCH 01/49] removed SrcType --- SimPEG/Survey.py | 3 +-- SimPEG/Tests/test_Survey.py | 40 ++++++++++++++++++------------------- 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 10e4c76f..05a2ad37 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -119,14 +119,13 @@ class BaseSrc(object): rxList = None #: SimPEG Receiver List rxPair = BaseRx - def __init__(self, loc, srcType, rxList, **kwargs): + def __init__(self, loc, rxList, **kwargs): assert type(rxList) is list, 'rxList must be a list' for rx in rxList: assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__ assert len(set(rxList)) == len(rxList), 'The rxList must be unique' self.loc = loc - self.srcType = srcType self.rxList = rxList Utils.setKwargs(self, **kwargs) diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Survey.py index f93395ce..7372a686 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Survey.py @@ -9,14 +9,14 @@ class DataAndFieldsTest(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0]) + Src0 = Survey.BaseSrc(srcLoc, [rxList0]) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1]) + Src1 = Survey.BaseSrc(srcLoc, [rxList1]) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2]) + Src2 = Survey.BaseSrc(srcLoc, [rxList2]) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc(srcLoc, [rxList3]) + Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) self.D = Survey.Data(survey) @@ -122,14 +122,14 @@ class FieldsTest_Alias(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0]) + Src0 = Survey.BaseSrc(srcLoc, [rxList0]) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1]) + Src1 = Survey.BaseSrc(srcLoc, [rxList1]) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2]) + Src2 = Survey.BaseSrc(srcLoc, [rxList2]) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc(srcLoc, [rxList3]) + Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) self.D = Survey.Data(survey) @@ -193,14 +193,14 @@ class FieldsTest_Time(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0]) + Src0 = Survey.BaseSrc(srcLoc, [rxList0]) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1]) + Src1 = Survey.BaseSrc(srcLoc, [rxList1]) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2]) + Src2 = Survey.BaseSrc(srcLoc, [rxList2]) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc(srcLoc, [rxList3]) + Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)]) @@ -295,14 +295,14 @@ class FieldsTest_Time_Aliased(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0]) + Src0 = Survey.BaseSrc(srcLoc, [rxList0]) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1]) + Src1 = Survey.BaseSrc(srcLoc, [rxList1]) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2]) + Src2 = Survey.BaseSrc(srcLoc, [rxList2]) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc(srcLoc, [rxList3]) + Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)]) From 1f58f5ef377179dc682ac40c950bb75811d0844a Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 6 May 2015 15:27:58 -0700 Subject: [PATCH 02/49] loc is now a kwarg --- SimPEG/Tests/test_Survey.py | 40 ++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Survey.py index 7372a686..f5109b71 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Survey.py @@ -9,14 +9,14 @@ class DataAndFieldsTest(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, [rxList0]) + Src0 = Survey.BaseSrc([rxList0], loc=srcLoc) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, [rxList1]) + Src1 = Survey.BaseSrc([rxList1], loc=srcLoc) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, [rxList2]) + Src2 = Survey.BaseSrc([rxList2], loc=srcLoc) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc([rxList3], loc=srcLoc) + Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) self.D = Survey.Data(survey) @@ -122,14 +122,14 @@ class FieldsTest_Alias(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, [rxList0]) + Src0 = Survey.BaseSrc([rxList0],loc=srcLoc) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, [rxList1]) + Src1 = Survey.BaseSrc([rxList1],loc=srcLoc) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, [rxList2]) + Src2 = Survey.BaseSrc([rxList2],loc=srcLoc) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc([rxList3],loc=srcLoc) + Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) self.D = Survey.Data(survey) @@ -193,14 +193,14 @@ class FieldsTest_Time(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, [rxList0]) + Src0 = Survey.BaseSrc([rxList0], loc=srcLoc) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, [rxList1]) + Src1 = Survey.BaseSrc([rxList1], loc=srcLoc) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, [rxList2]) + Src2 = Survey.BaseSrc([rxList2], loc=srcLoc) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc([rxList3], loc=srcLoc) + Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)]) @@ -295,14 +295,14 @@ class FieldsTest_Time_Aliased(unittest.TestCase): XYZ = Utils.ndgrid(x,x,np.r_[0.]) srcLoc = np.r_[0,0,0.] rxList0 = Survey.BaseRx(XYZ, 'exi') - Src0 = Survey.BaseSrc(srcLoc, [rxList0]) + Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc) rxList1 = Survey.BaseRx(XYZ, 'bxi') - Src1 = Survey.BaseSrc(srcLoc, [rxList1]) + Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc) rxList2 = Survey.BaseRx(XYZ, 'bxi') - Src2 = Survey.BaseSrc(srcLoc, [rxList2]) + Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc) rxList3 = Survey.BaseRx(XYZ, 'bxi') - Src3 = Survey.BaseSrc(srcLoc, [rxList3]) - Src4 = Survey.BaseSrc(srcLoc, [rxList0, rxList1, rxList2, rxList3]) + Src3 = Survey.BaseSrc( [rxList3],loc=srcLoc) + Src4 = Survey.BaseSrc( [rxList0, rxList1, rxList2, rxList3],loc=srcLoc) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)]) From 6f487e17d005a1973f3805965ae0298e34b8c218 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 6 May 2015 16:34:24 -0700 Subject: [PATCH 03/49] forgot to stage Survey.py --- SimPEG/Survey.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 05a2ad37..29bdb452 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -119,13 +119,12 @@ class BaseSrc(object): rxList = None #: SimPEG Receiver List rxPair = BaseRx - def __init__(self, loc, rxList, **kwargs): + def __init__(self, rxList, **kwargs): assert type(rxList) is list, 'rxList must be a list' for rx in rxList: assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__ assert len(set(rxList)) == len(rxList), 'The rxList must be unique' - self.loc = loc self.rxList = rxList Utils.setKwargs(self, **kwargs) From ce7b0b24dadf79e43a3428e1c85a66bcd7bd32b7 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 7 May 2015 14:44:57 -0700 Subject: [PATCH 04/49] changed check of out.shape, since if out is a vector (not an array with width 1), out.shape[1] does not exsist --- SimPEG/Problem.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index dff99a97..81ba0511 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -160,9 +160,8 @@ class Fields(object): assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' func = getattr(self, func) out = func(self._fields[alias][:,ind], srcII) - - if out.shape[1] == 1: - out = Utils.mkvc(out) + if not out.size > out.shape[0]: + out = Utils.mkvc(out) return out def __contains__(self, other): From 4bda4b60191b6137f228027e065978ea27aa629b Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 7 May 2015 15:29:41 -0700 Subject: [PATCH 05/49] a better way of checking sizes --- SimPEG/Problem.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 81ba0511..b629d723 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -160,8 +160,8 @@ class Fields(object): assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' func = getattr(self, func) out = func(self._fields[alias][:,ind], srcII) - if not out.size > out.shape[0]: - out = Utils.mkvc(out) + if out.shape[0] == out.size: + out = Utils.mkvc(out) return out def __contains__(self, other): From 24c754cea2cc1d8ba8201d138fb703416fed7432 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 8 May 2015 09:21:43 -0700 Subject: [PATCH 06/49] Working on making fields 1 wide in Problem.mp --- SimPEG/Problem.py | 6 +++--- SimPEG/Tests/test_Survey.py | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index b629d723..7b337977 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -139,7 +139,7 @@ class Fields(object): return self._getField(name, ind) def _setField(self, field, val, name, ind): - if isinstance(val, np.ndarray) and (field.shape[1] == 1 or val.ndim == 1): + if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1): val = Utils.mkvc(val,2) field[:,ind] = val @@ -160,8 +160,8 @@ class Fields(object): assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' func = getattr(self, func) out = func(self._fields[alias][:,ind], srcII) - if out.shape[0] == out.size: - out = Utils.mkvc(out) + if isinstance(out, np.ndarray) and (out.shape[0] == out.size or out.ndim == 1): + out = Utils.mkvc(out,2) return out def __contains__(self, other): diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Survey.py index f93395ce..2a506f8a 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Survey.py @@ -80,9 +80,9 @@ class DataAndFieldsTest(unittest.TestCase): b = np.random.rand(F.mesh.nF,1) F[self.Src0, 'b'] = b - self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b))) + self.assertTrue(np.all(F[self.Src0, 'b'] == b)) - b = np.random.rand(F.mesh.nF) + b = np.random.rand(F.mesh.nF,1) F[self.Src0, 'b'] = b self.assertTrue(np.all(F[self.Src0, 'b'] == b)) @@ -96,8 +96,8 @@ class DataAndFieldsTest(unittest.TestCase): b = np.random.rand(F.mesh.nF, 2) F[[self.Src0, self.Src1],'b'] = b - self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,)) - self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,)) + self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1)) + self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1)) self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0])) self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1])) @@ -158,7 +158,7 @@ class FieldsTest_Alias(unittest.TestCase): e = np.random.rand(F.mesh.nE,1) F[self.Src0, 'e'] = e - self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e))) + self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e)) def f(): F[self.Src0, 'b'] = F[self.Src0, 'b'] From 618f803ffa66acae3874f72ebd1db3924d81415c Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 8 May 2015 10:01:28 -0700 Subject: [PATCH 07/49] Fields returns matrices --- SimPEG/Tests/test_Survey.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Survey.py index e6a9e507..fd3dbf06 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Survey.py @@ -98,8 +98,8 @@ class DataAndFieldsTest(unittest.TestCase): F[[self.Src0, self.Src1],'b'] = b self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1)) self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1)) - self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0])) - self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1])) + self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2))) + self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2))) def test_assertions(self): freq = [self.Src0, self.Src1] From b8118464b87abd49e80ffbf8c1b3b949a05c8c88 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 8 May 2015 14:59:58 -0700 Subject: [PATCH 08/49] Time fields objects return arrays --- SimPEG/Problem.py | 6 ++++-- SimPEG/Tests/test_Survey.py | 12 ++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 7b337977..1787a472 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -216,6 +216,8 @@ class TimeFields(Fields): shape = nP, nSrc, nT if deflate: shape = tuple([s for s in shape if s > 1]) + if len(shape) == 1: + shape = shape + (1,) return shape def _setField(self, field, val, name, ind): @@ -260,8 +262,8 @@ class TimeFields(Fields): out = range(nT) for i, TIND_i in enumerate(timeII): fieldI = pointerFields[:,:,i] - if fieldI.ndim == 2 and fieldI.shape[1] == 1: - fieldI = Utils.mkvc(fieldI) + if fieldI.shape[0] == fieldI.size: + fieldI = Utils.mkvc(fieldI,2) out[i] = func(fieldI, srcII, TIND_i) if out[i].ndim == 1: out[i] = out[i][:,np.newaxis,np.newaxis] diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Survey.py index fd3dbf06..efc1e901 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Survey.py @@ -249,7 +249,7 @@ class FieldsTest_Time(unittest.TestCase): b = np.random.rand(F.mesh.nF,1,nT) F[self.Src0, 'b', 0] = b[:,:,0] - self.assertTrue(np.all(F[self.Src0, 'b', 0] == b[:,0,0])) + self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2))) phi = np.random.rand(F.mesh.nC,2,nT) F[[self.Src0,self.Src1], 'phi'] = phi @@ -265,10 +265,10 @@ class FieldsTest_Time(unittest.TestCase): self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT)) self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:])) self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:])) - self.assertTrue(np.all(F[self.Src0,'b',1] == b[:,0,1])) - self.assertTrue(np.all(F[self.Src1,'b',1] == b[:,1,1])) - self.assertTrue(np.all(F[self.Src0,'b',4] == b[:,0,4])) - self.assertTrue(np.all(F[self.Src1,'b',4] == b[:,1,4])) + self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2))) + self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2))) + self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2))) + self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2))) b = np.random.rand(F.mesh.nF, 2, nT) @@ -344,7 +344,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase): self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] )) self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] )) for t in range(nT): - self.assertTrue(np.all(F[self.Src1, 'e', t] == e[:,1,t] )) + self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) )) b = np.random.rand(F.mesh.nF,nT) F[self.Src0, 'b',:] = b From 0f703739ed1c475046c960b78d0da41028f5dfab Mon Sep 17 00:00:00 2001 From: Lindsey Date: Thu, 14 May 2015 14:12:52 -0700 Subject: [PATCH 09/49] Changes LRM to Curvilinear --- .../{LogicallyRectMesh.py => CurvilinearMesh.py} | 12 ++++++------ SimPEG/Mesh/__init__.py | 2 +- SimPEG/Tests/TestUtils.py | 6 +++--- ..._LogicallyRectMesh.py => test_CurvilinearMesh.py} | 6 +++--- SimPEG/Tests/test_innerProductDerivs.py | 4 ++-- docs/api_Mesh.rst | 4 ++-- docs/api_MeshCode.rst | 2 +- 7 files changed, 18 insertions(+), 18 deletions(-) rename SimPEG/Mesh/{LogicallyRectMesh.py => CurvilinearMesh.py} (98%) rename SimPEG/Tests/{test_LogicallyRectMesh.py => test_CurvilinearMesh.py} (97%) diff --git a/SimPEG/Mesh/LogicallyRectMesh.py b/SimPEG/Mesh/CurvilinearMesh.py similarity index 98% rename from SimPEG/Mesh/LogicallyRectMesh.py rename to SimPEG/Mesh/CurvilinearMesh.py index b2dc6f55..0fbf78a3 100644 --- a/SimPEG/Mesh/LogicallyRectMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -10,9 +10,9 @@ 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 LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): +class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): """ - LogicallyRectMesh is a mesh class that deals with logically rectangular meshes. + CurvilinearMesh is a mesh class that deals with logically rectangular meshes. Example of a logically rectangular mesh: @@ -21,7 +21,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): from SimPEG import Mesh, Utils X, Y = Utils.exampleLrmGrid([3,3],'rotate') - M = Mesh.LogicallyRectMesh([X, Y]) + M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ @@ -343,7 +343,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): from SimPEG import Mesh, Utils X, Y = Utils.exampleLrmGrid([3,3],'rotate') - M = Mesh.LogicallyRectMesh([X, Y]) + M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ @@ -435,9 +435,9 @@ if __name__ == '__main__': dee3 = True if dee3: X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False) - M = LogicallyRectMesh([X, Y, Z]) + M = CurvilinearMesh([X, Y, Z]) else: X, Y = Utils.ndgrid(h1, h2, vector=False) - M = LogicallyRectMesh([X, Y]) + M = CurvilinearMesh([X, Y]) print M.r(M.normals, 'F', 'Fx', 'V') diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index deb11321..9aadfcba 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,5 +1,5 @@ from TensorMesh import TensorMesh from CylMesh import CylMesh -from LogicallyRectMesh import LogicallyRectMesh +from CurvilinearMesh import CurvilinearMesh from TreeMesh import TreeMesh from BaseMesh import BaseMesh diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index dbcb92b6..57e17f06 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag, diagEst from SimPEG import Utils -from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh +from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh import numpy as np import scipy.sparse as sp import unittest @@ -126,10 +126,10 @@ class OrderTest(unittest.TestCase): raise Exception('Lom not supported for 1D') elif self.meshDimension == 2: X, Y = Utils.exampleLrmGrid([nc, nc], kwrd) - self.M = LogicallyRectMesh([X, Y]) + self.M = CurvilinearMesh([X, Y]) elif self.meshDimension == 3: X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd) - self.M = LogicallyRectMesh([X, Y, Z]) + self.M = CurvilinearMesh([X, Y, Z]) return 1./nc def getError(self): diff --git a/SimPEG/Tests/test_LogicallyRectMesh.py b/SimPEG/Tests/test_CurvilinearMesh.py similarity index 97% rename from SimPEG/Tests/test_LogicallyRectMesh.py rename to SimPEG/Tests/test_CurvilinearMesh.py index 5d223f68..132b6f03 100644 --- a/SimPEG/Tests/test_LogicallyRectMesh.py +++ b/SimPEG/Tests/test_CurvilinearMesh.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG.Mesh import TensorMesh, LogicallyRectMesh +from SimPEG.Mesh import TensorMesh, CurvilinearMesh from SimPEG.Utils import ndgrid @@ -13,10 +13,10 @@ class BasicLRMTests(unittest.TestCase): gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] X, Y = ndgrid(gridIt([a, b]), vector=False) self.TM2 = TensorMesh([a, b]) - self.LRM2 = LogicallyRectMesh([X, Y]) + self.LRM2 = CurvilinearMesh([X, Y]) X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) self.TM3 = TensorMesh([a, b, c]) - self.LRM3 = LogicallyRectMesh([X, Y, Z]) + self.LRM3 = CurvilinearMesh([X, Y, Z]) def test_area_3D(self): test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) diff --git a/SimPEG/Tests/test_innerProductDerivs.py b/SimPEG/Tests/test_innerProductDerivs.py index 4b7d63cd..27979445 100644 --- a/SimPEG/Tests/test_innerProductDerivs.py +++ b/SimPEG/Tests/test_innerProductDerivs.py @@ -9,7 +9,7 @@ class TestInnerProductsDerivs(unittest.TestCase): def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False): if meshType == 'LRM': hRect = Utils.exampleLrmGrid(h,'rotate') - mesh = Mesh.LogicallyRectMesh(hRect) + mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': mesh = Mesh.TreeMesh(h) elif meshType == 'Tensor': @@ -26,7 +26,7 @@ class TestInnerProductsDerivs(unittest.TestCase): def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False): if meshType == 'LRM': hRect = Utils.exampleLrmGrid(h,'rotate') - mesh = Mesh.LogicallyRectMesh(hRect) + mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': mesh = Mesh.TreeMesh(h) elif meshType == 'Tensor': diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index c698cf67..7bf398b1 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -29,7 +29,7 @@ the implementations. tM = Mesh.TensorMesh(sz) qM = Mesh.TreeMesh(sz) qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0) - rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) fig, axes = plt.subplots(1,3,figsize=(14,5)) opts = {} @@ -38,7 +38,7 @@ the implementations. qM.plotGrid(ax=axes[1], **opts) axes[1].set_title('TreeMesh') rM.plotGrid(ax=axes[2], **opts) - axes[2].set_title('LogicallyRectMesh') + axes[2].set_title('CurvilinearMesh') plt.show() diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 13ceffba..9fb52e53 100644 --- a/docs/api_MeshCode.rst +++ b/docs/api_MeshCode.rst @@ -21,7 +21,7 @@ Tree Mesh Logically Rectangular Mesh ========================== -.. automodule:: SimPEG.Mesh.LogicallyRectMesh +.. automodule:: SimPEG.Mesh.Curvilinear :show-inheritance: :members: :undoc-members: From 46b2e11ef8d57c4312216acd014d2154e525accb Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 14 May 2015 23:35:38 -0700 Subject: [PATCH 10/49] LRM --> Curv --- SimPEG/Mesh/CurvilinearMesh.py | 4 +- SimPEG/Mesh/InnerProducts.py | 16 ++-- SimPEG/Tests/TestUtils.py | 2 +- SimPEG/Tests/test_CurvilinearMesh.py | 104 ++++++++++----------- SimPEG/Tests/test_innerProduct.py | 4 +- SimPEG/Tests/test_innerProductDerivs.py | 116 ++++++++++++------------ SimPEG/Tests/test_operators.py | 6 +- docs/SimPEGFrameworkRevised.png | Bin 0 -> 31193 bytes docs/api_InnerProducts.rst | 4 +- docs/api_Utils.rst | 2 +- 10 files changed, 129 insertions(+), 129 deletions(-) create mode 100644 docs/SimPEGFrameworkRevised.png diff --git a/SimPEG/Mesh/CurvilinearMesh.py b/SimPEG/Mesh/CurvilinearMesh.py index 0fbf78a3..e4eeda3e 100644 --- a/SimPEG/Mesh/CurvilinearMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -27,7 +27,7 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): __metaclass__ = Utils.SimPEGMetaClass - _meshType = 'LRM' + _meshType = 'Curv' def __init__(self, nodes): assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" @@ -38,7 +38,7 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i) assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" - assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh." + assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh." BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 2cf48a7c..42a08702 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -328,7 +328,7 @@ class InnerProducts(object): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LRM': + if M._meshType == 'Curv': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') @@ -353,7 +353,7 @@ class InnerProducts(object): PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF)) - if M._meshType == 'LRM': + if M._meshType == 'Curv': I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]), getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy])) PXX = I2x2 * PXX @@ -376,7 +376,7 @@ class InnerProducts(object): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LRM': + if M._meshType == 'Curv': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') fN3 = M.r(M.normals, 'F', 'Fz', 'M') @@ -410,7 +410,7 @@ class InnerProducts(object): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]), getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]), getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ])) @@ -432,7 +432,7 @@ class InnerProducts(object): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LRM': + if M._meshType == 'Curv': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') @@ -452,7 +452,7 @@ class InnerProducts(object): PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]), getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j])) PXX = I2x2 * PXX @@ -466,7 +466,7 @@ class InnerProducts(object): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LRM': + if M._meshType == 'Curv': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') eT3 = M.r(M.tangents, 'E', 'Ez', 'M') @@ -495,7 +495,7 @@ class InnerProducts(object): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]), getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]), getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k])) diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index 57e17f06..c6aa5bc4 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -115,7 +115,7 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h - elif 'LRM' in self._meshType: + elif 'Curv' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' elif 'rotate' in self._meshType: diff --git a/SimPEG/Tests/test_CurvilinearMesh.py b/SimPEG/Tests/test_CurvilinearMesh.py index 132b6f03..42e3d877 100644 --- a/SimPEG/Tests/test_CurvilinearMesh.py +++ b/SimPEG/Tests/test_CurvilinearMesh.py @@ -4,7 +4,7 @@ from SimPEG.Mesh import TensorMesh, CurvilinearMesh from SimPEG.Utils import ndgrid -class BasicLRMTests(unittest.TestCase): +class BasicCurvTests(unittest.TestCase): def setUp(self): a = np.array([1, 1, 1]) @@ -13,91 +13,91 @@ class BasicLRMTests(unittest.TestCase): gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] X, Y = ndgrid(gridIt([a, b]), vector=False) self.TM2 = TensorMesh([a, b]) - self.LRM2 = CurvilinearMesh([X, Y]) + self.Curv2 = CurvilinearMesh([X, Y]) X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) self.TM3 = TensorMesh([a, b, c]) - self.LRM3 = CurvilinearMesh([X, Y, Z]) + self.Curv3 = CurvilinearMesh([X, Y, Z]) def test_area_3D(self): test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) - self.assertTrue(np.all(self.LRM3.area == test_area)) + self.assertTrue(np.all(self.Curv3.area == test_area)) def test_vol_3D(self): test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) - np.testing.assert_almost_equal(self.LRM3.vol, test_vol) + np.testing.assert_almost_equal(self.Curv3.vol, test_vol) self.assertTrue(True) # Pass if you get past the assertion. def test_vol_2D(self): test_vol = np.array([1, 1, 1, 2, 2, 2]) - t1 = np.all(self.LRM2.vol == test_vol) + t1 = np.all(self.Curv2.vol == test_vol) self.assertTrue(t1) def test_edge_3D(self): test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) - t1 = np.all(self.LRM3.edge == test_edge) + t1 = np.all(self.Curv3.edge == test_edge) self.assertTrue(t1) def test_edge_2D(self): test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) - t1 = np.all(self.LRM2.edge == test_edge) + t1 = np.all(self.Curv2.edge == test_edge) self.assertTrue(t1) def test_tangents(self): - T = self.LRM2.tangents - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy))) + T = self.Curv2.tangents + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy))) - T = self.LRM3.tangents - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx))) + T = self.Curv3.tangents + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz))) def test_normals(self): - N = self.LRM2.normals - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy))) + N = self.Curv2.normals + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy))) - N = self.LRM3.normals - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx))) + N = self.Curv3.normals + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz))) def test_grid(self): - self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC)) - self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN)) - self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx)) - self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy)) - self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx)) - self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy)) + self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC)) + self.assertTrue(np.all(self.Curv2.gridN == self.TM2.gridN)) + self.assertTrue(np.all(self.Curv2.gridFx == self.TM2.gridFx)) + self.assertTrue(np.all(self.Curv2.gridFy == self.TM2.gridFy)) + self.assertTrue(np.all(self.Curv2.gridEx == self.TM2.gridEx)) + self.assertTrue(np.all(self.Curv2.gridEy == self.TM2.gridEy)) - self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC)) - self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN)) - self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx)) - self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy)) - self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz)) - self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx)) - self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy)) - self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz)) + self.assertTrue(np.all(self.Curv3.gridCC == self.TM3.gridCC)) + self.assertTrue(np.all(self.Curv3.gridN == self.TM3.gridN)) + self.assertTrue(np.all(self.Curv3.gridFx == self.TM3.gridFx)) + self.assertTrue(np.all(self.Curv3.gridFy == self.TM3.gridFy)) + self.assertTrue(np.all(self.Curv3.gridFz == self.TM3.gridFz)) + self.assertTrue(np.all(self.Curv3.gridEx == self.TM3.gridEx)) + self.assertTrue(np.all(self.Curv3.gridEy == self.TM3.gridEy)) + self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz)) if __name__ == '__main__': diff --git a/SimPEG/Tests/test_innerProduct.py b/SimPEG/Tests/test_innerProduct.py index bbaf4b14..0a5ff809 100644 --- a/SimPEG/Tests/test_innerProduct.py +++ b/SimPEG/Tests/test_innerProduct.py @@ -7,7 +7,7 @@ from SimPEG import Utils class TestInnerProducts(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] + meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] meshDimension = 3 meshSizes = [16, 32] @@ -154,7 +154,7 @@ class TestInnerProducts(OrderTest): class TestInnerProducts2D(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] + meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] meshDimension = 2 meshSizes = [4, 8, 16, 32, 64, 128] diff --git a/SimPEG/Tests/test_innerProductDerivs.py b/SimPEG/Tests/test_innerProductDerivs.py index 27979445..6eb561c1 100644 --- a/SimPEG/Tests/test_innerProductDerivs.py +++ b/SimPEG/Tests/test_innerProductDerivs.py @@ -7,7 +7,7 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False): - if meshType == 'LRM': + if meshType == 'Curv': hRect = Utils.exampleLrmGrid(h,'rotate') mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': @@ -24,7 +24,7 @@ class TestInnerProductsDerivs(unittest.TestCase): return checkDerivative(fun, sig, num=5, plotIt=False) def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False): - if meshType == 'LRM': + if meshType == 'Curv': hRect = Utils.exampleLrmGrid(h,'rotate') mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': @@ -137,65 +137,65 @@ class TestInnerProductsDerivs(unittest.TestCase): - def test_FaceIP_2D_float_LRM(self): - self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM')) - def test_FaceIP_3D_float_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM')) - def test_FaceIP_2D_isotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM')) - def test_FaceIP_3D_isotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM')) - def test_FaceIP_2D_anisotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM')) - def test_FaceIP_3D_anisotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM')) - def test_FaceIP_2D_tensor_LRM(self): - self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM')) - def test_FaceIP_3D_tensor_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM')) + def test_FaceIP_2D_float_Curv(self): + self.assertTrue(self.doTestFace([10, 4],0, False, 'Curv')) + def test_FaceIP_3D_float_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Curv')) + def test_FaceIP_2D_isotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4],1, False, 'Curv')) + def test_FaceIP_3D_isotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Curv')) + def test_FaceIP_2D_anisotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4],2, False, 'Curv')) + def test_FaceIP_3D_anisotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Curv')) + def test_FaceIP_2D_tensor_Curv(self): + self.assertTrue(self.doTestFace([10, 4],3, False, 'Curv')) + def test_FaceIP_3D_tensor_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Curv')) - def test_FaceIP_2D_float_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM')) - def test_FaceIP_3D_float_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM')) - def test_FaceIP_2D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM')) - def test_FaceIP_3D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM')) - def test_FaceIP_2D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM')) - def test_FaceIP_3D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM')) + def test_FaceIP_2D_float_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],0, True, 'Curv')) + def test_FaceIP_3D_float_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Curv')) + def test_FaceIP_2D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],1, True, 'Curv')) + def test_FaceIP_3D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Curv')) + def test_FaceIP_2D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],2, True, 'Curv')) + def test_FaceIP_3D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Curv')) - def test_EdgeIP_2D_float_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM')) - def test_EdgeIP_3D_float_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM')) - def test_EdgeIP_2D_isotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM')) - def test_EdgeIP_3D_isotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM')) - def test_EdgeIP_2D_anisotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM')) - def test_EdgeIP_3D_anisotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM')) - def test_EdgeIP_2D_tensor_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM')) - def test_EdgeIP_3D_tensor_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM')) + def test_EdgeIP_2D_float_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, 'Curv')) + def test_EdgeIP_3D_float_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Curv')) + def test_EdgeIP_2D_isotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, 'Curv')) + def test_EdgeIP_3D_isotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Curv')) + def test_EdgeIP_2D_anisotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, 'Curv')) + def test_EdgeIP_3D_anisotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Curv')) + def test_EdgeIP_2D_tensor_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],3, False, 'Curv')) + def test_EdgeIP_3D_tensor_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Curv')) - def test_EdgeIP_2D_float_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM')) - def test_EdgeIP_3D_float_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM')) - def test_EdgeIP_2D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM')) - def test_EdgeIP_3D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM')) - def test_EdgeIP_2D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM')) - def test_EdgeIP_3D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM')) + def test_EdgeIP_2D_float_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],0, True, 'Curv')) + def test_EdgeIP_3D_float_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Curv')) + def test_EdgeIP_2D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],1, True, 'Curv')) + def test_EdgeIP_3D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Curv')) + def test_EdgeIP_2D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],2, True, 'Curv')) + def test_EdgeIP_3D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Curv')) diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 9a465fcd..416689bc 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -4,7 +4,7 @@ from TestUtils import OrderTest import matplotlib.pyplot as plt #TODO: 'randomTensorMesh' -MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] +MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -38,7 +38,7 @@ class TestCurl(OrderTest): curlE_ana = self.M.projectFaceVector(Fc) curlE = self.M.edgeCurl.dot(E) - if self._meshType == 'rotateLRM': + if self._meshType == 'rotateCurv': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2) @@ -229,7 +229,7 @@ class TestFaceDiv3D(OrderTest): divF = self.M.faceDiv.dot(F) divF_ana = call3(sol, self.M.gridCC) - if self._meshType == 'rotateLRM': + if self._meshType == 'rotateCurv': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2) diff --git a/docs/SimPEGFrameworkRevised.png b/docs/SimPEGFrameworkRevised.png new file mode 100644 index 0000000000000000000000000000000000000000..18ce7a304d14227626f99c011d32873007b071bb GIT binary patch literal 31193 zcmc$GcT`hd({B7nfcAkK4JH??=l?aI0^!R z7}Qj6-v@yXVL%{A7ac9|jfKydG7!k`z1rC+cz&AzN(DlwP{M(dS9GeFz8DlfO$QaZZVbm`gnDHJ;DoR3TF&q zAE#GHY&;?=jsNb8Il>BE`bq-~0{&|lf3cBUDXWwFb5ut$#!)_M&AuDkstvT4*bP+h?0wx=zNdG4_>kGlKi}IU@gD}B1lMaZ zkIC$5xUb&7L(fycL&sHbARDr+=Q}n=Yz#eLpx7IpAROhckWHoPCcOFi4F6Xl;UJo= zX7>lWX*zy-N0%B;Qh%enOt+MK`VsY`tt#ABnKIbRS?8M&zug_)Y)M(o#0D&vS=Qv+ z$+4ASt*}Gk=1{F+C@)GeM^?9Bmp>4d93QcZdWSMOiVH1T%yxnv#UU)yTl-x6X?^*e zb~PUQVvt4KAZ*8rRR$xvrZoRhc_JZf0KHQhzwDt=?nKQY*obGLvY9 zLaDb0YN6|)cfqLmK2rP#F|%zyS9v)`h0SsEhtv1P&25{RHl_3(g~;U5YqNrO{~+ViMOd6n(;oguobtk<|km8FD&lebDeO6ZfSIFFtp5^ zJUo^!!#I{T5n!r{|2V%_iJ*6v$r$RH-!#@~U(EDE_U+jqyCa?z#wqf46W+yjW)e|0 z-F3;RntL9RPLCyOW75KPN}&U7E0y!c>7A3LJ44eiJl{OwFetgP8`IXT;;lUKtepX1 z`EoL(HX$?9>fBmD$LXuNBN_~()MGD?YYFhZIwxT$vRfP-8?DyoUrr2j9Q0_se(|k5 z?i_|UrZ}!7=%dQFOKwn}jm;smOONM^A$wQmiNbYkJ9YD=lFjA1bW7$e-KqJKhy+M2 z`;I(HuhC_LH)=j=Uxjw@8>QG))ik;v8hm?_txBQvXVjpXcbv<4V8a5mzUOY9>pw(-gbB2s6n6|as8 zVxB0NI#<5ECuH{ygq6&+mWpk=50_H%VBEORH{6TD^^5Gvv|6Zny~zB)@j+g%bpGMa z;=|9ByOQk+3`YV==4iI8%Ik^<+h@UKQ;#6S5uM>mbGjv!j*{U!V*q_1c0TZ!{o88O zNfBxCt_T!X?oo7Wy6@PPs=??Z?~E+$ znVQS@WZ!_0$EC+^LH2Sj5hjfhp%UEs;>4b09Ny0D{D+~<)nw8&l0w$xRDE2#uYP8g zFk^C1RJmI5=@>d>mKToqKHS!1eGQcmrkm(UX5#IBtYh3WGahK^)=YG2H+Jfv} zKkhhxpfSX@87RBu+ddHDd~E;aj^=E)je=RSOR<^6)nF)-CSYv+RQrVIT@Qiq&sJkPuFa2;eZ&u zOAH}L&(ttjxqnn?o1>ZgmPqHS$LP8)vZkago5zDsI?)W{g6w^nGE#27@%tjiZSvN7 z>gK5D$H(JL1V7Yh^{;K2_C~+Ea4=v5ub#!cqZ`pYXWBVyDu1ZfxU^%YZ13TbuX#4k zMXbN}5;c-!d6Pp+v@-+sA&Ms;rklr?ZpUpg7nE8&#)GVo zq(*BN{4r7V@=Mn(Teza&TSVTVbMHjZJJ2D3Dcig7d_oX9ctFgA4|D&RTKHt73#plx(#7F7E z^Er8u%CP63XQdT#RxQcVkUbPrkG=P-&o4ZU9kPVjRvUQztdmbE^n%dC-T2NqHOO&* zfK+a92y;uEAeyRS3v5R=J?X&$(4l(NJJ79V&7A`0#Ci9fB4^b(rFQ?S2TuH!yV1c@1`(HCf_DsCZ%=poAY-CyHJcU9 zJVe?Xt9#Un%a6>%BnuUL_|VBlMKU3?r5X%=hCD8#E^is11;sTePJdnC>;9bB7Z}?I z>vO}s2ts_JwL~TJubajEC{*GC6i=AjZRN`mVa1DXQBqeHZ3z9V{2y5duXmp6CW$yr z_-hP}??`oTC}b|+j_u;sO-3hkcLRr4sHzI;f~Vi6{srx8E(pHA)}YHOzVZA1asMaj zt9zxCDWr|_UMs9}9ZMmgjo6gv_=1^~=-j-=5ryQ&4#geAVLe-=AW9Bw*`JI+l0!$$ z=jr!{@C8f9E^lk9E-6|B)=f8`KDhH$Dfefgt+ivW^BWu1#|{QI3`VANrptiIw)qsx zVQJsLZePHXH%5pB$Ci>*uTcp~b?eMLwl8+jSqu$3^()L_0gnI4-48 zcA!FPD|XEIM*iJ7fd0z&w-14u^x+kc8@nLTz?);E8aS6X_l6=bVl;$AmOvn}3qL?` zR%ZNaUFDB(hq}c4?H~zyL#TrHbe-1>WD|ZpeW`sNvB0Y8vD>OQvnHv%VjV#Mx$&)S6+UCgtA=Gx0al zexJqY8JnAGY^PjJQ=`}{n{@PZ8NNfa^;}0@=bCun{pN~T)|B`u^eGVhnEEvkKIfS0 z$&XUM@E+exUY6+(M|yQ?myhhOxdz9MyiQJWw{m~sR;+cG=?9leXwYq&f(vXiam~N ziuq%M}KO~2Dyepe;<}Wx^nxQyZ{?NIC z>C?(U{88}Ou_FP@Qr(39*H|ngABNA#M~-WLYDBw4dQ=y_3eErI-`0Mn5^?%|zyt7@ z_KJ+LN9o6tUD<)}Af{8TFRiLt}8M3rZ@#nnZ8co1e<)W9}4P?@sASntun`%dnvkA9gSM(}6qh z_==9wY$dlY<%oznK}?1G-repz9zFE&sO!^~&|{KT!*}30-8dJHm04{!)~vQz6H!0F zDeph@I(Si*MpkUM{8x>7QMwmXY0Aq1l*XKyCIY98XvQjpzF;+MwEV&td-ff6SxoYk zai=p2w(DF(XjXDd?kqDUksB5Zf^jL4gMvZLggj!_UXY3EUNPi`B0jbHzT?NX*i*np zH@)`(4EiXDaf!SI&;RwsZ6(V-(WCd!mVlGBMA*QDW8F*_0(41?EtJs2L)hXXzTl(T z8fHH1dy#&)4L+4ecq$%Q3^K%?f@%11BI+aR%0d^wf|x(+c#&qyngv_#B>(%Y4PUFl zo7rC3-hwWj9)*j#ZHS&a-W_6}Xy}6TaXyTbauvcWu4o{9Po@l=M+kD$R5XDjPB7t< zx)23aE%ZJS-$ zw)u@7mWdLwgdr6Ghq!Nd5#KyndSCMOhbwuChK;20$7&6yy9{=2mwoq)FcLfU;o4Ii zt6f^7(l@pvnG}_wzV)(&m&Gu4FSbshUnbacNxjk$N3C?z<0Q3`*_HR1vZ$51I=7UH zZZa}^FHtl;&2g&^E0koWvwg*E-l4js)IrS%%g{nd5`no$rrP#i={)fki1VqWQ2#+R zsi7a>r&%9SFZbp1Z+~wHgrQdVI=>@G!?=;1j#6B~zIT1&K0;nqtlLJCJJ8A9lZQc& z8=fC1^Mi;E9C;NtKNhr+B2V?Vudak)`Po5qOI=|o5&>47&2kEG*bM`X2+E9BhdJPI zl?8Ut#|7?(9R-G4eLq)jc`X182l48{DW1ctcc|W-y~-9A;B*q$Xf#_!hmD>WAeR;J_tDo|G^p_opkjYt0HZiwf$U-vfYiW#Q~)n))c^rokM1&<_U6!} zz#mEvopoAap$1LV1&O}yM^BvzxEU0u=10FDO#tj1;&>IE@MjSH<(ori#lZ*G=kE(( zl*3s_rj|M6EC$F_0Q;|7P#k<9BliD9$9qj60-}aLa0y{xT7U{FS2qtWHy3w>7lsVX zH!hfe=XK+CGjRZ@|3heMW|`=z-Mny@x=OExGDkp%{&05rC}i&%lsF?gVZefCKNdRD ztS(*zQ1SQAwY??ZNoC&y|K?P?6omXKa`i_zkTDqxj}ji+7&jNs4qtD#g&X$Rder$2 zhBMar0^;tMI#^vm973WPpKghiLEESvTbGzx1P&HErrjFt`84y-FayHdAf!`x0J_g=>ysUtS#n?o<6cDuuml^oB}gy-{hEw8I?e%qukPUqha zj78Plf;T^e=fCPNA7qs_05s@8>Ng;%j3l0-3pH-uBVkra(tCdA=ZR`babD85?7|6? z9U;WivFwEmRULC`-ssGdrE6=UUGc-=Fc%vbG7C!^rXTk^m$|*Y=PU-%O?=Q)<*Cd+ zy;~!KO=o%ZrWn*b*TQXllo~jK$yaXBk@z}~9g43@n3~BesMm=vZ{^R$Dp9gKq^G|+ zi}QpMkq?v4Vw%rk%rf0CiW*K0%Rt^Rv*5!#wb%)BoKzjTVZhMR&&?TK^J~)op6dpU z^!1|?9+hsdmit{@rRSsxw`1mLYv`vg0Z#RdC8sfZt}2E%Z3|j zdRO}6dN@B%@6`@B3n}~fw9ay`Kui!aYl=rYEAr0J_~mKohAq9Puyl`+!_rZVGkXd? z-~H6!yCRaWCfzTh9cHrI>pklQRB8vMEK1Xw2e;3-9B)#iA%;w!$iYZI*X!~ad+zy+xx=$+YIe^duvKus%pL1 z@QQ;wu@w2<45aqpIs?^m@f^Ee3a#cmp~&wcFx&|ZRr#i_tn7V$o;TYl?CZmgP1E)b zDT&V=2)*>S%xYotkQZ261J<*3ywT*@?QfzvKQpqN6w+Q0l{e`lv|~X8W^9c#bH4~! z!2pa(o*GW$kQvqfW=Q7g?I;TiA-0>7G(;h$@)l!dY@6wq z@4fdF+djxnCw)L&n0Pj4QUkJ}C{lchdz2-{4f0Nd1%;wL^B39E;VL?gs$&`qd6tj2 zDrvSf`sivAnQP80&t|JWWw%BMmZ&@$WDhXD_VaYf8m_&Z^ub7+S7PvIcuAsY>0}&7 z-ZLXD8_7>XlLhpyVq(2?FNOOn&L3Nv2q~#fZ*>rWa?XVGwL4aYT2dHFFH|q944YOP zJm^Ovf8puk){|4X61u${Py-vQd^Ela!!9?-y$^~6yS?`AJWa3RW5ZWQq&?ta>-KWG zG>?$4sa6ro@&L)kE{M1iT^!MZ-H0Q|1m+5S;c=vV9o zo&*pvE94P3wr7*bfkLM`O(@lWmh33-J8xoB#Z+c5{mZ7lI_1|>>_%E6t(8pDJ9e(^ zN#i1yditfueD8T5POb&Im0Et7R_QOj2Ycblx8`!exnh38=@W6PLx`jHr!E}vq)5pB9|AfXxq9xMCCA66jp@v{8;! z!NC#KER4~@GFcp3sQVnnNweJxax3!XoDFz}cszsnWOp{`gka8@@@E$!#p#xOwuke= zjJz#ENgqv7y)9Gz?O$@obp>BLA#2(MyfZ9$6HIQ`Ri`_8G2yG`SIf7Vr#bi~+32V^P}itmw{P~erpT1VkGe9APUsIPiXcFY2{7LKekVjBBR|7;ihj`= zw{G84?D16iULBE>fl8>>mOzH?=!$zixW7?(3!XBTUza{A{q&m_hC~?D+*vT22tJDe z+^g^LB+c&4;^<$6y!@CKAx#hJBihFXIE~LcWoKO3;d4R@uheXKy!3_Ht8Vr@QpH+w zygC&=)70>OFdaL471Ij!EYbCyRV)n+k!9X=qi(eIbknP%kmXS#eIm zcBLvd$jx$;$3!$K6cz{WTFsSu*&DLXh_mKGdYE6kCg+mJt&G-Wjy zey83o*IM)RuuOtej=$}ayOh+0SM)zT`xjErS6Z zr%Bd?M}+kE5E3?x=cTjI=v(5k2Gi(Q-HLsB5j;>W$F{Nk606?Z832vv-r? zCZ*b)`{%8@FGkaHoaV$$+EzVJg>YuXUO;4lLI3pVval<@rN>^y)v)dZFeY>rTM&Ig_aeGlAUd5)K@c_DmXd3k&$X}<=>kzjy_Dp z+FyPjEt6K0x~_MuBxA3cX$7sf3VWV=XCVpDpsQ_@#dk#kl8U9pt+i!Hk3SYV9@M|| zF=V~}VQ|^ZnMw6Tg^E{Rs8M~$jTlHOtBaW8y(c$0sc4ygDe4z=fD-vVD%Z_NQvcg> zVA8TNt~u#!P`;hW9k(Ju@tx$pq?XC8CmW+~bW88XrGk1&<~@HaI-H6kX%A4oFHRNw z*b%*ko|7oipWjT^Q=JkL#|8swa@+PsW+6Yupjtr`7G>c)c=rqWfj6l(#mpK}u{ZrB z(5I?0$U>~B-C@uh_mUJ9#)SWI@)1=YB0p%hj(G6R!w%mwf#%u`=#Ov}FRJuClFtJX zkb)O+6~@gAJrpKf9r?{Ui_@~sX%SKtE-L!@+0aNIoy#9jX5J)KZOt3BD~#EE zFf|KZNr^@(kuv?AO7G#?_$#WlVaIn7mq_d_1f#junZ`WnRyFsyY>sbD_BO zP`QWbJS-%aP>A5hdx7+wxcKI7`nb zHw}71S1p$=*|u>W)_1|?>D{4<{o2*g+Q!2OQ_nm_TA<>jWCY|i(t+?9Pybnq-`{RTd~n@Zo-i23@w>pIV)PkZ~w=u-tX zcGqU;`GpUDlupv$-p5Y&T%G~qOI5#Kzu8~d;KOGT!wswg6*2Aa=;Jp5s=!o?c5wM= z?=KAXgoX-m?rq;PhX8+oSv0PJX&(gbW2&Df_~%rJ_0O3K-h=;nc`z~4pM?Ye8;ziM zh*}x|o2;Lv0=xGBb8rc@zb+pHwPn`Q&|QB2hPwNKN9+SdDs^E1h-_nAVfjlwfL2h6 z|LYm2uKnBP|LKuHXm%h^(1V~AmV1JKEhW57Cn=5!9OIY)e+ZzAcf9+!6(AOfR)I$z zp<>Jfs{q{d&y_jK@o$&=69W_dG&zgOQHQ@|`glqj6#@YFV5AGpR?Em@qb=bb)A>_! zb(PxiyE<^RjLb{+)~yfnK95at2UuD`qNei_pNheA`lVO&OCc&Bn9j4+%s4AVhRmff zS-S_G$EeiOxf!U`G65FrQ7Y;+BMAOkDdmnFh{_N+1h~El&)}v8>je_e{_RS)1q49Y zAuy0`Fg}eR#*#hdf&{fKCSPg>QPT(&uYe>1z&)TEyrrP<9u7agSY9B(&?HY?vpBVb z6b6nwU8EidqGma)FJR0M1GO>W2FAfddoC~!a{(&TO0Pb@R-|++rV)H>$`F9WG+FQ& ztlR?98#Do^ma7#w3E28-a{Z*`fVuqJebJ0*Ly2Wxl^>R04`OzyLs|$4~JBRh$$1 zME#j#pOWSM`)(#p0(sRxL*mrLRR2X6KrQ@l-v0aC_u4hZr%Mw#+)2>U8D_PfkyqSN z0HFdL2OM06IG)XMrZ8I1%&|a%oieA@*=Wz5z#Bc?8 zvvJS>P;{{|h+~6sf|jbM``Q7&dI3|Ryn6Zftl!udFysd723BmeDVF0*U2-8sFl^DY z`IT_3u}IVLcAfY~A13&M!k``g3K{xbJGEz$ZJ$RXuRR1xVqU=3 zp}%~vyWqDVrc%dujCbbyq9IVlF#Cj00X9+_H=6l=FvhOMp9A7EY98+;P`P`01YZoK zbfgWD%&+U;Tf59V?>J_^0EycmkUXnbutGD1{vl{VHNpMeb+C`vAHvUzE5q15CHJpZ$E!KWwL`i<+=X4_ z3W!J%EPO~uN}t|KtHny(N}E56xw_ByAIlJW3ZJB@ZQa1vrKe1{RMka197)`%GllYI zZ3POpvXv0yOBb4aHsX;ts}mb)iwu+h!PcV-fYcO^%QlsQ(Ra;KNDo;Eb(+!>a}g_s zymz+i^h+0Gl~5KrKoD8%9rF|r2Vi*)Oa_-4+_030^mp4@jhp8O_V3X%|1vVb z+T^_1ax{XTcHB6bJ=aRVNKtR=7t<~Bn2JkxSmsL6oBWq$g2#UFZ*m|5!}t+uCTiO| z0&?-$7y6B5UJ?d-I=!!PZIq;R4iQ4G5t+!p zF}{1b-ABhc(@Q*<>>y{JMI0}|){{3an}}j;!@lndV~-8JSPH!0{KXl{JKbwl804*tKBme$w}I;v_ziQCykQE41v()2lik8#jHdk$AeNGJ10R&T@cxo=jvSd2jTC zoG>@5B%xEK4ZRom!R58bxLP~`s7w?SsZD5wCRq{39sVYM?UGUqCVGk7Muk=Suhk9}$V{bR|4C~BP_VrvJuMt&b zYkrrLVI><|*8FQ%PNu(mW1}0F^$Mla3{_k}vV6%=)P*O+_`+Jd_s(W*1`gorYtmMx z9jZSUr$>5;(=UXUefN)Vma~NEY!!+Sf!c}RSoETO)o5JK+zd|rv2BT`Qx3MG zzOPOUp0hT@HODQy~5oncVFOeN^UJ$Ubpv&>ms?V{Y0N&J!}MP5N7~0nD5DT;=tkYczF#Vtqdmm9 zn}@C(T{8HD&<(I_;S>a-{rwniEJ?rDBX?Zt#o zX={|(x@yEXa#Hv(|Ek(=tO8fXGEkoG@7g==S~UukRR@<3-uKNtB| zL%~JKK{pv5aCbZ&jn97QLZVPyY{e98I zY{YpVj%=5}Yr4vXwRr+Qv1ev=D#U4&{s^u#xm(Q!&L1_tGY+^&F87yujx6JPU0tW=-i6=_XA%=$ zlL*~)=_}rHQ(mJkD4t}Cq%xip3dDf(YqoJkLC6c=e0p4%H`lL;N@YA{d9)wcMtr<1 z6}%UZB?7T-gcCJJwk{Z0ynqqp0fTN@b5bJ~b_WBwqf1YW4};)M*Qo*ZmCgnN+WSb0 zOX{1ZQbQnea6Ci}P_410U`VZSIN;_6Xm7%S+ytc;86>$n-H%4CPbXOPHO_>aQ1HAr z>{#)XIbO;{=!zPE8jKoNRyrtKcVXAEI{qM)i$9RC?I&dMgajYTN`l2k!^|5PCEsF! z)agPpDBSC$bP`+617Kz1G(znF#y~&arZldsB)?u;YhEFeoeMz}B$~3ja+g|oy)#eW z=M2Q^06GoT_T^>lu zfKYk+w6|YnrDh@olK&061B(y-C+auo`tM1Ae@EnhfHP2I{bx?(&^^H^po;w8z&aIU z11k8R(gl~E{@YCd&1AT<6T!&Y7p!>R67s zWL`tbJ`1$iz19EeqeRb=30l3mnhKD9nhXKzOXJHP|9Z_cjG$P&eCE!Iv$%IG4Mo$y zzI*5kU!`e5#b^w+`CZW#C$vyzx&yo$Zs(uO5d}O%$$XLf+;|miDsGBopgyb?*YS?m zHA42HRKT);8Ftn=c|PmLdBAMUiaJq2IP?NIQa*i7(oH&gZ_{gScE>i#uS0p)4F|AHl)NZa#>-Nv*Ll}G`r$ldt&zYF;#C!FU z)@lrQ9ouv6#U9b+;yq2kR$XDnzld;(VS11ni>u8MAxL|b<+ z9Ch_#o$iH0k5F@_EL!l;llUW+Fc-T!j`vU<9Q0A$-qs)IZGN!tuTYv;U3uzrdvg$3B6O#IgwL_O*(1q+j;L z(xx5vwuA*GuD()Su`?{(E~&zpY&h4Ezjs*$=qjKf18BEY_@4OaRlP8E;AlR zIs#j|Fg4R3DRl|7ei~@HWyJ%LmA0}n(LZ#d@QbQ_>;(gclstDB&uoCpX^9M|`)y!l zl9EGip+z|E*B^D2oZQv>#LE$Q)r9CCP7}>s{V0S0#;hoR+s({@2|wFpDM{li(T=XKr4@!fnACdP&n7oA zP*c|%ger-KD9#8Ioox_!G9mgvkO!u04JPM8gEIj;BO{n#-QPqhbk5vz-R+nV{1HtVpW zF`VdKS=9M>7d7Z!draAsGL*YS&XwA%h zA2HP2(yqMXvz28F!Ll79xew=0>8A`?49h(pWFEhE1aPBEI&63`&Bt>^41CUJGTLBg zPFgbtSdc+a*Zf6HRfB@PU1}d-(3gx^@yJ}q#`L=Eus37gB9w}f1p}#KXHTZa>&AP! zkW>q~8dJl5PMWPAzU~Ki$=BgqDXEIzTqwddWNv-<>fHJ}@b0KH)2LWojiN_6r34D) zmsVn;CL*B{fO1-+MoLB;hwjETN{?*b8KAuB-eCs1NxoU>EPe63q4T=wD-cUNJ_tIo zUO#j32g(gFtJ|wt+0RxlXF0qE?BrE`jLFPM$Aq>4t)eGWO5Twrr{9Vb=4?^*U;C?^ zf`f=|oZ&+J7_I(Z!l9ykv*5yfkWMh!1C5j2<2iLz7%{QFVqWdGC2&p0+m{{OIhQp9 zZZ7@dnqI+x;5iID%Xf~r*MY85_6(N3tTr~o#ziIk17z>=JVLb8@ppse6RlHgxhc3hR@2%dv=baXDQNv7$GhK$QNCHcDMt z*=940LL^e-+IVvbsu8=;`~$rjB8Lw_eSHSisx>ot`T)7`9wX-}7bY z*T0J)P#G>RS~IN=tCiOEPooD6&q1ZaHY$>v7VZF18uk+?)RXnMgST1LjYlT>E;kbr z#htwaxvzkoYl@0|DvN-=34(k9Bkkl5%>hMD56KbfjlCuw{1kr+@7(qVKg^SK&=eua@%{coC={|+ zl3^Jmh2i;909%d*I%a%MBREs$yDD#pfjc=hKKpE#_PHkC_y7{mQQ5rT?k3N$jul`5 zmc5Tg8hxU{*Jw4+%VV}4zF?~PL>$&fzF)pDCptGqJPT-9Pj1)Vh=bW+s91h2;uP5vn>!;Thjg@_08Bgg{4ajY4N5B8j%E@pM&6l=X-AC^`WcX>jno#Usr!_r7M zl3ZPB`?7w%s^#6#d3?X0Dfv4`Hqn72iuGj~0}Ebm(0jGXQZEW~e}IxYQvAGM5aU-r zYohVd(Fm~XG0uICM5Xwgs3z;BaLv?^cmhyM1GsJJ8uUbcUG_qt z_lpeeWS0<<2X`6$l8riH5(F+PH!0k_#>LeZB-4(qwp`26g}2L2J!|LMXUDhlmCswO zq$AOsuQet2%En~}X}~+{XZ`F5e!+UPP8FAiMeWA-ndz`>Uv%(9eW!6@|IE2|e8gsh zLHwAFo9BJ9CJ?q{uv`R~X@91B$G zD?hHO{0gu4L^4LVh5f{A8kk^K`4g@9M(Ohc_~OUeUKamf>eGGU_=| zu<+13iz!>ZTg<1@U?lf##ccIuM3}Q7II_9w$FJ(0vo+>AJf$F?2ve?pj!Zu(GmmU# zv=mH=YhD;-P)dG$1Y9)x7OQNochCh1rR5lo&E@Mltrt45=xATtpweP*e{GtfB-e4* zi}hhgFqvt7Rk9f>;I!#Igf;!F`DI~-FH%M^J3D0g{f#6o9JYMYqoD;TUiJ&MSbDr#4gC{m!)yp-$6N?7bHm2Mw$=8kbrag>P?W)p6a!P7h zHW;(g)rM1Zs}?&d;&zji%?(qEG+q?JdbnppZ>`^X$5_$*rSLnY72MyOJ1_esk$c!w znI|4SdcTb#*8fwif2(L4C}~rT4T#3bJo1`uHV?nCeQe6C+beNxwtw^(PpHnX=ZBn; z>ZrWIxl^CQD`)nnNg4z+CP!ZCUER_Ro|x3@4TkkDP1r+npr^h5YIjymspXQhVkFjd zDg4k&z2!TY3;EOR3p7Z(kA%evH4@+!q2#$Zt0)&}Zx`t@IkxH9|BtTgXh06|vi(*9 zYD%N;CJjJ4wFO{bK>wE=*f*2-|FmfXJsi|3-a*s1-9E0RB0(U>&`nFN&I1Qps6PYq zw(WPB1Bmy3)>sZCYYw{4-)mB50u~eKgxC2~%5N%XvZ4V}EYyyD>}kBqMN2*v2~_1I z@0ornK9#+L?r4f5c&tnHq}8NV1zUyLO4J*maDP4sDBJ@@V^}%B0w|Sk=+b@#nOcwk z>yW|0{=EpC1N&R{R(Uw0~#? zpz_Z_K!9@KWCXBMKo9?ZUH+ge`{0Zb@cH228n8wO=PLGZ+&_dzy^(qliK^WPg8)_z z2JKH(c@{VXKz$6LK;9Rr!T{|bCfXmre{hC6e*ZAy{&)az|2=5`dDKDwfl#WbK?er{ z|Dpx-9|+zj;D1N3Dx7+T;9w>HuZa8Cp#Mr7zyXy2kZCjsNR$E}4wU|np#BqaK=$lE zNfSRX@Vl{3Skq9mXXwnHd;;*M1vXy)<QZ!eg=!y zk_!RPs3#5q&ZNztH?sPFCuT2PU6WG@3rOm%a?#tt*&yyPr^}89@$^p4Z7@oQ<~emz z@X(rAXO>CFtU&RI9(rot4g@TcbPa%)!Y|oQN~O*>rmV>+5F4w;cUS9t(@?ODL(U?m zYTJRS4~A2B=n0hO@lp*tTq*U`7*$N*rl)RFq!68l@nPqT*J*LmHejC3Rm%B4i|y)t zSo-x(n=SJ3$MaYAQqPfLS4nF~#PIW;r52-o5S7g28I5!n3T zUUzutOwHgLrE;-iiDcX_{n!gN(BwU2DEv(KG%kw-1Ku!kI05T#-rD^oUX+(oHyR!m z=~t6jPT}L8)810htSYjVoz|PrA|;DfjJPSv+3mjY7ZK0F?RWu&Z1uUSIHa5~+9H@J zpYbGpK8skPjPCa-jUKljJSL3J$L(sxZ_LdQFIu$oOxmZEHsZA2DiY>c z?W*f3eZyLB5#kOr#Kz>+@h3wz-(MIqNU`QkRK>+pVpi7Oww@_h1C```?G9zdN2v{c zkgnBwZrqCP*Hm}T#4h$Z`dkh4xxyN0Z=suMEw^sRQgK1GXRz8tO&ga(3e;A*k<|Lr zV5n(AQ_MDe$F?RK=QFcvXg69|R)~B#Uzj;{(a$etV^tt92(64V?~HUP9l)7?Mz7CT zE(+)lt05RT9@ZM2|MKFopZJm$GOohduQ;n;Gp=doE_Od$rS3Av~JRW~`_IZ9xKyC(N^=x2}97+DnPMw1c zaISKt^<_HH9nRKnu7gk%p{ZmU$kj9|I<Rj|B>P8mr7tJQF+YdxdppW;|P6VA#-M zNFl$`8__S)PCY>l?ASN2KncpR;{B?V^=rlorfc>)j-s(u>JMgp6CS&KV2Yli{pH8C zvH3wOi#_PG{gsKGuY2PS2o2l$f)*P-Lf#mIw97u@n6;oN;u}~87iflrIN0Gq7DsbLGrG7iHeYxjc z)U0#N$njDC^(x?uG+=(<2oa1I>mkuzx}`?f;)+U+MOnjFtHVHmvl9_;_ay!$ea!-lu}DIc9oIGHJ}Y)-o9T@Y74Yh zGmr=l;=Np}w|^^I*CwR-h%3QP&s!e$?Q$`9ez70l0igy)_%vP)KuqQKlZ1L0GPs_h zCyjU(<|Fah@5fm_j|>ubWv|Ca2LKl6sa^Y+IjYqx)2rKL+q3+j?zf+gf67nqq{BP+ zhkVGUZUPuD`CNVC$^tVpelKY_+fYUy8T`84Omt<^O&Um1Yt}2EY2|_ZH`T+Yo&l$8 zK{hU|S+HI;hFZKfyud#nBp9XzsZB*)dRd(8>5y|rDbG-*m*Eo7%Y_E;1n{;BGbs_^ z6(|92TN^enx|~8@lrs10j>3I;l0&t0UXnE%E=dCRIYyif!1LDCd1I>qqJGOocMNdc z_u<4GG_u#R)Sm_~3>KuuFkrkRovZ4KeqGeYt&YuJkuIOtSY&^TTP(=1R6YTj7FL(w z)-pGDnm$)h;$qaoec5T9q%vN8c)o4A#Yu#YEz3jhVNKZ$MoUVjHERZL)z=wAUNNNH zK( z95G9(#eE~EcUNzM3}d#+2Vgw`Z~)V>4!Ri-KLkcQGu3K0r)WHubo3F*vb|)2o|`n# z%b+Bux!COvCovg3m+3u=ksBl<1h*ratBgb}i10RIT4b#1G1rQ2RusbvA-n@gBa zC`45=pFkawE{TpebB;#cm6V9Snif2%#S-tMc3k6!nbvNS9aCM+;Po*fKH-49@f`~wXVhg6>nW6yG9V4ddp)wfR8C6&*XZp2n4gr3aG7gijz zwS)eQjyHuPvvGTYuR}A)XXe9aq$^7IL{~F~DaoA(;4KICS_;_5a}2Lvb2hqg^|*_* zLQNMRF3ZVl&-j_9!j4H2?!-+p+70U{ljH*DI_qu`PEp;mb<3Y3ItYI|m&RAkfEtFu ze^p9mf-K%Wp}zis`o-nQ2dsXcZ*f~A--AUC0j16tK(qmOn^FT_#w4B(0!h2&8ES~O z9HpMcy?nz^_Gt0z(+6Q4wz9VxSFC6y{wco##bTiA_NIt9 zZZ8V9c2^Iu8ds)k$x{j2XjAhsQoBXaAZv~|dFq>z1Ch9+xO10UUJhpZU;1v*5vB`W#wHx!AnU) z>s`XMqyJZPUmgy17X~_Gmpxgs=4(;Ol09otDp^{DD5D5v8M2#UrcX#hJ|WqXN@XWI zV<+o~Y*}XP$uc9$48zPlL-l?4-sj#w?q9d(@jU#_Iqx~|d*1z=^G2zppAv)Vg`rEu zDL|7>da)BrX0=Q#wpt0>W4#(s>Hvz!8a)*a{CTOfo{c3ZfncdOF^s&uY`5@Yp(#NR zjccc)ItiDd-?(FV8asPNs@B(@9r%Nw>*eYxeDXaEFk|x&zpap-lz|%LseJBjy{Bik4FR|=5Vr{D7$`yR?yu!J5eq033}Nhc zVtI9fk?92}4gxf?NGIg^-JMh$Cj;l{jQ!|BP7t7hAJ77Rr?&ui?(z}n-TibVfblcy z1|xPHTqfaNSoEs_;P%JPJoHWnYL7ep0cC_fgN8aOGhU8|GIW3|yg(O5IzZYK>F5Hj z^s5C#ggeuNUOUs{W2Oz|ORn9F5&ykhun=B+Z<}Gfu1q20e^(8L{s$m9aQ?!1w4EM% zu+V>=Up#$-e((2?O~f5!3Rk6DHr*8$NAc5?_&%pe%k`JbzC!!sJ+_ z!=E`?Owm~cog*=Ux27Bg;D~}}IYNh5JU88(YpJkz(>6kBriTi)`tz}6Tfozmz1x2F z%;>TH%Cy$NYXKjXRazI!3@v&qu$qr+_E;q{{U);?6SKPE6lImBRz1juL)hy&%wvxA z2R;J_3f$`z<{fwEaVMrocMBHDm1BEw2rJ*BNsWQOYBxY7E*qm^RFXJsWM*++bkW_( z^-Sj}()E&2GOhU2Qmo793z(9L{C}5!2`%%m>Izcw+)w`!BONcllawuk2FySYsWQZA zth4`C)(rykJ81#q#)tpd;;Sg`yp}rIK18~tGFRNxj>pa;G7CDudh*i%i)FaL0$ggg}UEuVU}q$41hEr^3PjB!@GF`)>%GBzM1+agY~ zUi&l#=fxL?f_H{@bbTLJ-A`J$uc1sjp#L9|HZuFAH!zoWP>xnta{Wa=oH0XU;#~aelMeMT|WGYl$`Y{?mv8GhMomgMPf3? z@$OScUubg8yF8j1ic7i!-VP^NMfpozgpaQ4NPB3*vWvVQ+aZ?r^*XWq-dnN@8?n!w z=szeUJ3QKf8U7*c9L9@(a_L?g!FfuBd(k2I#;fkstVr#&G^h6QVv z3T$w8_l8}-<})NKCG~89E390(Hj}Fg{I)=}+6sS*GXoM$nWONGQxIOhqU%Dfj=mTT z{9s-ff#0c<>sal|z2OAwk`$h%+#^+Oj?tLz!lko@)+05(P-g%|k*=(u49HL>zjVuY zKhC9O-%u5rT8s~zd5wTyfaV7H;YwXl#RllvBCG@t@5=phYm<|7wvssI(5s@|{3vbj zFMY#+OlA=CPl^99RpRP1MHy0+{HknnCz=rEf+zwWZKO+LCD&dX<7dU!7HA$$sLq~^ z?lim%cUr2?6Fyr2kA32En)=Df1b25#nC1DVQUN3CS*U!tnsw>56VyiJkP&4>tN`bn zD~0>0Rq>h)oC9!zrJuinwurz@Tx~Bbfc`vgt=hI^@Qs@Bdb!Xl8qDe;2Z!5=fjFW~oz=z5daHgYj$FCOW%62{#kkbru2w(vvDgGDy3 z$Abp~bYAg!Ro@;Z-z3Dwd)sZwX%>TLE+^8|Ins5HXI7NpjMwzD2PWR9h^nmBElPT6 zMfa@h?!_A{l_XmLHBMK`l8HlRU*tl9^$ULzzJ-R!-i-`b73_V$p)Up#V-P_2gGaOa z*}{1xVpuO9ty`gzEq3A9xiLc7cHv)`#Qwk+{gD&%{)K4+PN^4bpb&|?NJ_{ zwZ+=pAGVWaX$M0G-N33aMY*H&2UI{=N6zRJKt&O=X(3cLD1@nO%I3SO-ZAMeA}rE; zb6*rKk)GP4yG8()za=g#T(3zMiH0@Q-JXT8xp$fQ>Zg?!-kZ>oN|e#aFKTQeHK=^FqyET?4CNmiRO>qLPbC9oS*KoYvt!t{0a@wx^D<~H_szdUQ);QY?uD!#fM)s`l{McrzX@`jhzY;Zepe@1hQ`={1- z{i60>GMqtN*bTuUt&b4Yo;9oz+`fI}o_g;6(E0NqNaaYUA=k3p-CyS?aQ&0$S$sGr zj1S)~-XmE)aG;0fw;3)2Ee>fa~h#xhU3L zYuf#%`ef2eq7KwvIh{2&qn=!O@2`)8E84wSCKmoaY+IIZ!tQINCUN%36G9hZkEuIgP+$D zDLR$V?zJIMfb4+Dm;mA)9OSq(9Du@2snkld)(+dwmZEiezpv?Uc@;#DrRZtLgE!B# zuv;vJ)fH(YScHX_o7W#d)6MgL;k%*J+MkxgH<%;Ew5iTac3_2a>z+nuG@5%VBRPU$ z?ECJXGC94b${Ex+tGJTZ3YPRj)(*~WSmkBolvNer@MP$hteL;DG=b> zUk`i)1xsYWvHb?u!#?Oec}Kh4fTxrBMAkXWJieV|Ke(i&Opspx)|hB8sZe~2&el@u z3F5<(QUjgRX5UwouHu-i06X-=O&D2iG=R zx(%Ah6%MpmqVw!GNcYi}=2QAfdBrD4e?47teg=eqPLj+=i&1=O)lYqqnWjWbwcY|3 zi94fBCuAcGKYxHLTf{m-57c^(bVMyX|C8~7Wdj*2ODVZseiJ?aMuU2bt6-a#8*}`_ z-u%I9r!?#71yBi_fs!v_zGl9?=^9beTx$E51ZKMo+u=SKa2eUg(I4^F0jH*zjvD$0 zwPOT$@z(|~&+YZ7e=^Pgd@VQj#v>^aUf?jhuebfYKd2jB>-NuE>5Grr*pnt4jS2ZV zOl0elobwkSOWdV|{{+|8pP@j(rMCtb(J+T$M1nna68%Evta9_p^O?Bk=Lv+2ht3zqB(;Pl`?t>dfopHkO@jKN(*pBJ<5K(Sz8 z$fA!IIIiDtLLJ6Z1nW*PlkFlcXY2jaZ*PBaCyvhY7+gk|HW1aiWIUc$H^D{X!#4!Z zFB&YFoN7KL_C)SzPn*`m=! zFFW4M4dwIP_6jq(M-SBKLD~*~J6;8~G3Gn^%&78vma}xSEM@G~On&kq!qxU18*>7V zY18JF(Doau@}N#zQAk;HxZ2Viztlh*5~jfrjl0`IQygXTOBU}fU6`pe(HI+aPB9Jm z*;D@Bh@nk!H#}n_IsyOJWkZcC4t3Zs4sCETHvQ#az4kb7gp|LRogItQ+9x7)Q?eauKD zZ_qEfg3{KLugQOX#f}MwEIEAqx+2wLVZ~H7r(GBKJ@6i1^@nc4MfzI=VZzNSx|YUVs%qgX8~d~~D&qoSYU$S& z=CVWObVb~aHawd>RT#EP?vA#)=q8CXsiTUoSreS)+6FY(e1#{!xY@uroxB9!P^Q(rD7j%U4gVsL2O8_tkaW?yTz2z&gx@nRiL7E4D zDP8;D)ev5jmBTIfdy1D9shs}#09p+$1^L700TAwrCA0SG`2)4Xx5+j{S%&N{oFH7? zG2=wC&A!s8vJ_z2=w)pOz{5KeL;VEHUD4V}3jjwJf26fN>@-d+KPQXP-n*OkyIML; zBLh4;te*M4E!As`R(S~p=1y%#|8&PWDkdUkS5D3PUYHr9j8JN3E~e&T1Nmz6y1$3n zj#|Fosxw1~ltmMd_00Ja#s*JAa97`NJ+zK185Y78)~Cbl4xz=B8Q1Uf)-nh)Kk&Vr8(`7CdiYR-w6XmA2NY*TXRDsm5XYJ%zw?nTsiSVC5T9<&!Gtpq4qGs? z%p4+rTx?5joOMs$*0*%t8j^#=HF(@Jqnd;2C<<@%+>IHE%uRDyMFRHGnr%o+^@*ZU z;;{iGru<-sS_~(HCCJAJ)*^n~gnQ2DG$i6H_s9pnzfSHN%T!EQ%;x*YbRTz+dtcYo zdZm(wc#uwMI3zW#l0Y24p>||=8|hFf^>@OQo6P@&NbW~@p=8zvSE>q3F1VR8c%wpw zbuAB{t^9G;fo+Der>|!7d(58Y|Bi7NoilV82N^_#p{~~~dAdc!o=oj9+AW{Mv_+F3 zo6gr$Iz8Q-I3c9*1B$En#PW*{tLP6!QSR{?Rg&t+Z>#-r%E=%sMkn)WKrldN7S>L5;QobJfcZi6HX|4JODIJVbaVN_W?Oqk$ri=} zJ%MIcU=Viym&)j1hkSX)9a^!7DF+Ac-`O$+GV3=T_FcfdwfI;2|Mc0Jg5C#!s@s7{ zc>}1-W=lIA{EO@yR?7dG0%|19oPW~cK3DgC*5T$C$2JLBpubTKWQ-tx#_rwJ!kU`! zClwk-q?d(;*<%asSmnTrQ)5&e+Bz|QsuWSqZ!b}>$a}#k7;`?vVUt-QrQwnz6vDky zP`EQMlLE++3f3M#^VNc^AAQdazyawlaiE-#0|4PqW!GAYq_;Y}#_SmfS8hK&*(`hw zRbNm5hI#owlDJb`AP61$SG$*zfLNc(;W;{-AQzi0IL-!Jw}VEB69g>z6|=0!X@G;%nFtnMfD`}^r8fgO z^o{=i2@K*mr~&Zm93ulbp8sEv14s$h(y`^O73890OYcX=7W9LF{&XVHn|~AGe*?3} zmQDn)ktfj5|2tKF6OB$RI?({Jc8CT>p#}H-S_H$OB=i4)s=vqoZvHj?p{kFgoR4Ut zjj8Ak@n4OexThFIYEo+NW!c#zr}NRue6=J*#V!#R@7wW1tGarbcN;B=IOMB8mB&C* z=jRx9(EL<1%GF5FPd&0jgwvIk?;k?OX724adi^au&qVyoly7$~ewRGH8dJZ3UC=3a zf+8Vfb91{q%JHC36p^eLsE68y*cyt3FLEz<)?`iJG>p* z0JbJ@FP$=x!LGGB>~!hnuoZJ*T)()(G~6Xn%AuVNWgVt>1ki$mUGq&D=$l(tr#dJN zr85)WR&*hD2Q~M6E-;kNSi3`CF*7Ud{Kj-smf_^oGBw41`A=tNcMWG&)exTS&9Ht(zlCvw4BxrKzAncN=2Ip1AJaxPy(H1x z$cg_h^f^GGFm+MVd76T>$g7?)y{IkqgvJB*!?6Qppx+&S=HZ&ynB)IS!mGy`J zPFpj%lfQnCXC(JT-^%&31`QoFwq_y2`cJ>wM^6~n7aI6#Cx>R36L0}3J9N7>TJNWm z(xHfA`tj|hmW8#p0ASj8rAMErveR{)NB*cyg`%^>>iJi>+O3|&E!S?RF`b4)E1^+Z zU!YZmzQb}m1hNLQ0>p#O7W@T_oIC$aqmu51Jb=YS^kX0+_{D3pD0!PNcjlBu7nm}f zZeb;Udsq;*v%bZC20w@Kw?!!{O1ep+rf{v>x-cQ^oIuH`gPxUkKSpTT9o~*+L*^C_ zV8au05sqgf7dt<@!VP}<<%Q4I?-}5U&eJSMIzJfQP$|KE-j83f8%x!i3Gq%1TRn;G z(2}2<{Fdwj52Xu7l)5cJ*bzc&I?GDLcW<`ra4(o&9dYAvp629G8X~oa)X>$+jE`-- z3(NMT`5qoXROsPD+>syaw?j6Ndne=c6st`hudqHS9gg#VdRv?F(FJ81*PTm}z}JaH zpBizvlA?o`BiK)CN5&#{zCgYxK*l!ubR$Ec=}P}`^ZHJI#^0< zOR*+9jxljq9;mL<6c%DzpCs96ZzcA5`}B}X32i~sVx;gE*iH9*?F}gMKzL(fV!Pw? zUn{o~+Xd}|B(9eemg>_Tov#TAaBKG_q4ux6uCFV~=ILj{E07Kkq@`&^C6>386x2^z z#g!wigQ3aQoz(&oEB|!ZaeHWGHm6ndVm69~hZPeD^aUUQ@7j#T#Pm>|-mI>qVq{Uz ztFz6HYj@+fwN|2hahpjyoWI!0foIJc338j!7|v?EQ(QCB>(_;{Y|wAT(kd=*2j&Tx zvlCn(x^ZU%Y~V8*eN@=MXzR;TRpDCHEp@Wu^o{E!{yve9jZrr9PdDsB=12lY48myD zJ^0H0yt$em^HO>dSN!3O1)Cp^5C}uQ7|Ow6X1WJE9Hxs_8yGYwPB_Wq;HInVJG4+= zN5VN^>-#k|lj1N~a_@Bpa_;d`YxDO(=9Q_2(5n2p*<%-Acec;Q=^~29l&UpXP@F)~ z+wNg<*TRNfb0wsNKnK5fgcy`jiq!I`* z5cRsc^+r~A3l3;M+aU?e`IpcA*k5Ah`?`U%7$G|>9vLA z$(ZT(rdTXy`xPg|Z~!D~W5P#ie0LP1MhG|`>c+rWVDtXztw#_O)b_)rB~va?lD#Go z67~AM`jy*bE2*+Dx--UI5{-m7!rZFf=)JLdUP)cAPbX-O_0}%9Nuh&@J+os|>v0&L zqsjZ1+zz1CcC?JEc%|s%Mrq8^fbq~va^*Q!7y3Ey3HpB1ZCFmOJQ}5y6W^gBA-D_N z1;RR)*q_o$=rvX5CX9EYK3{)cjt!B`&P;@K?TMaTm@3e|s=hk^)S-_Zk2xZ-(yuqo zgWp>VRM!z{pGw2GJ-E=zBLPhUwCXQ|uFB}%p`PBA^T7?G7;9-XPBpc-u@v^qfY6pB zLpU|Lq}0|rAhT_wRg<-&h$X=sEO6DwGSxU0L~~hQmmLqG^wzV8?*@#1xA#`{+2gVh zw^JZCJ&kslQg>oyBnQKw4twpmgeziVR4^y%T3m|3*&U$k5`>I z6?p1F4xyQYfHl@Vw;5r~LR`{YQ|gZEr{-!v#Y&dV(y12vIy<#y5KQ@jd(o`|Ew=D- z%dIc=z-kk!f3{IrB19#yEbciiA`|LCRncLli8kJgA8~p6Dz7v6!9u*@@^OZ!+#SqK zO(3i!0dCBjGMf%v0q?D-S7p7W52c_ zjB2-Ic6d+_(o8DueY3aY%$<`+;r`~Sx+VRkEM-XYLVrp~#fkz6> zZSp=;>mca`$zEcmq{DOHOtT1TeFyUfkJ{e8wI>#V=dP@#6OUnKS=|fBMlz+8@8z>@ z7ss*?w$l`gppw49oV~S(d2ONT2OCcpF1-b&&OI~rS6R=xulZn5F~mf2MO&q$T#bFqvXdc%|$t^lj-gO+e83 z0x~+M_MXEwl_HKuwW=ux=*BVlxNqZeIhgxvsUA#39_zrUCruObs-pEk`A^!40z;wD z)OdTd!U4|Dy2xmeW3%07fqd-=Cs_25i0dR0O)o3@gZbxBo`#B_FJ)?w>WpQ1d4#p) zwHl%Rf~}cpwQBM|3%Ku;cfiX2^099~9QDZQhy383R4wXu(WJjoB-0}CJCP@MFraDa z2m_h5Zz<36MD8AfepK5RC84c;X339Ul!si5(f7_oG}o!xrZK&W3`roFtg1>MM8mkp zy&GsFfwxGxEDYr3j!6nJ!tsOJI;}_EX+l?Y_Y(Pb)g?D_<^cgtNq@MXp?XlVlb9pRpc;x4vD3iaBf;x7Wmn=7E;o zI0d+BJ7x|bV{64sluZ$=y`@=!N#?;{AIRNI{?l(KVqsPQO93pABKzbutA;EZwY?6b z{r#tp!fmM7`SFnt!nKzOH#b=~AT~k0Kdqbh*P7n{mtc|7=9%{)vy;AmGmt%i(EdWM z7a}JZacTH4x}~zy&luMGI|!>)=E3`k5_pD(6cmVqQMtLUzG-}{cKfP%VL>je%tX5L z`3pn_5SH&9WL27($o?VS{sB@BZ)Ncj%+g9jLqkg`0+L|(HX6iyzjjWa(a*`r$pN_` zJ~uRq9I+V~t2Hj4^F3( zSc!w*igEW;?Hj5{Z8Eg7xZ%T&kM!BF$&VpH6hc*C!#87RnpXDV)G+5;rSNRRum>mY zE{7u;+xMF;r;|h_B!cYxOyNa*h`t|v2rDkU5u2owD!S6g%QeE_^=!Q>g+f8SjY9`v z2N9T)>KbHwe(09qA9{H`KKf*-(x>kuRz)N)iG_uFrwS@I!yh3O?VztNnWmqiuX~|u zDQM6wYRdTgKuY_Gj*VQGrcCyx0+m@_#gW4-&x^zOerp~DRWyO(6YkW@Lp8w zz4cd)ZIZ|}*rqT30Y%Z(mM1X9aU}A(%iM;ZSUv6X;)8z|^S+>r%3?s3$EK~GkqAU9#U7iFhR9EF+*|XnZI5&4nL#P(!ddN>6Gd^J*t5k*h7v9UY2Mh0FME~+x8Y$1yULUlB zk>b^==Y&iUllO}1lLLNys1b%yH@!?W=H94HNmrk^sb*AmYtUfOe9&QVz31&En5dpA ze-o|pDlxTaC7E>C(GT@+WNgLP^;04xb;vkaEpk-ZYcfsJvz}dg>nd^xb<=Oeqeg7z zxt(dEIBdPyia$5jMt4Nj<)DiU*Oc-TN(ylSNP>P2H5#?CEbKdVnRJc8i$0Vf=R7QOG}dtc5PKRH?ITro!%Gh zYS-~q1oFZq&k#0qQx@5)R3&BEi|0BcTCa-fV)tOfMxbdu(AX&oN=@kK7-x>4oro?Ex8@xCNp8 z%|@+5^Lf&ccZOS`5H6{gm`|se?u%fgZ1g|Bt3&)8i4|(Rf-u}Fx;8m1i(IMlrybKA zcN{NtPjMgRZ5}VVEqdEp>P;^%cT7quU+Qr_daixZHH__dDtn$CebMW4dulIfGnBO9 zfZghUg5xce?$n;^lycMmE~h^JsCNQKv)ZDRR{2x6&fT7X)u?QpxG9K5ai2^Mb$>kT zZmyopmn%>=7v}JmFWLXpCkV(GL(m6BYT=!%#MYKHWIz@)S(ihs%hGg)@{n$6Kmk}> zmWV&Q*X9_G_q~*v%sF@J;ch9P-it66P+TE53By!jv-M8$iNI^EM*qiN@(%?;78R1C z{xOEFE0w(~%$;p;fz|e{WKdk%7VGW1IErYOT$t;C#tP!igawFNLH4hxO1N$#LElkd zeOu(sLfD<@8ir%rZ5h0!L#-%*tD7nphQ4>AHa*v~;UNjyk}8R&QNVqo^0T5?s+MTQoU~?UU6OXMroijnmD#>z_1L^R=OPSl#5cT&Nlsh+ zx!U1Y<|^kCoqD`33#11D$YOO9#`2g`_m3Oo0+t5kyl1+D76cpx07}pk%Agm}Kj055 z_cqx8Bl*AeQIzI;wZ$-QAemdD3r0R== zkk0d@ETNo!_v%9T}X*@0=YxJ?Q}E!1wp*^K_94Y{9huum3(p a(0JzV9{BsI0cXLI#z|Nb9eB+RM+ literal 0 HcmV?d00001 diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index c3fbf4d6..bc9426c6 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -75,7 +75,7 @@ We multiply by square-root of volume on each side of the tensor conductivity to .. math:: \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\ - \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM} + \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv} Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above. We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like: @@ -114,7 +114,7 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of .. math:: - \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)} + \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)} .. note:: diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 7903cdc4..f893097e 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -23,7 +23,7 @@ Solver Utilities :members: :undoc-members: -LRM Utilities +Curv Utilities ============= .. automodule:: SimPEG.Utils.lrmutils From e9fea3bad2c49e7da15e1474f6864f20e55ab7c8 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 14 May 2015 23:43:32 -0700 Subject: [PATCH 11/49] lrm --> curv --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/{lrmutils.py => curvutils.py} | 0 docs/api_Utils.rst | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename SimPEG/Utils/{lrmutils.py => curvutils.py} (100%) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 6b04e042..6637a138 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,7 +1,7 @@ from matutils import * from codeutils import * from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile -from lrmutils import volTetra, faceInfo, indexCube +from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate from CounterUtils import * diff --git a/SimPEG/Utils/lrmutils.py b/SimPEG/Utils/curvutils.py similarity index 100% rename from SimPEG/Utils/lrmutils.py rename to SimPEG/Utils/curvutils.py diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index f893097e..042aef67 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -26,7 +26,7 @@ Solver Utilities Curv Utilities ============= -.. automodule:: SimPEG.Utils.lrmutils +.. automodule:: SimPEG.Utils.curvutils :members: :undoc-members: From 8a36cbab3b0c04b14f9b50a44fbd87387ee73211 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 15 May 2015 01:18:25 -0700 Subject: [PATCH 12/49] DCfwd example: - TensorMesh - CurvilinearMesh --- SimPEG/Examples/DCfwd.py | 71 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 SimPEG/Examples/DCfwd.py diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/DCfwd.py new file mode 100644 index 00000000..49fbebcd --- /dev/null +++ b/SimPEG/Examples/DCfwd.py @@ -0,0 +1,71 @@ +from SimPEG import Mesh, Utils, np, SolverLU +import matplotlib.pyplot as plt +import matplotlib +from matplotlib.mlab import griddata + +## 2D DC forward modeling example with Tensor and Curvilinear Meshes + +# Step1: Generate Tensor and Curvilinear Mesh +sz = [40,40] +# Tensor Mesh +tM = Mesh.TensorMesh(sz) +# Curvilinear Mesh +rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + +# Step2: Direct Current (DC) operator +def DCfun(mesh, pts): + D = mesh.faceDiv + G = D.T + sigma = 1e-2*np.ones(mesh.nC) + Msigi = mesh.getFaceInnerProduct(1./sigma) + MsigI = Utils.sdInv(Msigi) + A = D*MsigI*G + A[-1,-1] /= mesh.vol[-1] # Remove null space + rhs = np.zeros(mesh.nC) + txind = Utils.meshutils.closestPoints(mesh, pts) + rhs[txind] = np.r_[1,-1] + return A, rhs + +pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5])) + +#Step3: Solve DC problem (LU solver) +AtM, rhstM = DCfun(tM, pts) +AinvtM = SolverLU(AtM) +phitM = AinvtM*rhstM + +ArM, rhsrM = DCfun(rM, pts) +AinvrM = SolverLU(ArM) +phirM = AinvrM*rhsrM + +#Step4: Making Figure +fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) +label = ["(a)", "(b)"] +opts = {} +vmin, vmax = phitM.min(), phitM.max() +dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True) + +#TODO: At the moment Curvilinear Mesh do not have plotimage + +Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F') +Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F') +PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear') +axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax) + +cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)") +cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)") + +tM.plotGrid(ax=axes[0], **opts) +axes[0].set_title('TensorMesh') +rM.plotGrid(ax=axes[1], **opts) +axes[1].set_title('CurvilinearMesh') +for i in range(2): + axes[i].set_xlim(0.025, 0.975) + axes[i].set_ylim(0.025, 0.975) + axes[i].text(0., 1.0, label[i], fontsize=20) + if i==0: + axes[i].set_ylabel("y") + else: + axes[i].set_ylabel(" ") + axes[i].set_xlabel("x") + +plt.show() From 2be920800a57813313d3848328810d8c0bed6ec8 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 May 2015 09:02:26 -0700 Subject: [PATCH 13/49] LogicallyRect --> Curvilinear --- SimPEG/Examples/DCfwd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/DCfwd.py index 49fbebcd..385babf6 100644 --- a/SimPEG/Examples/DCfwd.py +++ b/SimPEG/Examples/DCfwd.py @@ -10,7 +10,7 @@ sz = [40,40] # Tensor Mesh tM = Mesh.TensorMesh(sz) # Curvilinear Mesh -rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) +rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) # Step2: Direct Current (DC) operator def DCfun(mesh, pts): From 6380add4ee59f246859d98c958de3e9c05b18e35 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 May 2015 09:08:05 -0700 Subject: [PATCH 14/49] Logically Rectangular --> Curvilinear --- docs/api_MeshCode.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 9fb52e53..a3813c13 100644 --- a/docs/api_MeshCode.rst +++ b/docs/api_MeshCode.rst @@ -18,8 +18,8 @@ Tree Mesh :undoc-members: -Logically Rectangular Mesh -========================== +Curvilinear Mesh +================ .. automodule:: SimPEG.Mesh.Curvilinear :show-inheritance: From d9d9d6828bf34c07c968218465458d6625aff391 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 May 2015 09:09:04 -0700 Subject: [PATCH 15/49] updated documentation inside of CurvilinearMesh.py --- SimPEG/Mesh/CurvilinearMesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Mesh/CurvilinearMesh.py b/SimPEG/Mesh/CurvilinearMesh.py index e4eeda3e..f8b0bcd2 100644 --- a/SimPEG/Mesh/CurvilinearMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -12,9 +12,9 @@ normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): """ - CurvilinearMesh is a mesh class that deals with logically rectangular meshes. + CurvilinearMesh is a mesh class that deals with curvilinear meshes. - Example of a logically rectangular mesh: + Example of a curvilinear mesh: .. plot:: :include-source: From d8b04b08592064f9491f2848d11b494333bb10b2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 May 2015 09:21:54 -0700 Subject: [PATCH 16/49] LOM --> Curv --- SimPEG/Mesh/View.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index b5ede21f..015b21aa 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -576,11 +576,11 @@ class CylView(object): def plotImage(self, *args, **kwargs): return self._plotCylTensorMesh('plotImage', *args, **kwargs) -class LomView(object): +class CurvView(object): """ - Provides viewing functions for LogicallyOrthogonalMesh + Provides viewing functions for CurvilinearMesh - This class is inherited by LogicallyOrthogonalMesh + This class is inherited by CurvilinearMesh """ def __init__(self): @@ -594,8 +594,8 @@ class LomView(object): :include-source: from SimPEG import Mesh, Utils - X, Y = Utils.exampleLomGird([3,3],'rotate') - M = Mesh.LogicallyOrthogonalMesh([X, Y]) + X, Y = Utils.exampleCurvGird([3,3],'rotate') + M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ From abe0a9affe60f0d12d69414524ca887af887c27a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 15 May 2015 11:35:06 -0700 Subject: [PATCH 17/49] =?UTF-8?q?Bump=20version:=200.1.1=20=E2=86=92=200.1?= =?UTF-8?q?.2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- SimPEG/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 189947df..e4b0ba03 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.1 +current_version = 0.1.2 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index ea08c8f0..61b91377 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.1' +__version__ = '0.1.2' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 161c663e..75319788 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.1' +version = '0.1.2' # The full version, including alpha/beta/rc tags. -release = '0.1.1' +release = '0.1.2' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index d1ad32e5..b1f7b742 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.1", + version = "0.1.2", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13' From 4a1772bbb109cefb2bf0a510ec65ed26ad69e8a8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 15 May 2015 11:46:57 -0700 Subject: [PATCH 18/49] updates to website. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index b1f7b742..f1e62dfc 100644 --- a/setup.py +++ b/setup.py @@ -44,7 +44,7 @@ setup( long_description = LONG_DESCRIPTION, license = "MIT", keywords = "geophysics inverse problem", - url = "http://simpeg.3ptscience.com/", + url = "http://simpeg.xyz/", download_url = "http://github.com/simpeg/simpeg", classifiers=CLASSIFIERS, platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], From 0c432e1d821aae11ac336a47a15e8d9b4c354368 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 15 May 2015 12:36:26 -0700 Subject: [PATCH 19/49] =?UTF-8?q?Bump=20version:=200.1.2=20=E2=86=92=200.1?= =?UTF-8?q?.3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- SimPEG/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index e4b0ba03..f739d080 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.2 +current_version = 0.1.3 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 61b91377..369134e7 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.2' +__version__ = '0.1.3' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 75319788..3bb33621 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.2' +version = '0.1.3' # The full version, including alpha/beta/rc tags. -release = '0.1.2' +release = '0.1.3' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index f1e62dfc..31c07dd5 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.2", + version = "0.1.3", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13' From 369694335afb5c44f88faeb817d2cef7bafc2e63 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 May 2015 13:47:33 -0700 Subject: [PATCH 20/49] return dobs on makeSyntheticData call --- SimPEG/Survey.py | 1 + 1 file changed, 1 insertion(+) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 29bdb452..8abb547a 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -358,3 +358,4 @@ class BaseSurvey(object): noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape) self.dobs = self.dtrue+noise self.std = self.dobs*0 + std + return self.dobs From ec7ed8a585f60c63a942e9d53f9416172ef7e114 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 May 2015 08:39:49 -0700 Subject: [PATCH 21/49] consistent size checking --- SimPEG/Problem.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 1787a472..e729ee10 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -160,7 +160,7 @@ class Fields(object): assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' func = getattr(self, func) out = func(self._fields[alias][:,ind], srcII) - if isinstance(out, np.ndarray) and (out.shape[0] == out.size or out.ndim == 1): + if out.shape[0] == out.size or out.ndim == 1: out = Utils.mkvc(out,2) return out @@ -263,7 +263,7 @@ class TimeFields(Fields): for i, TIND_i in enumerate(timeII): fieldI = pointerFields[:,:,i] if fieldI.shape[0] == fieldI.size: - fieldI = Utils.mkvc(fieldI,2) + fieldI = Utils.mkvc(fieldI,1) out[i] = func(fieldI, srcII, TIND_i) if out[i].ndim == 1: out[i] = out[i][:,np.newaxis,np.newaxis] From 14ee13fadb340ddbbc56e663aca1bc7c6b9e04f1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 10:18:26 -0700 Subject: [PATCH 22/49] remove redundant mkvc option --- SimPEG/Problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index e729ee10..87111b31 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -263,7 +263,7 @@ class TimeFields(Fields): for i, TIND_i in enumerate(timeII): fieldI = pointerFields[:,:,i] if fieldI.shape[0] == fieldI.size: - fieldI = Utils.mkvc(fieldI,1) + fieldI = Utils.mkvc(fieldI) out[i] = func(fieldI, srcII, TIND_i) if out[i].ndim == 1: out[i] = out[i][:,np.newaxis,np.newaxis] From 7e171ede05a06844b6323e09711ba6be9b4ec95d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 10:26:53 -0700 Subject: [PATCH 23/49] Move Fields Objects to their own file. --- SimPEG/Fields.py | 273 +++++++++++++++++++++++++++++++++++++++++++++ SimPEG/Problem.py | 275 +--------------------------------------------- 2 files changed, 274 insertions(+), 274 deletions(-) create mode 100644 SimPEG/Fields.py diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py new file mode 100644 index 00000000..f81cfb67 --- /dev/null +++ b/SimPEG/Fields.py @@ -0,0 +1,273 @@ +import Utils, numpy as np, scipy.sparse as sp + +class Fields(object): + """Fancy Field Storage + + u[:,'phi'] = phi + print u[src0,'phi'] + + """ + + knownFields = None #: Known fields, a dict with locations, e.g. {"e": "E", "phi": "CC"} + aliasFields = None #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]} + dtype = float #: dtype is the type of the storage matrix. This can be a dictionary. + + def __init__(self, mesh, survey, **kwargs): + self.survey = survey + self.mesh = mesh + Utils.setKwargs(self, **kwargs) + self._fields = {} + + if self.knownFields is None: + raise Exception('knownFields cannot be set to None') + if self.aliasFields is None: + self.aliasFields = {} + + allFields = [k for k in self.knownFields] + [a for a in self.aliasFields] + assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.' + self.startup() + + def startup(self): + pass + + @property + def approxSize(self): + """The approximate cost to storing all of the known fields.""" + sz = 0.0 + for f in self.knownFields: + loc =self.knownFields[f] + sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2) + return "%e MB"%sz + + def _storageShape(self, loc): + nSrc = self.survey.nSrc + + nP = {'CC': self.mesh.nC, + 'N': self.mesh.nN, + 'F': self.mesh.nF, + 'E': self.mesh.nE}[loc] + + return (nP, nSrc) + + def _initStore(self, name): + if name in self._fields: + return self._fields[name] + + assert name in self.knownFields, 'field name is not known.' + + loc = self.knownFields[name] + + if type(self.dtype) is dict: + dtype = self.dtype[name] + else: + dtype = self.dtype + field = np.zeros(self._storageShape(loc), dtype=dtype) + + self._fields[name] = field + + return field + + def _srcIndex(self, srcTestList): + if type(srcTestList) is slice: + ind = srcTestList + else: + if type(srcTestList) is not list: + srcTestList = [srcTestList] + for srcTest in srcTestList: + if srcTest not in self.survey.srcList: + raise KeyError('Invalid Source, not in survey list.') + + ind = np.in1d(self.survey.srcList, srcTestList) + return ind + + def _nameIndex(self, name, accessType): + + if type(name) is slice: + assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.' + name = None + + if name is None: + return + if accessType=='set' and name not in self.knownFields: + if name in self.aliasFields: + raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name) + else: + raise KeyError('Invalid field name (%s) for setter'%name) + + elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields): + raise KeyError('Invalid field name (%s) for getter'%name) + return name + + def _indexAndNameFromKey(self, key, accessType): + if type(key) is not tuple: + key = (key,) + if len(key) == 1: + key += (None,) + + assert len(key) == 2, 'must be [Src, fieldName]' + + srcTestList, name = key + name = self._nameIndex(name, accessType) + ind = self._srcIndex(srcTestList) + return ind, name + + def __setitem__(self, key, value): + ind, name = self._indexAndNameFromKey(key, 'set') + if name is None: + freq = key + assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.' + newFields = value + elif name in self.knownFields: + newFields = {name: value} + else: + raise Exception('Unknown setter') + + for name in newFields: + field = self._initStore(name) + self._setField(field, newFields[name], name, ind) + + def __getitem__(self, key): + ind, name = self._indexAndNameFromKey(key, 'get') + if name is None: + out = {} + for name in self._fields: + out[name] = self._getField(name, ind) + return out + return self._getField(name, ind) + + def _setField(self, field, val, name, ind): + if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1): + val = Utils.mkvc(val,2) + field[:,ind] = val + + def _getField(self, name, ind): + if name in self._fields: + out = self._fields[name][:,ind] + else: + # Aliased fields + alias, loc, func = self.aliasFields[name] + + srcII = np.array(self.survey.srcList)[ind] + if isinstance(srcII, np.ndarray): + srcII = srcII.tolist() + if len(srcII) == 1: + srcII = srcII[0] + + if type(func) is str: + assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' + func = getattr(self, func) + out = func(self._fields[alias][:,ind], srcII) + if out.shape[0] == out.size or out.ndim == 1: + out = Utils.mkvc(out,2) + return out + + def __contains__(self, other): + if other in self.aliasFields: + other = self.aliasFields[other][0] + return self._fields.__contains__(other) + + +class TimeFields(Fields): + """Fancy Field Storage for time domain problems + + u[:,'phi', timeInd] = phi + print u[src0,'phi'] + + """ + + def _storageShape(self, loc): + nP = {'CC': self.mesh.nC, + 'N': self.mesh.nN, + 'F': self.mesh.nF, + 'E': self.mesh.nE}[loc] + nSrc = self.survey.nSrc + nT = self.survey.prob.nT + 1 + return (nP, nSrc, nT) + + def _indexAndNameFromKey(self, key, accessType): + if type(key) is not tuple: + key = (key,) + if len(key) == 1: + key += (None,) + if len(key) == 2: + key += (slice(None,None,None),) + + assert len(key) == 3, 'must be [Src, fieldName, times]' + + srcTestList, name, timeInd = key + + name = self._nameIndex(name, accessType) + srcInd = self._srcIndex(srcTestList) + + return (srcInd, timeInd), name + + def _correctShape(self, name, ind, deflate=False): + srcInd, timeInd = ind + if name in self.knownFields: + loc = self.knownFields[name] + else: + loc = self.aliasFields[name][1] + nP, total_nSrc, total_nT = self._storageShape(loc) + nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum() + nT = np.ones(total_nT, dtype=bool)[timeInd].sum() + shape = nP, nSrc, nT + if deflate: + shape = tuple([s for s in shape if s > 1]) + if len(shape) == 1: + shape = shape + (1,) + return shape + + def _setField(self, field, val, name, ind): + srcInd, timeInd = ind + shape = self._correctShape(name, ind) + if Utils.isScalar(val): + field[:,srcInd,timeInd] = val + return + if val.size != np.array(shape).prod(): + raise ValueError('Incorrect size for data.') + correctShape = field[:,srcInd,timeInd].shape + field[:,srcInd,timeInd] = val.reshape(correctShape, order='F') + + def _getField(self, name, ind): + srcInd, timeInd = ind + + if name in self._fields: + out = self._fields[name][:,srcInd,timeInd] + else: + # Aliased fields + alias, loc, func = self.aliasFields[name] + if type(func) is str: + assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' + func = getattr(self, func) + pointerFields = self._fields[alias][:,srcInd,timeInd] + pointerShape = self._correctShape(alias, ind) + pointerFields = pointerFields.reshape(pointerShape, order='F') + + timeII = np.arange(self.survey.prob.nT + 1)[timeInd] + srcII = np.array(self.survey.srcList)[srcInd] + if isinstance(srcII, np.ndarray): + srcII = srcII.tolist() + if len(srcII) == 1: + srcII = srcII[0] + + if timeII.size == 1: + pointerShapeDeflated = self._correctShape(alias, ind, deflate=True) + pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F') + out = func(pointerFields, srcII, timeII) + else: #loop over the time steps + nT = pointerShape[2] + out = range(nT) + for i, TIND_i in enumerate(timeII): + fieldI = pointerFields[:,:,i] + if fieldI.shape[0] == fieldI.size: + fieldI = Utils.mkvc(fieldI) + out[i] = func(fieldI, srcII, TIND_i) + if out[i].ndim == 1: + out[i] = out[i][:,np.newaxis,np.newaxis] + elif out[i].ndim == 2: + out[i] = out[i][:,:,np.newaxis] + out = np.concatenate(out, axis=2) + + shape = self._correctShape(name, ind, deflate=True) + return out.reshape(shape, order='F') + diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 87111b31..75e9d4bf 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -1,280 +1,7 @@ import Utils, Survey, Models, numpy as np, scipy.sparse as sp Solver = Utils.SolverUtils.Solver import Maps, Mesh - - -class Fields(object): - """Fancy Field Storage - - u[:,'phi'] = phi - print u[src0,'phi'] - - """ - - knownFields = None #: Known fields, a dict with locations, e.g. {"e": "E", "phi": "CC"} - aliasFields = None #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]} - dtype = float #: dtype is the type of the storage matrix. This can be a dictionary. - - def __init__(self, mesh, survey, **kwargs): - self.survey = survey - self.mesh = mesh - Utils.setKwargs(self, **kwargs) - self._fields = {} - - if self.knownFields is None: - raise Exception('knownFields cannot be set to None') - if self.aliasFields is None: - self.aliasFields = {} - - allFields = [k for k in self.knownFields] + [a for a in self.aliasFields] - assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.' - self.startup() - - def startup(self): - pass - - @property - def approxSize(self): - """The approximate cost to storing all of the known fields.""" - sz = 0.0 - for f in self.knownFields: - loc =self.knownFields[f] - sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2) - return "%e MB"%sz - - def _storageShape(self, loc): - nSrc = self.survey.nSrc - - nP = {'CC': self.mesh.nC, - 'N': self.mesh.nN, - 'F': self.mesh.nF, - 'E': self.mesh.nE}[loc] - - return (nP, nSrc) - - def _initStore(self, name): - if name in self._fields: - return self._fields[name] - - assert name in self.knownFields, 'field name is not known.' - - loc = self.knownFields[name] - - if type(self.dtype) is dict: - dtype = self.dtype[name] - else: - dtype = self.dtype - field = np.zeros(self._storageShape(loc), dtype=dtype) - - self._fields[name] = field - - return field - - def _srcIndex(self, srcTestList): - if type(srcTestList) is slice: - ind = srcTestList - else: - if type(srcTestList) is not list: - srcTestList = [srcTestList] - for srcTest in srcTestList: - if srcTest not in self.survey.srcList: - raise KeyError('Invalid Source, not in survey list.') - - ind = np.in1d(self.survey.srcList, srcTestList) - return ind - - def _nameIndex(self, name, accessType): - - if type(name) is slice: - assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.' - name = None - - if name is None: - return - if accessType=='set' and name not in self.knownFields: - if name in self.aliasFields: - raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name) - else: - raise KeyError('Invalid field name (%s) for setter'%name) - - elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields): - raise KeyError('Invalid field name (%s) for getter'%name) - return name - - def _indexAndNameFromKey(self, key, accessType): - if type(key) is not tuple: - key = (key,) - if len(key) == 1: - key += (None,) - - assert len(key) == 2, 'must be [Src, fieldName]' - - srcTestList, name = key - name = self._nameIndex(name, accessType) - ind = self._srcIndex(srcTestList) - return ind, name - - def __setitem__(self, key, value): - ind, name = self._indexAndNameFromKey(key, 'set') - if name is None: - freq = key - assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.' - newFields = value - elif name in self.knownFields: - newFields = {name: value} - else: - raise Exception('Unknown setter') - - for name in newFields: - field = self._initStore(name) - self._setField(field, newFields[name], name, ind) - - def __getitem__(self, key): - ind, name = self._indexAndNameFromKey(key, 'get') - if name is None: - out = {} - for name in self._fields: - out[name] = self._getField(name, ind) - return out - return self._getField(name, ind) - - def _setField(self, field, val, name, ind): - if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1): - val = Utils.mkvc(val,2) - field[:,ind] = val - - def _getField(self, name, ind): - if name in self._fields: - out = self._fields[name][:,ind] - else: - # Aliased fields - alias, loc, func = self.aliasFields[name] - - srcII = np.array(self.survey.srcList)[ind] - if isinstance(srcII, np.ndarray): - srcII = srcII.tolist() - if len(srcII) == 1: - srcII = srcII[0] - - if type(func) is str: - assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' - func = getattr(self, func) - out = func(self._fields[alias][:,ind], srcII) - if out.shape[0] == out.size or out.ndim == 1: - out = Utils.mkvc(out,2) - return out - - def __contains__(self, other): - if other in self.aliasFields: - other = self.aliasFields[other][0] - return self._fields.__contains__(other) - - -class TimeFields(Fields): - """Fancy Field Storage for time domain problems - - u[:,'phi', timeInd] = phi - print u[src0,'phi'] - - """ - - def _storageShape(self, loc): - nP = {'CC': self.mesh.nC, - 'N': self.mesh.nN, - 'F': self.mesh.nF, - 'E': self.mesh.nE}[loc] - nSrc = self.survey.nSrc - nT = self.survey.prob.nT + 1 - return (nP, nSrc, nT) - - def _indexAndNameFromKey(self, key, accessType): - if type(key) is not tuple: - key = (key,) - if len(key) == 1: - key += (None,) - if len(key) == 2: - key += (slice(None,None,None),) - - assert len(key) == 3, 'must be [Src, fieldName, times]' - - srcTestList, name, timeInd = key - - name = self._nameIndex(name, accessType) - srcInd = self._srcIndex(srcTestList) - - return (srcInd, timeInd), name - - def _correctShape(self, name, ind, deflate=False): - srcInd, timeInd = ind - if name in self.knownFields: - loc = self.knownFields[name] - else: - loc = self.aliasFields[name][1] - nP, total_nSrc, total_nT = self._storageShape(loc) - nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum() - nT = np.ones(total_nT, dtype=bool)[timeInd].sum() - shape = nP, nSrc, nT - if deflate: - shape = tuple([s for s in shape if s > 1]) - if len(shape) == 1: - shape = shape + (1,) - return shape - - def _setField(self, field, val, name, ind): - srcInd, timeInd = ind - shape = self._correctShape(name, ind) - if Utils.isScalar(val): - field[:,srcInd,timeInd] = val - return - if val.size != np.array(shape).prod(): - raise ValueError('Incorrect size for data.') - correctShape = field[:,srcInd,timeInd].shape - field[:,srcInd,timeInd] = val.reshape(correctShape, order='F') - - def _getField(self, name, ind): - srcInd, timeInd = ind - - if name in self._fields: - out = self._fields[name][:,srcInd,timeInd] - else: - # Aliased fields - alias, loc, func = self.aliasFields[name] - if type(func) is str: - assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' - func = getattr(self, func) - pointerFields = self._fields[alias][:,srcInd,timeInd] - pointerShape = self._correctShape(alias, ind) - pointerFields = pointerFields.reshape(pointerShape, order='F') - - timeII = np.arange(self.survey.prob.nT + 1)[timeInd] - srcII = np.array(self.survey.srcList)[srcInd] - if isinstance(srcII, np.ndarray): - srcII = srcII.tolist() - if len(srcII) == 1: - srcII = srcII[0] - - if timeII.size == 1: - pointerShapeDeflated = self._correctShape(alias, ind, deflate=True) - pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F') - out = func(pointerFields, srcII, timeII) - else: #loop over the time steps - nT = pointerShape[2] - out = range(nT) - for i, TIND_i in enumerate(timeII): - fieldI = pointerFields[:,:,i] - if fieldI.shape[0] == fieldI.size: - fieldI = Utils.mkvc(fieldI) - out[i] = func(fieldI, srcII, TIND_i) - if out[i].ndim == 1: - out[i] = out[i][:,np.newaxis,np.newaxis] - elif out[i].ndim == 2: - out[i] = out[i][:,:,np.newaxis] - out = np.concatenate(out, axis=2) - - shape = self._correctShape(name, ind, deflate=True) - return out.reshape(shape, order='F') - - +from Fields import Fields, TimeFields class BaseProblem(object): """ From 116f7620a6beeee411d53a6ac68b75bd686fc528 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 10:32:48 -0700 Subject: [PATCH 24/49] Split up the tests --- .../Tests/{test_Survey.py => test_Fields.py} | 31 +++---------- SimPEG/Tests/test_SurveyAndData.py | 45 +++++++++++++++++++ 2 files changed, 51 insertions(+), 25 deletions(-) rename SimPEG/Tests/{test_Survey.py => test_Fields.py} (95%) create mode 100644 SimPEG/Tests/test_SurveyAndData.py diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Fields.py similarity index 95% rename from SimPEG/Tests/test_Survey.py rename to SimPEG/Tests/test_Fields.py index efc1e901..22189a15 100644 --- a/SimPEG/Tests/test_Survey.py +++ b/SimPEG/Tests/test_Fields.py @@ -1,7 +1,8 @@ import unittest from SimPEG import * -class DataAndFieldsTest(unittest.TestCase): + +class FieldsTest(unittest.TestCase): def setUp(self): mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) @@ -26,25 +27,6 @@ class DataAndFieldsTest(unittest.TestCase): self.mesh = mesh self.XYZ = XYZ - def test_overlappingFields(self): - self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey, - knownFields={'b':'F'}, - aliasFields={'b':['b',(lambda F, b, ind: b)]}) - - def test_data(self): - V = [] - for src in self.D.survey.srcList: - for rx in src.rxList: - v = np.random.rand(rx.nD) - V += [v] - self.D[src, rx] = v - self.assertTrue(np.all(v == self.D[src, rx])) - V = np.concatenate(V) - self.assertTrue(np.all(V == Utils.mkvc(self.D))) - - D2 = Survey.Data(self.D.survey, V) - self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D))) - def test_contains(self): F = self.F nSrc = F.survey.nSrc @@ -55,10 +37,10 @@ class DataAndFieldsTest(unittest.TestCase): self.assertTrue('b' not in F) self.assertTrue('e' in F) - def test_uniqueSrcs(self): - srcs = self.D.survey.srcList - srcs += [srcs[0]] - self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs) + def test_overlappingFields(self): + self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey, + knownFields={'b':'F'}, + aliasFields={'b':['b',(lambda F, b, ind: b)]}) def test_SetGet(self): F = self.F @@ -132,7 +114,6 @@ class FieldsTest_Alias(unittest.TestCase): Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc) srcList = [Src0,Src1,Src2,Src3,Src4] survey = Survey.BaseSurvey(srcList=srcList) - self.D = Survey.Data(survey) self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]}) self.Src0 = Src0 self.Src1 = Src1 diff --git a/SimPEG/Tests/test_SurveyAndData.py b/SimPEG/Tests/test_SurveyAndData.py new file mode 100644 index 00000000..bb6b4645 --- /dev/null +++ b/SimPEG/Tests/test_SurveyAndData.py @@ -0,0 +1,45 @@ +import unittest +from SimPEG import * + +class TestData(unittest.TestCase): + + def setUp(self): + mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0.]) + srcLoc = np.r_[0,0,0.] + rxList0 = Survey.BaseRx(XYZ, 'exi') + Src0 = Survey.BaseSrc([rxList0], loc=srcLoc) + rxList1 = Survey.BaseRx(XYZ, 'bxi') + Src1 = Survey.BaseSrc([rxList1], loc=srcLoc) + rxList2 = Survey.BaseRx(XYZ, 'bxi') + Src2 = Survey.BaseSrc([rxList2], loc=srcLoc) + rxList3 = Survey.BaseRx(XYZ, 'bxi') + Src3 = Survey.BaseSrc([rxList3], loc=srcLoc) + Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc) + srcList = [Src0,Src1,Src2,Src3,Src4] + survey = Survey.BaseSurvey(srcList=srcList) + self.D = Survey.Data(survey) + + def test_data(self): + V = [] + for src in self.D.survey.srcList: + for rx in src.rxList: + v = np.random.rand(rx.nD) + V += [v] + self.D[src, rx] = v + self.assertTrue(np.all(v == self.D[src, rx])) + V = np.concatenate(V) + self.assertTrue(np.all(V == Utils.mkvc(self.D))) + + D2 = Survey.Data(self.D.survey, V) + self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D))) + + def test_uniqueSrcs(self): + srcs = self.D.survey.srcList + srcs += [srcs[0]] + self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs) + + +if __name__ == '__main__': + unittest.main() From 59fcd3925f1df6ca74ec9672cde0b57b2d247390 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 11:03:47 -0700 Subject: [PATCH 25/49] getSourceIndex --- SimPEG/Survey.py | 14 ++++++++++++-- SimPEG/Tests/test_SurveyAndData.py | 10 ++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 8abb547a..b6a17f28 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -1,4 +1,4 @@ -import Utils, numpy as np, scipy.sparse as sp +import Utils, numpy as np, scipy.sparse as sp, uuid class BaseRx(object): @@ -13,6 +13,7 @@ class BaseRx(object): storeProjections = True #: Store calls to getP (organized by mesh) def __init__(self, locs, rxType, **kwargs): + self.uid = str(uuid.uuid4()) self.locs = locs self.rxType = rxType self._Ps = {} @@ -124,7 +125,7 @@ class BaseSrc(object): for rx in rxList: assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__ assert len(set(rxList)) == len(rxList), 'The rxList must be unique' - + self.uid = str(uuid.uuid4()) self.rxList = rxList Utils.setKwargs(self, **kwargs) @@ -144,6 +145,7 @@ class Data(object): """Fancy data storage by Src and Rx""" def __init__(self, survey, v=None): + self.uid = str(uuid.uuid4()) self.survey = survey self._dataDict = {} for src in self.survey.srcList: @@ -225,6 +227,14 @@ class BaseSurvey(object): assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__ assert len(set(value)) == len(value), 'The srcList must be unique' self._srcList = value + self._sourceOrder = dict() + [self._sourceOrder.setdefault(src.uid, ii) for ii, src in enumerate(self._srcList)] + + def getSourceIndex(self, sources): + inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources) + if None in inds: + raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds)) + return inds @property def prob(self): diff --git a/SimPEG/Tests/test_SurveyAndData.py b/SimPEG/Tests/test_SurveyAndData.py index bb6b4645..6feccc71 100644 --- a/SimPEG/Tests/test_SurveyAndData.py +++ b/SimPEG/Tests/test_SurveyAndData.py @@ -40,6 +40,16 @@ class TestData(unittest.TestCase): srcs += [srcs[0]] self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs) + def test_sourceIndex(self): + survey = self.D.survey + srcs = survey.srcList + assert survey.getSourceIndex([srcs[1],srcs[0]]) == [1,0] + assert survey.getSourceIndex([srcs[1],srcs[2],srcs[2]]) == [1,2,2] + SrcNotThere = Survey.BaseSrc(srcs[0].rxList, loc=np.r_[0,0,0]) + self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere]) + self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere]) + + if __name__ == '__main__': unittest.main() From de27c4e4ec94bcd18b5cb27bed81bf3d86aa3347 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 11:17:56 -0700 Subject: [PATCH 26/49] fixes #99 --- SimPEG/Fields.py | 8 +------- SimPEG/Survey.py | 5 +++++ SimPEG/Tests/test_SurveyAndData.py | 2 -- 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py index f81cfb67..edd4cd92 100644 --- a/SimPEG/Fields.py +++ b/SimPEG/Fields.py @@ -71,13 +71,7 @@ class Fields(object): if type(srcTestList) is slice: ind = srcTestList else: - if type(srcTestList) is not list: - srcTestList = [srcTestList] - for srcTest in srcTestList: - if srcTest not in self.survey.srcList: - raise KeyError('Invalid Source, not in survey list.') - - ind = np.in1d(self.survey.srcList, srcTestList) + ind = self.survey.getSourceIndex(srcTestList) return ind def _nameIndex(self, name, accessType): diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index b6a17f28..0fdb0cd1 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -231,6 +231,11 @@ class BaseSurvey(object): [self._sourceOrder.setdefault(src.uid, ii) for ii, src in enumerate(self._srcList)] def getSourceIndex(self, sources): + if type(sources) is not list: + sources = [sources] + for src in sources: + if getattr(src,'uid',None) is None: + raise KeyError('Source does not have a uid: %s'%str(src)) inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources) if None in inds: raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds)) diff --git a/SimPEG/Tests/test_SurveyAndData.py b/SimPEG/Tests/test_SurveyAndData.py index 6feccc71..f02f14e8 100644 --- a/SimPEG/Tests/test_SurveyAndData.py +++ b/SimPEG/Tests/test_SurveyAndData.py @@ -49,7 +49,5 @@ class TestData(unittest.TestCase): self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere]) self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere]) - - if __name__ == '__main__': unittest.main() From 2827e85330312d11b01cb4daf0350a997e9a82a5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 May 2015 11:36:56 -0700 Subject: [PATCH 27/49] fix mkvc to return an (n,1) array for consistency --- SimPEG/Fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py index edd4cd92..1801bedc 100644 --- a/SimPEG/Fields.py +++ b/SimPEG/Fields.py @@ -254,7 +254,7 @@ class TimeFields(Fields): for i, TIND_i in enumerate(timeII): fieldI = pointerFields[:,:,i] if fieldI.shape[0] == fieldI.size: - fieldI = Utils.mkvc(fieldI) + fieldI = Utils.mkvc(fieldI, 2) out[i] = func(fieldI, srcII, TIND_i) if out[i].ndim == 1: out[i] = out[i][:,np.newaxis,np.newaxis] From a953a52ccc528726ea03ae33c514ce0fec5d1c69 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 31 May 2015 10:43:23 -0700 Subject: [PATCH 28/49] 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 29/49] 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 30/49] 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 31/49] 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 32/49] 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 33/49] 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 34/49] 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 35/49] 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 36/49] 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) From 658d481dd6d49a4e013d7e5a3da7f8d922ccdcc0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 15:41:25 -0700 Subject: [PATCH 37/49] PropMap Bug fix with unmapped derivs. --- SimPEG/PropMaps.py | 5 ++++- SimPEG/Tests/test_PropMaps.py | 8 +++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 3bc11d3c..65b735b8 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -68,7 +68,10 @@ class Property(object): if mapping is None and prop.propertyLink is not None: linkName, linkMapClass = prop.propertyLink - linkMap = linkMapClass(None) * getattr(self, '%sMap'%linkName) + linkedMap = getattr(self, '%sMap'%linkName) + if linkedMap is None: + return None + linkMap = linkMapClass(None) * linkedMap m = getattr(self, '%s'%linkName) return linkMap.deriv( m ) diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index bb5788bc..75c747cf 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -10,7 +10,8 @@ class MyPropMap(Maps.PropMap): 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) + mu = Maps.Property("Mu", defaultVal=mu_0, propertyLink=('mui', Maps.ReciprocalMap)) + mui = Maps.Property("Mu", defaultVal=1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) class TestPropMaps(unittest.TestCase): @@ -102,6 +103,11 @@ class TestPropMaps(unittest.TestCase): assert pm.sigmaDeriv is not None assert pm.rhoDeriv is not None + assert pm.muMap is None + assert pm.muDeriv is None + assert pm.muiMap is None + assert pm.muiDeriv is None + PM = MyReciprocalPropMap([('rho', iMap)]) pm = PM(np.r_[1,2.,3]) # print pm.sigma From 4fa4ef643d7fbe64118867c0a1300b7a45b38004 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 15:58:26 -0700 Subject: [PATCH 38/49] updates to recursion issues in reciprocal maps --- SimPEG/PropMaps.py | 2 ++ SimPEG/Tests/test_PropMaps.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 65b735b8..6ce61ccb 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -52,6 +52,8 @@ class Property(object): if mapping is None and prop.propertyLink is not None: linkName, linkMapClass = prop.propertyLink linkMap = linkMapClass(None) + if getattr(self, '%sMap'%linkName, None) is None: + return prop.defaultVal m = getattr(self, '%s'%linkName) return linkMap * m diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index 75c747cf..c2726d71 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -103,6 +103,8 @@ class TestPropMaps(unittest.TestCase): assert pm.sigmaDeriv is not None assert pm.rhoDeriv is not None + assert pm.mu == mu_0 + assert pm.mui == 1.0/mu_0 assert pm.muMap is None assert pm.muDeriv is None assert pm.muiMap is None From f90637509e4f203677ca28491de8b498179749ee Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 16:07:05 -0700 Subject: [PATCH 39/49] Fields takes srcList --- SimPEG/Fields.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py index 1801bedc..d9da317d 100644 --- a/SimPEG/Fields.py +++ b/SimPEG/Fields.py @@ -141,11 +141,8 @@ class Fields(object): # Aliased fields alias, loc, func = self.aliasFields[name] - srcII = np.array(self.survey.srcList)[ind] - if isinstance(srcII, np.ndarray): - srcII = srcII.tolist() - if len(srcII) == 1: - srcII = srcII[0] + srcII = np.array(self.survey.srcList)[ind] + srcII = srcII.tolist() if type(func) is str: assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.' From 4df383ccec5c745c81036dc7e407c648fd035c4c Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 1 Jun 2015 16:17:05 -0700 Subject: [PATCH 40/49] fixed test_Fields --- SimPEG/Tests/test_Fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Tests/test_Fields.py b/SimPEG/Tests/test_Fields.py index 22189a15..bd75d712 100644 --- a/SimPEG/Tests/test_Fields.py +++ b/SimPEG/Tests/test_Fields.py @@ -147,7 +147,7 @@ class FieldsTest_Alias(unittest.TestCase): def test_aliasFunction(self): def alias(e, ind): - self.assertTrue(ind is self.Src0) + self.assertTrue(ind[0] is self.Src0) return self.F.mesh.edgeCurl * e F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]}) e = np.random.rand(F.mesh.nE,1) From 5caf2371212f27adea43b65a44db400bf6d8c27d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 16:27:31 -0700 Subject: [PATCH 41/49] updates to the time fields side of things as well --- SimPEG/Fields.py | 7 ++----- SimPEG/Tests/test_Fields.py | 2 +- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py index d9da317d..9baa32f8 100644 --- a/SimPEG/Fields.py +++ b/SimPEG/Fields.py @@ -235,11 +235,8 @@ class TimeFields(Fields): pointerFields = pointerFields.reshape(pointerShape, order='F') timeII = np.arange(self.survey.prob.nT + 1)[timeInd] - srcII = np.array(self.survey.srcList)[srcInd] - if isinstance(srcII, np.ndarray): - srcII = srcII.tolist() - if len(srcII) == 1: - srcII = srcII[0] + srcII = np.array(self.survey.srcList)[srcInd] + srcII = srcII.tolist() if timeII.size == 1: pointerShapeDeflated = self._correctShape(alias, ind, deflate=True) diff --git a/SimPEG/Tests/test_Fields.py b/SimPEG/Tests/test_Fields.py index bd75d712..629189b8 100644 --- a/SimPEG/Tests/test_Fields.py +++ b/SimPEG/Tests/test_Fields.py @@ -343,7 +343,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase): count = [0] def alias(e, srcInd, timeInd): count[0] += 1 - self.assertTrue(srcInd is self.Src0) + self.assertTrue(srcInd[0] is self.Src0) return self.F.mesh.edgeCurl * e F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]}) e = np.random.rand(F.mesh.nE,1,nT) From 165afb958f3164c957f82e477950045b020f8c3c Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 22:10:00 -0700 Subject: [PATCH 42/49] Make reciprocal test more reliable. --- SimPEG/Tests/test_maps.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py index 437e51fa..a80bc40b 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", "ReciprocalMap"] -MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "ReciprocalMap"] +MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"] +MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"] class MapTests(unittest.TestCase): @@ -30,8 +30,8 @@ class MapTests(unittest.TestCase): self.assertTrue(maps.test()) - def test_transforms_logMap(self): - # Note that log maps can be kinda finicky, so we are being explicit about the random seed. + def test_transforms_logMap_reciprocalMap(self): + # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed. v2 = np.r_[ 0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519, 0.79754415] dv2 = np.r_[ 0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762, 0.44252488] v3 = np.r_[ 0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109, 0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864, 0.84130763, 0.22123854] @@ -41,6 +41,11 @@ class MapTests(unittest.TestCase): maps = Maps.LogMap(self.mesh3) self.assertTrue(maps.test(v3, dx=dv3)) + maps = Maps.ReciprocalMap(self.mesh2) + self.assertTrue(maps.test(v2, dx=dv2)) + maps = Maps.ReciprocalMap(self.mesh3) + self.assertTrue(maps.test(v3, dx=dv3)) + def test_Mesh2MeshMap(self): maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2]) self.assertTrue(maps.test()) From 6b4dede7a4467d9127d54a016e6f3b80db43e368 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 1 Jun 2015 22:11:14 -0700 Subject: [PATCH 43/49] remove redundant code --- SimPEG/PropMaps.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 6ce61ccb..d7bd2ba2 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -151,10 +151,6 @@ class _PropMapMetaClass(type): 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) From 1b7ad56e94f9478d865e0a119d9e4ed5eadafa03 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 3 Jun 2015 15:46:38 -0700 Subject: [PATCH 44/49] updates to #104 --- SimPEG/PropMaps.py | 3 +++ SimPEG/Tests/test_PropMaps.py | 22 ++++++++++++++++++++++ 2 files changed, 25 insertions(+) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index d7bd2ba2..bddd3df7 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -238,3 +238,6 @@ class PropMap(object): def __call__(self, vec): return self.PropModel(self, vec) + def __contains__(self, val): + activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None] + return val in activeMaps diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index c2726d71..6861f09b 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -36,6 +36,10 @@ class TestPropMaps(unittest.TestCase): assert PM.muMap is None assert PM.muIndex is None + assert 'sigma' in PM + assert 'mu' not in PM + assert 'mui' not in PM + m = PM(np.r_[1.,2,3]) assert m.mu == mu_0 assert m.muModel is None @@ -67,6 +71,10 @@ class TestPropMaps(unittest.TestCase): assert pm.nP == 6 + assert 'sigma' in PM + assert 'mu' in PM + assert 'mui' not in PM + 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]) @@ -83,6 +91,10 @@ class TestPropMaps(unittest.TestCase): assert pm.nP == 3 + assert 'sigma' in PM + assert 'mu' in PM + assert 'mui' not in PM + 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]) @@ -103,6 +115,11 @@ class TestPropMaps(unittest.TestCase): assert pm.sigmaDeriv is not None assert pm.rhoDeriv is not None + assert 'sigma' in PM + assert 'rho' not in PM + assert 'mu' not in PM + assert 'mui' not in PM + assert pm.mu == mu_0 assert pm.mui == 1.0/mu_0 assert pm.muMap is None @@ -121,6 +138,11 @@ class TestPropMaps(unittest.TestCase): assert pm.sigmaDeriv is not None assert pm.rhoDeriv is not None + assert 'sigma' not in PM + assert 'rho' in PM + assert 'mu' not in PM + assert 'mui' not in PM + self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)]) self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)]) From 838ee6e09b47f7b716a359bbf55cb2cff544b4bb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 3 Jun 2015 15:59:20 -0700 Subject: [PATCH 45/49] updates to adding projections to the models --- SimPEG/PropMaps.py | 18 ++++++++++++++++-- SimPEG/Tests/test_PropMaps.py | 14 ++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index bddd3df7..c0913cae 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -1,4 +1,4 @@ -import Utils, Maps, numpy as np +import Utils, Maps, numpy as np, scipy.sparse as sp class Property(object): @@ -87,10 +87,23 @@ class Property(object): mapping = getattr(self, '%sMap'%prop.name) if mapping is None: return None - index = getattr(self.propMap, '_%sIndex'%prop.name, slice(None)) + index = getattr(self.propMap, '%sIndex'%prop.name) return self.vector[index] return property(fget=fget) + def _getModelProjProperty(self): + prop = self + def fget(self): + mapping = getattr(self, '%sMap'%prop.name) + if mapping is None: + return None + inds = getattr(self.propMap, '%sIndex'%prop.name) + if type(inds) is slice: + inds = range(*inds.indices(self.nP)) + nI, nP = len(inds),self.nP + return sp.csr_matrix((np.ones(nI), (range(nI), inds) ), shape=(nI, nP)) + return property(fget=fget) + def _getModelMapProperty(self): prop = self def fget(self): @@ -167,6 +180,7 @@ class _PropMapMetaClass(type): attrs[attr ] = prop._getProperty() attrs[attr + 'Map' ] = prop._getModelMapProperty() + attrs[attr + 'Proj' ] = prop._getModelProjProperty() attrs[attr + 'Model'] = prop._getModelProperty() attrs[attr + 'Deriv'] = prop._getModelDerivProperty() diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index 6861f09b..93ef4129 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -100,6 +100,20 @@ class TestPropMaps(unittest.TestCase): assert np.all(pm.muModel == [1,2,3]) assert np.all(pm.mu == [1,2,3]) + def test_Projections(self): + m = Mesh.TensorMesh((3,)) + iMap = Maps.IdentityMap(m) + PM = MyReciprocalPropMap([('sigma', iMap)]) + v = np.r_[1,2.,3] + pm = PM(v) + + assert pm.sigmaProj is not None + assert pm.rhoProj is None + assert pm.muProj is None + assert pm.muiProj is None + + assert np.all(pm.sigmaProj * v == pm.sigmaModel) + def test_Links(self): m = Mesh.TensorMesh((3,)) expMap = Maps.ExpMap(m) From 9ec2fa5e794a46b05849b69d9d869211e02b69fd Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 3 Jun 2015 16:05:14 -0700 Subject: [PATCH 46/49] ask 'mu' in self.curModel ? --- SimPEG/PropMaps.py | 3 +++ SimPEG/Tests/test_PropMaps.py | 25 +++++++++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index c0913cae..995216f7 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -132,6 +132,9 @@ class PropModel(object): self._nP = len(set(inds)) return self._nP + def __contains__(self, val): + return val in self.propMap + _PROPMAPCLASSREGISTRY = {} diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py index 93ef4129..b4c0eb8d 100644 --- a/SimPEG/Tests/test_PropMaps.py +++ b/SimPEG/Tests/test_PropMaps.py @@ -41,6 +41,11 @@ class TestPropMaps(unittest.TestCase): assert 'mui' not in PM m = PM(np.r_[1.,2,3]) + + assert 'sigma' in m + assert 'mu' not in m + assert 'mui' not in m + assert m.mu == mu_0 assert m.muModel is None assert m.muMap is None @@ -75,6 +80,10 @@ class TestPropMaps(unittest.TestCase): assert 'mu' in PM assert 'mui' not in PM + assert 'sigma' in pm + assert 'mu' in pm + assert 'mui' not in pm + 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]) @@ -95,6 +104,10 @@ class TestPropMaps(unittest.TestCase): assert 'mu' in PM assert 'mui' not in PM + assert 'sigma' in pm + assert 'mu' in pm + assert 'mui' not in pm + 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]) @@ -134,6 +147,12 @@ class TestPropMaps(unittest.TestCase): assert 'mu' not in PM assert 'mui' not in PM + + assert 'sigma' in pm + assert 'rho' not in pm + assert 'mu' not in pm + assert 'mui' not in pm + assert pm.mu == mu_0 assert pm.mui == 1.0/mu_0 assert pm.muMap is None @@ -157,6 +176,12 @@ class TestPropMaps(unittest.TestCase): assert 'mu' not in PM assert 'mui' not in PM + + assert 'sigma' not in pm + assert 'rho' in pm + assert 'mu' not in pm + assert 'mui' not in pm + self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)]) self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)]) From 686598cc8f5dbb77a93b87658b45918206d7bcdc Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 3 Jun 2015 23:37:36 -0700 Subject: [PATCH 47/49] Add ignore values --- SimPEG/Mesh/View.py | 1 + 1 file changed, 1 insertion(+) diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 015b21aa..38432d3c 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -328,6 +328,7 @@ class TensorView(object): v = getattr(np,view)(v) # e.g. np.real(v) if clim is None: clim = [v.min(),v.max()] + v = np.ma.masked_where(np.isnan(v), v) out += (ax.pcolormesh(self.vectorNx, self.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) elif view in ['vec']: U, V = self.r(v.reshape((self.nC,-1), order='F'), 'CC', 'CC', 'M') From 08c9013fd11585351f7547b54fb228c0ac99500e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 5 Jun 2015 10:29:11 -0700 Subject: [PATCH 48/49] updates to inversion directives --- SimPEG/Directives.py | 51 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 10 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e84a7d39..eb9830b5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -146,21 +146,52 @@ class BetaSchedule(InversionDirective): -class SaveModelEveryIteration(InversionDirective): - """SaveModelEveryIteration""" +class _SaveEveryIteration(InversionDirective): + @property + def name(self): + if getattr(self, '_name', None) is None: + self._name = 'InversionModel' + return self._name + @name.setter + def name(self, value): + self._name = value @property - def modelName(self): - if getattr(self, '_modelName', None) is None: + def fileName(self): + if getattr(self, '_fileName', None) is None: from datetime import datetime - self._modelName = 'inversionModel-%s'%datetime.now().strftime('%Y-%m-%d') - return self._modelName - @modelName.setter - def modelName(self, value): - self._modelName = value + self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M')) + return self._fileName + @fileName.setter + def fileName(self, value): + self._fileName = value + + +class SaveModelEveryIteration(_SaveEveryIteration): + """SaveModelEveryIteration""" + + def initialize(self): + print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName def endIter(self): - np.save('%03d-%s' % (self.opt.iter, self.modelName), self.opt.xc) + np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc) + + +class SaveOutputEveryIteration(_SaveEveryIteration): + """SaveModelEveryIteration""" + + def initialize(self): + print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName + f = open(self.fileName+'.txt', 'w') + f.write(" # beta phi_d phi_m f\n") + f.close() + + def endIter(self): + f = open(self.fileName+'.txt', 'a') + f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f)) + f.close() + + # class UpdateReferenceModel(Parameter): From 6d3d8d78b601a11c52e90f87367f7b40e4d537a5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 5 Jun 2015 15:50:00 -0700 Subject: [PATCH 49/49] updates to target misfit --- SimPEG/Directives.py | 15 +++++++++++++++ SimPEG/Optimization.py | 4 ++++ 2 files changed, 19 insertions(+) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index eb9830b5..3ec9e96e 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,6 +144,21 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor +class TargetMisfit(InversionDirective): + + @property + def target(self): + if getattr(self, '_target', None) is None: + self._target = self.survey.nD + return self._target + @target.setter + def target(self, val): + self._target = val + + def endIter(self): + if self.invProb.phi_d < self.target: + self.opt.stopNextIteration = True + class _SaveEveryIteration(InversionDirective): diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 549cb8ae..4f2cb062 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -97,6 +97,8 @@ class Minimize(object): tolG = 1e-1 #: Tolerance on gradient norm eps = 1e-5 #: Small value + stopNextIteration = False #: Stops the optimization program nicely. + debug = False #: Print debugging information debugLS = False #: Print debugging information for the line-search @@ -186,6 +188,7 @@ class Minimize(object): xt, caught = self.modifySearchDirectionBreak(p) if not caught: return self.xc self.doEndIteration(xt) + if self.stopNextIteration: break self.printDone() self.finish() @@ -210,6 +213,7 @@ class Minimize(object): self.iter = 0 self.iterLS = 0 + self.stopNextIteration = False x0 = self.projection(x0) # ensure that we start of feasible. self.x0 = x0