From f20fcb4504ef3f229d951982d6bc39a09c58da6a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Tue, 24 May 2016 08:45:12 -0700 Subject: [PATCH] Minor type error based upon numpy version Cross gradient? --- SimPEG/EM/Static/DC/ProblemDC.py | 2 +- SimPEG/EM/Static/IP/ProblemIP.py | 2 +- SimPEG/EM/Static/SIP/Regularization.py | 116 ++++++++++++++++++++----- 3 files changed, 98 insertions(+), 22 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 0dd4f259..51b528a5 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -77,7 +77,7 @@ class BaseDCProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += df_dmT + du_dmT + Jtv += (df_dmT + du_dmT).astype(float) return Utils.mkvc(Jtv) diff --git a/SimPEG/EM/Static/IP/ProblemIP.py b/SimPEG/EM/Static/IP/ProblemIP.py index 4ff81990..a637516a 100644 --- a/SimPEG/EM/Static/IP/ProblemIP.py +++ b/SimPEG/EM/Static/IP/ProblemIP.py @@ -89,7 +89,7 @@ class BaseIPProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += df_dmT + du_dmT + Jtv += (df_dmT + du_dmT).astype(float) # Conductivity ((d u / d log sigma).T) if self._formulation is 'EB': return -Utils.mkvc(Jtv) diff --git a/SimPEG/EM/Static/SIP/Regularization.py b/SimPEG/EM/Static/SIP/Regularization.py index a327ce4d..3c20e6fe 100644 --- a/SimPEG/EM/Static/SIP/Regularization.py +++ b/SimPEG/EM/Static/SIP/Regularization.py @@ -89,6 +89,11 @@ class MultiRegularization(Simple): self._W = sp.vstack(wlist) return self._W + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + @Utils.timeIt def _evalSmall(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) @@ -102,27 +107,98 @@ class MultiRegularization(Simple): r = self.Wsmooth * ( self.mapping * m) return 0.5 * r.dot(r) - # TODO: Implement Cross Gradients.. - # @Utils.timeIt - # def _evalCross(self, m): - # if self.crossgrad == False: - # return 0. - # elif self.crossgrad == True: - # M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + def cross(a,b): + ax, ay, az = a[0], a[1], a[2] + bx, by, bz = b[0], b[1], b[2] + cx = ay*bz - az*by + cy = az*bx - ax*bz + cz = ax*by - ay*bx + return [cx, cy, cz] - # for imodel in range(self.nModels): - # ux.append(self.regmesh.aveFx2CC*self.regmesh.wx[imodel]*M[:,imodel]) - # uy.append(self.regmesh.aveFy2CC*self.regmesh.wy[imodel]*M[:,imodel]) - # uz.append(self.regmesh.aveFz2CC*self.regmesh.wz[imodel]*M[:,imodel]) + # TODO: Implement Cross Gradients.. + @Utils.timeIt + def _evalCross(self, m): + if self.crossgrad == False: + return 0. + elif self.crossgrad == True: + M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + + ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0] + ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0] + az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0] + bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ab + out_ab = cross([ax, ay, az], [bx, by, bz]) + r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross) + + if self.nModels == 3: + cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ac + out_ac = cross([ax, ay, az], [cx, cy, cz]) + #bc + out_bc = cross([bx, by, bz], [cx, cy, cz]) + r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)] + + return 0.5 * r.dot(r) + + @Utils.timeIt + def evalDeriv(self, m): + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + if self.crossgrad==True: + deriv += self._evalCrossDeriv(m) + return deriv + + @Utils.timeIt + def _evalCrossDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @Utils.timeIt + def eval2Deriv(self, m, v=None): + """ + Second derivative + + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :rtype: scipy.sparse.csr_matrix or numpy.ndarray + :return: WtW or WtW*v + + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the second derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W} + + """ + mD = self.mapping.deriv(m - self.mref) + if v is None: + return mD.T * self.W.T * self.W * mD + + return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - # ax, ay, az = ux[0], uy[0], uz[0] - # for imodel in range(1,self.nModels): - # bx, by, bz = ux[imodel], uy[imodel], uz[imodel] - # cx = ay*bz - az*by - # cy = az*bx - ax*bz - # cz = ax*by - ay*bx - # ax, ay, az = cx.copy(), cy.copy(), cz.copy() - # r = np.r_[ax, ay, az]*np.sqrt(self.betacross) - # return 0.5 * r.dot(r)