Support ndarray B in 3D-Var (ref #53 & email by Burba Mareike)

dapper/da_methods/baseline.py

@@ -95,13 +95,16 @@ class Var3D:
     def assimilate(self, HMM, xx, yy):
         Dyn, Obs, chrono, X0, stats = HMM.Dyn, HMM.Obs, HMM.t, HMM.X0, self.stats
 
-        if self.B in (None, 'clim'):
-            # Use climatological cov, ...
-            B = np.cov(xx.T)  # ... estimated from truth
+        if isinstance(self.B, np.ndarray):
+            # compare ndarray 1st to avoid == error for ndarray
+            B = self.B.astype(float)
+        elif self.B in (None, 'clim'):
+            # Use climatological cov, estimated from truth
+            B = np.cov(xx.T)
         elif self.B == 'eye':
             B = np.eye(HMM.Nx)
         else:
-            B = self.B
+            raise ValueError("Bad input B.")
         B *= self.xB
 
         # ONLY USED FOR DIAGNOSTICS, not to change the Kalman gain.

dapper/da_methods/variational.py

@@ -264,13 +264,16 @@ def assimilate(self, HMM, xx, yy):
         Nx = Dyn.M
 
         # Set background covariance. Note that it is static (compare to iEnKS).
-        if self.B in (None, 'clim'):
-            # Use climatological cov, ...
-            B = np.cov(xx.T)  # ... estimated from truth
+        if isinstance(self.B, np.ndarray):
+            # compare ndarray 1st to avoid == error for ndarray
+            B = self.B.astype(float)
+        elif self.B in (None, 'clim'):
+            # Use climatological cov, estimated from truth
+            B = np.cov(xx.T)
         elif self.B == 'eye':
             B = np.eye(Nx)
         else:
-            B = self.B
+            raise ValueError("Bad input B.")
         B *= self.xB
         B12 = CovMat(B).sym_sqrt
 

dapper/xp_launch.py

@@ -16,6 +16,7 @@
 from textwrap import dedent
 
 import dill
+import numpy as np
 import struct_tools
 import tabulate as _tabulate
 from tabulate import tabulate
@@ -316,9 +317,31 @@ def _aggregate_keys():
 
         for key in _aggregate_keys():
 
-            # Want to distinguish actual None's from empty ("N/A").
-            # => Don't use getattr(obj,key,None)
-            vals = [getattr(xp, key, "N/A") for xp in self]
+            def _getattr(xp, key):
+                # Don't use None, to avoid mixing with actual None's
+                # TODO 4: use an object yet more likely to be unique.
+                missing = "N/A"
+                a = getattr(xp, key, missing)
+
+                # Replace ndarray by its id, since o/w it will
+                # complain that you must use all().
+                # Alternative: replace all == (and !=) below by "is".
+                #     Tabulation with multi-line params actually works,
+                #     (though it's still likely to take up too much space,
+                #     unless we set np.printoptions...).
+                #     However, then python (since 3.8) will complain about
+                #     comparison to literal.
+                if isinstance(a, np.ndarray):
+                    shorten = 6
+                    a = f"arr(<id {id(a)//10**shorten}>)"
+                # TODO 3: leave formatting to sub() below?
+                # TODO 4: do similar formatting for other non-trivial params?
+                # TODO 4: document alternative way to specify non-trivial params:
+                #         use key to be looked up in some globally accessible dct.
+                #         Advantage: names are meaningful, rather than ids.
+                return a
+
+            vals = [_getattr(xp, key) for xp in self]
 
             # Sort (assign dct) into distinct, redundant, common
             if struct_tools.flexcomp(key, *nomerge):