diff --git a/SimPEG/Data.py b/SimPEG/Data.py
index bc874de2..dc9239e1 100644
--- a/SimPEG/Data.py
+++ b/SimPEG/Data.py
@@ -1,5 +1,5 @@
 import Utils
-
+import numpy as np
 
 def requiresProblem(f):
     """
@@ -13,12 +13,13 @@ def requiresProblem(f):
     If a problem is not bound an Exception will be raised, and an nice error message printed.
     """
     extra = """
-        This function requires that a problem be bound to the data.
+        To use data.%s(), SimPEG requires that a problem be bound to the data.
         If a problem has not been bound, an Exception will be raised.
         To bind a problem to the Data object::
 
             data.setProblem(myProblem)
-    """
+
+    """ % f.__name__
     from functools import wraps
     @wraps(f)
     def requiresProblemWrapper(self,*args,**kwargs):
@@ -37,11 +38,11 @@ class BaseData(object):
 
     __metaclass__ = Utils.Save.Savable
 
+    prob = None      #: The geophysical problem that explains this data, use data.setProblem(prob)
     std = None       #: Estimated Standard Deviations
     dobs = None      #: Observed data
     dtrue = None     #: True data, if data is synthetic
     mtrue = None     #: True model, if data is synthetic
-    prob = None      #: The geophysical problem that explains this data
 
     counter = None   #: A SimPEG.Utils.Counter object
 
@@ -49,19 +50,39 @@ class BaseData(object):
         Utils.setKwargs(self, **kwargs)
 
     def setProblem(self, prob):
+        # Bind these two instances together using pointers
         self.prob = prob
+        prob.data = self
 
     @Utils.count
     @requiresProblem
     def dpred(self, m, u=None):
         """
-            Projection matrix.
+            Create the projected data from a model.
+            The field, u, (if provided) will be used for the predicted data
+            instead of recalculating the fields (which may be expensive!).
+
+            .. math::
+                d_\\text{pred} = P(u(m))
+
+            Where P is a projection of the fields onto the data space.
+        """
+        if u is None: u = self.prob.field(m)
+        return self.projectField(u)
+
+
+    @Utils.count
+    def projectField(self, u):
+        """
+            This function projects the fields onto the data space.
+
 
             .. math::
                 d_\\text{pred} = Pu(m)
         """
-        if u is None: u = self.prob.field(m)
-        return self.P*u
+        return u
+
+    #TODO: def projectFieldDeriv(self, u):  Does this need to be made??!
 
     @Utils.count
     def residual(self, m, u=None):
@@ -122,7 +143,7 @@ class BaseData(object):
         """
             Source matrix.
         """
-        return self._RHS
+        return getattr(self, '_RHS', None)
     @RHS.setter
     def RHS(self, value):
         self._RHS = value
diff --git a/SimPEG/Examples/DC.py b/SimPEG/Examples/DC.py
index fa50cd1c..a05e3e3c 100644
--- a/SimPEG/Examples/DC.py
+++ b/SimPEG/Examples/DC.py
@@ -10,9 +10,10 @@ class DCData(Data.BaseData):
 
     """
 
-    def __init__(self, mesh, model, **kwargs):
-        problem.BaseProblem.__init__(self, mesh, model)
-        self.mesh.setCellGradBC('neumann')
+    P = None #: projection
+
+    def __init__(self, **kwargs):
+        Data.BaseData.__init__(self, **kwargs)
         Utils.setKwargs(self, **kwargs)
 
     def reshapeFields(self, u):
@@ -20,18 +21,14 @@ class DCData(Data.BaseData):
             u = u.reshape([-1, self.RHS.shape[1]], order='F')
         return u
 
-    def dpred(self, m, u=None):
+    def projectField(self, u):
         """
             Predicted data.
 
             .. math::
                 d_\\text{pred} = Pu(m)
         """
-        if u is None:
-            u = self.field(m)
-
         u = self.reshapeFields(u)
-
         return Utils.mkvc(self.P*u)
 
 
@@ -47,7 +44,7 @@ class DCProblem(Problem.BaseProblem):
     dataPair = DCData
 
     def __init__(self, mesh, model, **kwargs):
-        problem.BaseProblem.__init__(self, mesh, model)
+        Problem.BaseProblem.__init__(self, mesh, model)
         self.mesh.setCellGradBC('neumann')
         Utils.setKwargs(self, **kwargs)
 
@@ -75,7 +72,7 @@ class DCProblem(Problem.BaseProblem):
     def field(self, m):
         A = self.createMatrix(m)
         solve = Solver(A)
-        phi = solve.solve(self.RHS)
+        phi = solve.solve(self.data.RHS)
         return Utils.mkvc(phi)
 
     def J(self, m, v, u=None):
@@ -103,9 +100,9 @@ class DCProblem(Problem.BaseProblem):
         if u is None:
             u = self.field(m)
 
-        u = self.reshapeFields(u)
+        u = self.data.reshapeFields(u)
 
-        P = self.P
+        P = self.data.P
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
         A = self.createMatrix(m)
@@ -128,10 +125,10 @@ class DCProblem(Problem.BaseProblem):
         if u is None:
             u = self.field(m)
 
-        u = self.reshapeFields(u)
-        v = self.reshapeFields(v)
+        u = self.data.reshapeFields(u)
+        v = self.data.reshapeFields(v)
 
-        P = self.P
+        P = self.data.P
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
         A = self.createMatrix(m)
@@ -186,7 +183,7 @@ if __name__ == '__main__':
     # Create the mesh
     h1 = np.ones(20)
     h2 = np.ones(100)
