From fc4294eb4dc6fe3f4f43c1c413ea5ccfc5011bcc Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@visiblegeology.com>
Date: Thu, 24 Oct 2013 15:11:48 -0700
Subject: [PATCH] Added A sample Linear problem that runs!

---
 SimPEG/forward/LinearProblem.py         | 91 +++++++++++++++++++++++++
 SimPEG/forward/Problem.py               | 22 +++---
 SimPEG/inverse/Inversion.py             | 11 ++-
 SimPEG/inverse/Optimize.py              |  2 +-
 SimPEG/regularization/Regularization.py | 20 +++---
 5 files changed, 124 insertions(+), 22 deletions(-)
 create mode 100644 SimPEG/forward/LinearProblem.py

diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py
new file mode 100644
index 00000000..d8f3621e
--- /dev/null
+++ b/SimPEG/forward/LinearProblem.py
@@ -0,0 +1,91 @@
+import numpy as np
+from SimPEG.mesh import TensorMesh
+from SimPEG.forward import Problem
+from SimPEG.regularization import Regularization
+from SimPEG.inverse import *
+import matplotlib.pyplot as plt
+
+N = 100
+h = np.ones(N)/N
+M = TensorMesh([h])
+
+nk = 20
+jk = np.linspace(1.,20.,nk)
+p = -0.25
+q = 0.25
+
+
+
+g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx)
+
+G = np.empty((nk, M.nC))
+
+for i in range(nk):
+    G[i,:] = g(i)
+
+
+
+plt.figure(1)
+for i in range(nk):
+    plt.plot(G[i,:])
+
+
+m_true = np.zeros(M.nC)
+m_true[M.vectorCCx > 0.3] = 1.
+m_true[M.vectorCCx > 0.45] = -0.5
+m_true[M.vectorCCx > 0.6] = 0
+
+
+d_true = G.dot(m_true)
+noise = 0.1 * np.random.rand(d_true.size)
+
+d_obs = d_true + noise
+
+# plt.figure(3)
+# plt.plot(d_true,'-o')
+# plt.plot(d_obs,'r-o')
+
+
+
+class LinearProblem(Problem):
+    """docstring for LinearProblem"""
+
+    def dpred(self, m, u=None):
+        return self.G.dot(m)
+
+    def J(self, m, v, u=None):
+        return G.dot(v)
+
+    def Jt(self, m, v, u=None):
+        return G.T.dot(v)
+
+    def modelTransform(self, m):
+        return m
+
+    def modelTransformDeriv(self, m):
+        return sp.eye(m.size)
+
+prob = LinearProblem(M)
+prob.G = G
+prob.dobs = d_obs
+prob.std = np.ones_like(d_obs)*0.1
+
+reg = Regularization(M)
+
+opt = InexactGaussNewton(maxIter=10)
+
+inv = Inversion(prob,reg,opt)
+
+m0 = np.zeros_like(m_true)
+
+mrec = inv.run(m0)
+
+
+plt.figure(2)
+
+plt.plot(M.vectorCCx, m_true, 'b-')
+plt.plot(M.vectorCCx, mrec, 'r-')
+
+
+
+plt.show()
diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py
index 0333b8b4..32dad675 100644
--- a/SimPEG/forward/Problem.py
+++ b/SimPEG/forward/Problem.py
@@ -82,6 +82,17 @@ class Problem(object):
     def dobs(self, value):
         self._dobs = value
 
+    def dpred(self, m, u=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pu(m)
+        """
+        if u is None:
+            u = self.field(m)
+        return self.P*u
+
     def misfit(self, m, u=None):
         """
             :param numpy.array m: geophysical model
@@ -152,17 +163,6 @@ class Problem(object):
         """
         pass
 
-    def dpred(self, m, u=None):
-        """
-            Predicted data.
-
-            .. math::
-                d_\\text{pred} = Pu(m)
-        """
-        if u is None:
-            u = self.field(m)
-        return self.P*u
-
     def modelTransform(self, m):
         """
             :param numpy.array m: model
diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py
index fd166318..cf96676f 100644
--- a/SimPEG/inverse/Inversion.py
+++ b/SimPEG/inverse/Inversion.py
@@ -25,8 +25,15 @@ class Inversion(object):
 
     @property
     def phi_d_target(self):
+        """
+        target for phi_d
+
+        By default this is the number of data.
+
+        Note that we do not set the target if it is None, but we return the default value.
+        """
         if getattr(self, '_phi_d_target', None) is None:
-            return self.prob.dobs.size
+            return self.prob.dobs.size #
         return self._phi_d_target
     @phi_d_target.setter
     def phi_d_target(self, value):
@@ -43,7 +50,7 @@ class Inversion(object):
         return m
 
     def getBeta(self):
-        return 1e2
+        return 1e-2
 
     def stoppingCriteria(self):
         self._STOP = np.zeros(2,dtype=bool)
diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py
index 2545a381..4ac96d6a 100644
--- a/SimPEG/inverse/Optimize.py
+++ b/SimPEG/inverse/Optimize.py
@@ -22,7 +22,7 @@ class Minimize(object):
     tolF = 1e-4
     tolX = 1e-4
     tolG = 1e-4
-    eps = 1e-16
+    eps = 1e-5
 
     printIter = [] # push to here if you want to print these on iter
 
diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py
index 7e211df6..6f5970e6 100644
--- a/SimPEG/regularization/Regularization.py
+++ b/SimPEG/regularization/Regularization.py
@@ -13,29 +13,33 @@ class Regularization(object):
     def mref(self, value):
         self._mref = value
 
+    @property
+    def Ws(self):
+        if getattr(self,'_Ws', None) is None:
+            self._Ws = sdiag(self.mesh.vol)
+        return self._Ws
+
     @property
     def Wx(self):
         if getattr(self, '_Wx', None) is None:
-            self._Wx = self.mesh.cellGradx
+            a = self.mesh.r(self.mesh.area,'F','Fx','V')
+            self._Wx = sdiag(a)*self.mesh.cellGradx
         return self._Wx
 
     @property
     def Wy(self):
         if getattr(self, '_Wy', None) is None:
-            self._Wy = self.mesh.cellGrady
+            a = self.mesh.r(self.mesh.area,'F','Fy','V')
+            self._Wy = sdiag(a)*self.mesh.cellGrady
         return self._Wy
 
     @property
     def Wz(self):
         if getattr(self, '_Wz', None) is None:
-            self._Wz = self.mesh.cellGradz
+            a = self.mesh.r(self.mesh.area,'F','Fz','V')
+            self._Wz = sdiag(a)*self.mesh.cellGradz
         return self._Wz
 
-    @property
-    def Ws(self):
-        if getattr(self,'_Ws', None) is None:
-            self._Ws = sdiag(self.mesh.vol)
-        return self._Ws
 
 
     def __init__(self, mesh):