Merge branch 'master' into abhinab

Author: Abhinab Bhattacharjee

.zenodo.json

@@ -0,0 +1,34 @@
+{
+    "description": "<p>Initial release</p>", 
+    "license": "other-open", 
+    "title": "ComputationalScienceLaboratory/DATools: v0.0.1", 
+    "version": "v0.0.1", 
+    "upload_type": "software", 
+    "publication_date": "2022-06-16", 
+    "creators": [
+        {
+            "affiliation": "Virginia Tech", 
+            "name": "Andrey A Popov"
+        }, 
+        {
+            "affiliation": "Virginia Tech", 
+            "name": "Abhinab Bhattacharjee"
+        }, 
+        {
+            "affiliation": "Virginia Tech", 
+            "name": "Amit N Subrahmanya"
+        }, 
+        {
+            "affiliation": "Virginia Tech", 
+            "name": "Adrian Sandu"
+        }
+    ], 
+    "access_right": "open", 
+    "related_identifiers": [
+        {
+            "scheme": "url", 
+            "identifier": "https://github.com/ComputationalScienceLaboratory/DATools/tree/v0.0.1", 
+            "relation": "isSupplementTo"
+        }
+    ]
+}

README.md

@@ -1,3 +1,5 @@
+[![DOI](https://zenodo.org/badge/269520040.svg)](https://zenodo.org/badge/latestdoi/269520040)
+
 A set of Data assimilation tools. Filter smoothers, ensemble and variational methods, localization, etc.
 
-An active work in progress.
+An active work in progress.

src/+datools/+statistical/+ensemble/ETKF.m

@@ -13,12 +13,13 @@ function analysis(obj, R, y)
 
             xfm = mean(xf, 2);
             Af = xf - repmat(xfm, 1, ensN);
-            Af = inflation * Af ./ sqrt(ensN-1);
-            xf = repmat(xfm, 1, ensN) + Af;
+            Af = inflation * Af/sqrt(ensN - 1);
+            xf = repmat(xfm, 1, ensN) + Af*sqrt(ensN - 1);
 
             Hxf = obj.Observation.observeWithoutError(tc, xf);
             Hxfm = mean(Hxf, 2);
             HAf = Hxf - repmat(Hxfm, 1, ensN);
+            HAf = HAf/sqrt(ensN - 1);
 
             dS = decomposition((HAf * HAf.')+R, 'chol');
             dR = decomposition(R, 'chol');
@@ -31,7 +32,6 @@ function analysis(obj, R, y)
             xa = sqrt(ensN-1) .* Aa + repmat(xam, 1, ensN);
 
             obj.Ensemble = xa;
-            %obj.Model.update(0, obj.BestEstimate);
 
             obj.Weights = ones(ensN, 1) / ensN;
 

src/+datools/+statistical/+ensemble/EnGMF.m

@@ -0,0 +1,55 @@
+classdef EnGMF < datools.statistical.ensemble.EnF
+
+    methods
+
+        function analysis(obj, R, y)
+
+            tc = obj.Model.TimeSpan(1);
+
+            xf = obj.Ensemble;
+            ensN = obj.NumEnsemble;
+            n = size(xf, 1);
+
+            beta = ((4/(n + 2))^(2/(n + 4)))*((ensN)^(-2/(n + 4)));
+
+            xfm = mean(xf, 2);
+            Af = xf - repmat(xfm, 1, ensN);
+            Af = sqrt(beta) * Af/sqrt(ensN - 1);
+
+            Hxf = obj.Observation.observeWithoutError(tc, xf);
+            Hxfm = mean(Hxf, 2);
+            HAf = Hxf - repmat(Hxfm, 1, ensN);
+            HAf = sqrt(beta)*HAf/sqrt(ensN - 1);
+
+            dS = decomposition((HAf * HAf.') + R, 'chol');
+
+            xtilde = xf - (Af*HAf.')*(dS\(Hxf - y));
+            Aa = xtilde - mean(xtilde, 2);
+            Aa = sqrt(beta)*Aa/sqrt(ensN - 1);
+
+            t0 = (Hxf - y);
+
+            as = (-0.5 * sum(t0.*(dS \ t0), 1)).';
+            m = max(as);
+            w = exp(as-(m + log(sum(exp(as-m))))).';
+            w = w/sum(w);
+
+            defensivefactor = 1;
+            w = defensivefactor*w + (1 - defensivefactor)/ensN;
+
+            % resample
+            xa = xf;
+            for  k = 1:ensN
+                ind = find(rand <= cumsum(w), 1, 'first');
+                xa(:, k) = xtilde(:, ind) + Aa*randn(ensN, 1);
+            end
+
+            obj.Ensemble = xa;
+
+            obj.Weights = ones(ensN, 1) / ensN;
+
+        end
+
+    end
+
+end