added rejection and localization of samples

da/baseclasses/optimizer.py

@@ -45,9 +45,6 @@ class Optimizer(object):
         self.nmembers           = dims[1]
         self.nparams            = dims[2]
         self.nobs               = dims[3]
-        self.localization       = False
-        self.localization       = False
-
         self.CreateMatrices()
 
         return None
@@ -68,6 +65,11 @@ class Optimizer(object):
         self.obs             =  np.zeros( (self.nobs,), float)
         # covariance of observations
         self.R               =  np.zeros( (self.nobs,self.nobs,), float)
+        # localization of obs
+        self.may_localize       = np.zeros(self.nobs,int)
+        # rejection of obs
+        self.may_reject         = np.zeros(self.nobs,int)
+
         # Total covariance of fluxes and obs in units of obs [H P H^t + R]
         self.HPHR            =  np.zeros( (self.nobs,self.nobs,), float)
         # Kalman Gain matrix
@@ -91,6 +93,10 @@ class Optimizer(object):
 
             if members[0].ModelSample != None:
 
+                self.rejection_threshold                        = members[0].ModelSample.rejection_theshold
+                self.may_reject[:]                              = float(members[0].ModelSample.Data.getvalues('may_reject'))
+                self.may_localize[:]                            = float(members[0].ModelSample.Data.getvalues('may_localize'))
+
                 allobs.extend(                                    members[0].ModelSample.Data.getvalues('obs') )
                 allsamples.extend(                                members[0].ModelSample.Data.getvalues('simulated') )
                 allmdm.extend(                                    members[0].ModelSample.Data.getvalues('mdm')  )
@@ -286,21 +292,28 @@ class Optimizer(object):
         import numpy.linalg as la
 
         tvalue=1.97591
+
         for n in range(self.nobs):
 
-            if self.flags[n] == 1: continue
+            res                         = self.obs[n]-self.Hx[n]
+
+            # Screen for outliers greather than 3x model-data mismatch, only apply if obs may be rejected
+
+            if self.may_reject == 1:
+                if np.abs(res) > self.rejection_threshold*np.sqrt(self.R[n,n]):
+                    self.flags[n] == 1
+                    continue
 
             PHt                         = 1./(self.nmembers-1)*np.dot(self.X_prime,self.HX_prime[n,:])
             self.HPHR[n,n]              = 1./(self.nmembers-1)*(self.HX_prime[n,:]*self.HX_prime[n,:]).sum()+self.R[n,n]
 
             self.KG[:,n]                = PHt/self.HPHR[n,n]
 
-            dummy                       = self.Localize(n)
+            if self.may_localize[n] == 1:
+                dummy                   = self.Localize(n)
 
             alpha                       = np.double(1.0)/(np.double(1.0)+np.sqrt( (self.R[n,n])/self.HPHR[n,n] ) )
 
-            res                         = self.obs[n]-self.Hx[n]
-
             self.x[:]                   = self.x + self.KG[:,n]*res
 
             for r in range(self.nmembers):