From e1ffa48e3ec189986258b3e2a3c5a602dd046ab0 Mon Sep 17 00:00:00 2001
From: Wouter Peters <wouter.peters@wur.nl>
Date: Fri, 30 Jul 2010 15:09:13 +0000
Subject: [PATCH] Added the least squares minimization algebra for serial
solutions, still to be tested
---
da/ct/optimizer.py | 31 +++++++++++++++++++++++++++++--
da/tools/general.py | 3 ++-
2 files changed, 31 insertions(+), 3 deletions(-)
diff --git a/da/ct/optimizer.py b/da/ct/optimizer.py
index da3d7204..375cf8d7 100755
--- a/da/ct/optimizer.py
+++ b/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
diff --git a/da/tools/general.py b/da/tools/general.py
index bd2cbf84..d07aaa67 100755
--- a/da/tools/general.py
+++ b/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
--
GitLab