diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index de7e0b6b..ec77b7a1 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -356,10 +356,29 @@ class update_lin_PreCond(InversionDirective):
             PC     = Utils.sdiag(diagA**-1.)
             self.opt.approxHinv = PC
             print 'Updated pre-cond'
+
+class update_Wj(InversionDirective):
+    """
+        Create approx-sensitivity base weighting
+    """
+    k = None
     
-#==============================================================================
-#          import pylab as plt
-#          plt.figure()
-#          ax = plt.subplot(221)
-#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
-#==============================================================================
+    def endIter(self):
+        
+        if self.opt.iter == 2:
+            m = self.invProb.curModel
+            if self.k is None:
+                self.k = int(self.survey.nD/10)
+                
+            def JtJv(v):
+
+                Jv = self.prob.Jvec(m, v)
+            
+                return self.prob.Jtvec(m,Jv)
+            
+            JtJdiag = Utils.diagEst(JtJv,6,k=100)
+            JtJdiag = JtJdiag
+            JtJdiag = JtJdiag / max(JtJdiag)
+        
+
+            self.reg.wght = JtJdiag
diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py
index 5f71d87c..f6b24f10 100644
--- a/SimPEG/Maps.py
+++ b/SimPEG/Maps.py
@@ -759,15 +759,27 @@ class PolyMap(IdentityMap):
 
             m = [\sigma_1, \sigma_2, c]
 
-    """
-    def __init__(self, mesh, order, logSigma=True, normal='X'):
+    """    
+    def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
         IdentityMap.__init__(self, mesh)
         self.logSigma = logSigma
         self.order = order
-        self.normal = normal
-
-    slope = 1e4
-
+        self.normal = normal  
+        self.actInd = actInd
+          
+        if getattr(self, 'actInd', None) is None:
+            self.actInd = range(self.mesh.nC)
+            self.nC = self.mesh.nC
+            
+        else:
+            self.nC = len(self.actInd)            
+            
+    slope = 1e4   
+    
+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+        
     @property
     def nP(self):
         if np.isscalar(self.order):
@@ -785,8 +797,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
             elif self.normal =='Y':
@@ -795,9 +807,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
+            Z = self.mesh.gridCC[self.actInd,2]
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
             elif self.normal =='Y':
@@ -806,6 +818,7 @@ class PolyMap(IdentityMap):
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
+                            
         else:
             raise(Exception("Only supports 2D"))
 
@@ -819,8 +832,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
 
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
@@ -832,9 +845,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
+            Z = self.mesh.gridCC[self.actInd,2]
 
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X