-    M = mesh.TensorMesh([h1,h2])
+    M = Mesh.TensorMesh([h1,h2])
 
     # Create some parameters for the model
     sig1 = np.log(1)
@@ -198,7 +195,7 @@ if __name__ == '__main__':
     condVals = [sig1, sig2]
     mSynth = Utils.ModelBuilder.defineBlockConductivity(p0,p1,M.gridCC,condVals)
     plt.colorbar(M.plotImage(mSynth))
-    plt.show()
+    # plt.show()
 
     # Set up the projection
     nelec = 50
@@ -211,29 +208,29 @@ if __name__ == '__main__':
     q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
     P = Q.T
 
-    # Create some data
-    problem = DCProblem(M)
-    problem.P = P
-    problem.RHS = q
-    data = problem.createSyntheticData(mSynth, std=0.05)
+    model = Model.LogModel()
+    prob = DCProblem(M, model)
 
-    u = problem.field(mSynth)
-    u = problem.reshapeFields(u)
+    # Create some data
+    data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
+
+    u = prob.field(mSynth)
+    u = data.reshapeFields(u)
     M.plotImage(u[:,10])
     # plt.show()
 
-    # Now set up the problem to do some minimization
-    # problem.dobs = dobs
-    # problem.std = dobs*0 + 0.05
+    # Now set up the prob to do some minimization
+    # prob.dobs = dobs
+    # prob.std = dobs*0 + 0.05
     m0 = M.gridCC[:,0]*0+sig2
 
-    opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-    reg = inverse.Regularization(M)
-    inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4)
+    opt = Inverse.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+    reg = Inverse.Regularization(M)
+    inv = Inverse.Inversion(prob, reg, opt, data, beta0=1e4)
 
     # Check Derivative
     derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)]
-    tests.checkDerivative(derChk, mSynth)
+    # Tests.checkDerivative(derChk, mSynth)
 
 
 
diff --git a/SimPEG/Inverse/Inversion.py b/SimPEG/Inverse/Inversion.py
index 1679c789..2e7d7f50 100644
--- a/SimPEG/Inverse/Inversion.py
+++ b/SimPEG/Inverse/Inversion.py
@@ -201,7 +201,7 @@ class BaseInversion(object):
         phi_d = self.dataObj(m, u)
         phi_m = self.reg.modelObj(m)
 
-        self.dpred = self.prob.dpred(m, u=u)  # This is a cheap matrix vector calculation.
+        self.dpred = self.data.dpred(m, u=u)  # This is a cheap matrix vector calculation.
         self.phi_d = phi_d
         self.phi_m = phi_m
 
@@ -245,7 +245,7 @@ class BaseInversion(object):
             u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
         """
         # TODO: ensure that this is a data is vector and Wd is a matrix.
-        R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
+        R = self.data.residualWeighted(m, u=u)
         R = Utils.mkvc(R)
         return 0.5*np.vdot(R, R)
 
@@ -285,9 +285,9 @@ class BaseInversion(object):
         if u is None:
             u = self.prob.field(m)
 
-        R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
+        R = self.data.residualWeighted(m, u=u)
 
-        dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
+        dmisfit = self.prob.Jt(m, self.data.Wd * R, u=u)
 
         return dmisfit
 
@@ -330,11 +330,11 @@ class BaseInversion(object):
         if u is None:
             u = self.prob.field(m)
 
-        R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
+        R = self.data.residualWeighted(m, u=u)
 
         # TODO: abstract to different norms a little cleaner.
         #                                        \/ it goes here. in l2 it is the identity.
-        dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u)
+        dmisfit = self.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.prob.J_approx(m, v, u=u), u=u)
 
         return dmisfit
 
diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index e3e5fed7..18f7aba4 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -130,7 +130,7 @@ class BaseProblem(object):
         """
         pass
 
-    def createSyntheticData(self, m, std=0.05, u=None):
+    def createSyntheticData(self, m, std=0.05, u=None, **geometry_kwargs):
         """
             Create synthetic data given a model, and a standard deviation.
 
@@ -142,11 +142,13 @@ class BaseProblem(object):
             Returns the observed data with random Gaussian noise
             and Wd which is the same size as data, and can be used to weight the inversion.
         """
-        dtrue = self.dpred(m,u=u)
-        noise = std*abs(dtrue)*np.random.randn(*dtrue.shape)
-        dobs = dtrue+noise
-        stdev = dobs*0 + std
-        return self.dataPair(dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m)
+        data = self.dataPair(mtrue=m, **geometry_kwargs)
+        data.setProblem(self)
+        data.dtrue = self.data.dpred(m,u=u)
+        noise = std*abs(data.dtrue)*np.random.randn(*data.dtrue.shape)
+        data.dobs = data.dtrue+noise
+        data.std = data.dobs*0 + std
+        return data
 
 
 
diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py
index be172b6b..ec9b61e8 100644
--- a/SimPEG/Utils/ModelBuilder.py
+++ b/SimPEG/Utils/ModelBuilder.py
@@ -196,7 +196,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
 
 if __name__ == '__main__':
 
-    from SimPEG.mesh import TensorMesh
+    from SimPEG.Mesh import TensorMesh
     from matplotlib import pyplot as plt
 
     # Define the mesh
diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py
index fc1c1d0f..b48a856e 100644
--- a/SimPEG/Utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -7,6 +7,7 @@ from ipythonutils import easyAnimate as animate
 from Solver import Solver
 import Save
 import Geophysics
+import ModelBuilder
 
 import types
 import time