Added the least squares minimization algebra for serial solutions, still to be tested

da/ct/optimizer.py

@@ -77,8 +77,35 @@ class CtOptimizer():
 
         return None
 
-    def MinimumLeastSquares(self):
-        """ Make minimum least squares solution"""
+    def SerialMinimumLeastSquares(self):
+        """ Make minimum least squares solution by looping over obs"""
+        import numpy as np
+        import numpy .linalg as la
+
+        for n in range(self.nobs):
+
+            df_row                      = self.HX_prime[n,:]
+            QHt                         = 1./(self.nmembers-1)*np.dot(self.X_prime,df_row)
+            HQHR                        = 1./(self.nmembers-1)*(df_row*df_row).sum()+self.R[n,n]
+            K                           = QHt/HQHR
+
+            alpha                       = 1.0/(1.0+np.sqrt( (self.R[n,n])/HQHR ) )
+
+            res                         = self.obs[n]-self.Hx[n]
+
+            self.x[:]                   = self.x + K*res
+
+            for r in range(self.nmembers):
+                self.X_prime[:,r]       = self.X_prime[:,r]-alpha*K*(self.HX_prime[n,r])
+
+            for m in range(self.nobs):
+                res                     = self.obs[n]-self.Hx[n]
+                fac                     = 1.0/(self.nmembers-1)*(self.HX_prime[n,:]*self.HX_prime[m,:]).sum()/HQHR
+                self.Hx[m]              = self.Hx[m] + fac * res
+                self.HX_prime[m,:]      = self.HX_prime[m,:] - alpha*fac*self.HX_prime[n,:]
+
+    def BulkMinimumLeastSquares(self):
+        """ Make minimum least squares solution by solving matrix equations"""
         import numpy as np
         import numpy.linalg as la
 

da/tools/general.py

@@ -344,7 +344,8 @@ def Optimize(CycleInfo, StateVector):
 
     dummy     = optimizer.FillMatrices(StateVector )
 
-    dummy     = optimizer.MinimumLeastSquares()
+    dummy     = optimizer.BulkMinimumLeastSquares()
+    dummy     = optimizer.SerialMinimumLeastSquares()
 
     StateVector.isOptimized = True