diff --git a/SimPEG/Data.py b/SimPEG/Data.py
index 17eebaeb..6c94811a 100644
--- a/SimPEG/Data.py
+++ b/SimPEG/Data.py
@@ -7,7 +7,6 @@ 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
@@ -18,10 +17,32 @@ class BaseData(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    def setProblem(self, prob):
-        # Bind these two instances together using pointers
-        self.prob = prob
-        prob.data = self
+    @property
+    def prob(self):
+        """
+        The geophysical problem that explains this data, use::
+
+            data.setProblem(prob)
+        """
+        return getattr(self, '_prob', None)
+
+    def pair(self, p):
+        """Bind a problem to this data instance using pointers"""
+        assert hasattr(p, 'dataPair'), "Problem must have an attribute 'dataPair'."
+        assert isinstance(self, p.dataPair), "Problem requires data object must be an instance of a %s class."%(p.dataPair.__name__)
+        if p.ispaired:
+            raise Exception("The problem object is already paired to a data. Use prob.unpair()")
+        self._prob = p
+        p._data = self
+
+    def unpair(self):
+        """Unbind a problem from this data instance"""
+        if not self.ispaired: return
+        self.prob._data = None
+        self._prob = None
+
+    @property
+    def ispaired(self): return self.prob is not None
 
     @Utils.count
     @Utils.requires('prob')
@@ -117,6 +138,7 @@ class BaseData(object):
     def RHS(self, value):
         self._RHS = value
 
+    @property
     def isSynthetic(self):
         "Check if the data is synthetic."
         return (self.mtrue is not None)
diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py
index 65abb6eb..0bef2d92 100644
--- a/SimPEG/Inversion.py
+++ b/SimPEG/Inversion.py
@@ -2,6 +2,7 @@ import SimPEG
 from SimPEG import Utils, sp, np
 from Optimization import Remember, IterationPrinters, StoppingCriteria
 
+
 class BaseInversion(object):
     """BaseInversion(prob, reg, opt, data, **kwargs)
     """
diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index 18f7aba4..501b1c45 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -1,5 +1,6 @@
-from SimPEG import Utils, np, sp, Data
-norm = np.linalg.norm
+import Utils, Data
+import scipy.sparse as sp
+import numpy as np
 
 
 class BaseProblem(object):
@@ -41,12 +42,35 @@ class BaseProblem(object):
 
     dataPair = Data.BaseData
 
-
     def __init__(self, mesh, model, *args, **kwargs):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
         self.model = model
 
+    @property
+    def data(self):
+        """
+        The data object for this problem. Data
+        """
+        return getattr(self, '_data', None)
+
+    def pair(self, d):
+        """Bind a data to this problem instance using pointers."""
+        assert isinstance(d, self.dataPair), "Data object must be an instance of a %s class."%(self.dataPair.__name__)
+        if d.ispaired:
+            raise Exception("The data object is already paired to a problem. Use data.unpair()")
+        self._data = d
+        d._prob = self
+
+    def unpair(self):
+        """Unbind a data from this problem instance."""
+        if not self.ispaired: return
+        self.data._prob = None
+        self._data = None
+
+    @property
+    def ispaired(self): return self.data is not None
+
     @Utils.timeIt
     def J(self, m, v, u=None):
         """
@@ -143,8 +167,8 @@ class BaseProblem(object):
             and Wd which is the same size as data, and can be used to weight the inversion.
         """
         data = self.dataPair(mtrue=m, **geometry_kwargs)
-        data.setProblem(self)
-        data.dtrue = self.data.dpred(m,u=u)
+        data.pair(self)
+        data.dtrue = 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