bug fix in residual calculation in serial_least_squares

da/stilt/optimizer.py

@@ -109,6 +109,8 @@ class CO2Optimizer(Optimizer):
                 logging.debug('Skipping observation (%s,%i) because of flag value %d' % (self.sitecode[n], self.obs_ids[n], self.flags[n]))
                 continue
 
+            res = self.obs[n] - self.Hx[n]
+
 
             logging.debug('Proceeding to assimilate observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
 
@@ -120,11 +122,11 @@ class CO2Optimizer(Optimizer):
             if 'surface' in self.sitecode[n]: 
                 self.KG[-1]=0.
                 self.KG[3078]=0
-                logging.debug("Surface observations do not update BC")
+                logging.debug('BC KG value set to zero for %s' %(self.sitecode[n]))
             if 'aircraft' in self.sitecode[n]: 
                 self.KG[0:3078]=0.
                 self.KG[3079:-1]=0.
-                logging.debug("Aircraft observations do not update Fluxes")
+                logging.debug('Flux KG values set to zero for %s' %(self.sitecode[n]))
 
             if self.may_localize[n]:
                 logging.debug('Trying to localize observation %s, %i' % (self.sitecode[n], self.obs_ids[n]))
@@ -134,7 +136,13 @@ class CO2Optimizer(Optimizer):
 
             alpha = np.double(1.0) / (np.double(1.0) + np.sqrt((self.R[n]) / self.HPHR[n]))
 
+
             self.x[:] = self.x + self.KG[:] * res
+            logging.debug('Residual %s'%res)
+            logging.debug('New self.KG BC1 %s' %self.KG[3078])
+            logging.debug('New self.KG BC2 %s' %self.KG[-1])
+            logging.debug('New self.x BC1 %s' %self.x[3078])
+            logging.debug('New self.x BC2 %s' %self.x[-1])
 
             for r in range(self.nmembers):
                 self.X_prime[:, r] = self.X_prime[:, r] - alpha * self.KG[:] * (self.HX_prime[n, r])