Intial EnKF and ETKF commit

Author: Abhinab Bhattacharjee

src/+datools/+error/Error.m

@@ -1,4 +1,6 @@
 classdef Error < handle
+%ERROR base class for all kind of errors
+%   Curently supports Gaussian, Laplace, Logistics, Tent and Triangle.
 
     methods
 

src/+datools/+error/Gaussian.m

@@ -1,8 +1,10 @@
 classdef Gaussian < datools.error.Error
+%GAUSSIAN Gaussian child class for adding gaussina noise with a user
+%   defined deviation and bias
 
     properties
-        CovarianceSqrt
-        Bias
+        CovarianceSqrt    % Covariance(sqrt) for the noise
+        Bias              % Bias
     end
 
     methods

src/+datools/+examples/l63_enkf.m

@@ -30,7 +30,7 @@
 modelODE.TimeSpan = [0, Delta_t];
 
 % Propogate
-[tt, yy] = ode45(natureODE.Rhs.F, [0, 10], nature0);
+[tt, yy] = ode45(natureODE.RHS.F, [0, 10], nature0);
 natureODE.Y0 = yy(end, :).';
 
 % initialize model

src/+datools/+observation/Indexed.m renamed to src/+datools/+observation/+operator/Indexed.m

@@ -1,7 +1,8 @@
-classdef Indexed < datools.observation.Observation
+classdef Indexed < datools.observation.operator.Observation
+%INDEXED Defines the indexed observation operator
 
     properties
-        Indices
+        Indices    % Indices of observed state
     end
 
     methods
@@ -14,7 +15,7 @@
 
             s = p.Results;
 
-            [email protected](nvars, p.Unmatched);
+            [email protected].operator.Observation(nvars, p.Unmatched);
 
             obj.Indices = s.Indices;
 

src/+datools/+observation/Linear.m renamed to src/+datools/+observation/+operator/Linear.m

@@ -1,7 +1,9 @@
-classdef Linear < datools.observation.Observation
+classdef Linear < datools.observation.operator.Observation
+%LINEAR Defines the linear observation operator
+%   H is already linearized
 
     properties
-        H
+        H    % Observation operator
     end
 
     methods
@@ -14,7 +16,7 @@
 
             s = p.Results;
 
-            [email protected](nvars, p.Unmatched);
+            [email protected].operator.Observation(nvars, p.Unmatched);
 
             obj.H = s.H;
 
@@ -30,4 +32,4 @@
 
     end
 
-end
+end

src/+datools/+observation/Nonlinear.m renamed to src/+datools/+observation/+operator/Nonlinear.m

@@ -1,9 +1,10 @@
-classdef Nonlinear < datools.observation.Observation
+classdef Nonlinear < datools.observation.operator.Observation
+%NONLINEAR Defines the non linear observation operator
 
     properties
         F
-        %J
-        Indicies
+        Jacobian    % Jacobian of non-linear Observation operator
+        Indicies    % Indices of observed state
     end
 
     methods
@@ -12,16 +13,16 @@
             p = inputParser;
             p.KeepUnmatched = true;
             addParameter(p, 'F', @(~, x) x);
-            %addParameter(p, 'J', @(~, x) speye(numel(x)));
+            addParameter(p, 'Jacobian', @(~, x) speye(numel(x)));
             addParameter(p, 'Indicies', 1);
             parse(p, varargin{:});
 
             s = p.Results;
 
-            [email protected](nvars, p.Unmatched);
+            [email protected].operator.Observation(nvars, p.Unmatched);
 
             obj.F = s.F;
-            %obj.J = s.J;
+            obj.Jacobian = s.Jacobian;
             obj.Indicies = s.Indicies;
 
         end
@@ -32,8 +33,10 @@
         end
 
         function H = linearization(obj, t, x)
-            %H = obj.J(t, x);
-            H = 1;
+            H = obj.Jacobian(t, x);
+            if size(H, 1) == size(H, 2)
+                H = H(obj.Indices, :);
+            end
         end
 
     end

src/+datools/+observation/Observation.m renamed to src/+datools/+observation/+operator/Observation.m

@@ -1,8 +1,10 @@
 classdef Observation < handle
+%OBSERVATION Base class for every observation operator
+%   other details 
 
     properties
-        NumVars
-        ErrorModel
+        NumVars      % number of state variables
+        ErrorModel   % observation error model
     end
 
     methods

src/+datools/+observation/StateObservation.m

@@ -0,0 +1,49 @@
+classdef StateObservation < handle
+    %STATEOBSERVATION this is the class that handles the observed state
+    %   Observed state are noisy snapshots of reality (truth)
+
+    properties
+        Y         % Current Observation of States (maybe sparse)
+        R         % Observation Error Covariance
+        Noise     % Noise
+        NumObs    % Number of Observations(Make sure > 0)
+    end
+    
+    methods
+        function obj = StateObservation(varargin)
+            p = inputParser;
+            p.KeepUnmatched = true;
+            addParameter(p, 'Y', []);
+            addParameter(p, 'Noise', datools.error.Error);
+            validationFunction = @(x) (x>0);
+            addParameter(p, 'NumObs', 1, validationFunction);
+            
+            parse(p, varargin{:});
+
+            s = p.Results;
+            
+            
+            obj.Y = s.Y;
+            obj.Noise = s.Noise;
+            obj.NumObs = s.NumObs;
+            
+            unMatched = p.Unmatched;
+            
+            p = inputParser;
+            addParameter(p, 'R', speye(obj.NumObs));
+            
+            parse(p, unMatched);
+            s = p.Results;
+            obj.R = s.R;
+        end
+        
+        function updateY(obj, y)
+            obj.Y = y;
+        end
+        
+        function addNoise(obj, t, ~)
+            obj.Y = obj.Noise.adderr(t, obj.Y);
+        end
+    end
+
+end

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

@@ -2,33 +2,35 @@
 
     methods
 
-        function analysis(obj, R, y)
+        function analysis(obj)
 
             inflation = obj.Inflation;
 
             tc = obj.Model.TimeSpan(1);
 
             xf = obj.Ensemble;
             ensN = obj.NumEnsemble;
+            
+            R = obj.Observation.R;
 
             xfm = mean(xf, 2);
             Af = xf - repmat(xfm, 1, ensN);
             Af = inflation * Af/sqrt(ensN - 1);
             xf = repmat(xfm, 1, ensN) + Af*sqrt(ensN - 1);
 
-            Hxf = obj.Observation.observeWithoutError(tc, xf);
+            Hxf = obj.ObservationOperator.observeWithoutError(tc, xf);
             Hxfm = mean(Hxf, 2);
             HAf = Hxf - repmat(Hxfm, 1, ensN);
             HAf = HAf/sqrt(ensN - 1);
 
-            dS = decomposition((HAf * HAf.')+R, 'chol');
+            dS = decomposition((HAf * HAf.')+ R, 'chol');
             dR = decomposition(R, 'chol');
 
             temp = dS \ HAf;
             T = sqrtm(eye(ensN)-(HAf.' * temp));
 
             Aa = Af * T;
-            xam = xfm + ((Aa * (HAf * T).') * (dR \ (y - Hxfm)));
+            xam = xfm + ((Aa * (HAf * T).') * (dR \ (obj.Observation.Y - Hxfm)));
             xa = sqrt(ensN-1) .* Aa + repmat(xam, 1, ensN);
 
             obj.Ensemble = xa;

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

@@ -1,34 +1,35 @@
-classdef EnF < handle
-    % This is the base class for all statistical methods
-    % Derive from this class and implement methods/functions as required
-    % Deriving from handle base class allows an object of this class to be
-    % passed by reference.
+classdef EnF < datools.DABase
+%ENF This is the base class for all statistical methods
+%   Derive from this class and implement methods/functions as required
+%   Deriving from handle base class allows an object of this class to be
+%   passed by reference.
 
     properties
-        Model % type of ODE solver (ode45/Runge Kutta) and the model (eg: Lorenz63)
-        ModelError % type err
-        Observation % type of obervation
-        Ensemble % current ensemble values for all the states
-        Weights % Weight of each particle/ensemble
-        Inflation % inflation constant
-        Rejuvenation % rejuvenation constant
-        Localization % A boolean if localization needs to be used
-        Parallel % A boolean if parallel threads are to be implemented
-        RankHistogram % State variables for which RH is needed
-        RankValue % store the RH values for each state variables required
-        ResamplingThreshold % threshold below which resampling needs to be done
+        %Model                 % type of ODE solver and the model 
+        %ModelError            % type err
+        %Observation           % type of obervation
+        Ensemble               % current ensemble values for all the states
+        Anomalies              % Anomalies
+        Weights                % Weight of each particle/ensemble
+        Inflation              % inflation constant
+        Rejuvenation           % rejuvenation constant
+        Localization           % A boolean if localization needs to be used
+        Parallel               % A boolean if parallel threads are used
+        RankHistogram          % State variables for which RH is needed
+        RankValue              % store the RH values for each state variables required
+        ResamplingThreshold    % threshold below which resampling needs to be done
     end
 
     properties (Dependent)
-        BestEstimate % Current estimate of the particles/ensembles
-        NumEnsemble % Number of ensemble members
+        BestEstimate        % Current estimate of the particles/ensembles
+        NumEnsemble         % Number of ensemble members
     end
 
     methods (Abstract)
         % A method that will be implemented by child  classes to make
         % approximate inference on ensembles of states by combining
         % prior forecast/background data with noisy observations
-        analysis(obj, R, y)
+        analysis(obj)
     end
 
 
@@ -42,8 +43,9 @@
 
             p = inputParser;
             p.KeepUnmatched = true;
-            addRequired(p, 'Model', @(x) isa(x, 'datools.Model'));
-            addParameter(p, 'ModelError', datools.error.Error);
+            
+            %addRequired(p, 'Model', @(x) isa(x, 'datools.Model'));
+            %addParameter(p, 'ModelError', datools.error.Error);
             addParameter(p, 'NumEnsemble', 1);
             addParameter(p, 'Inflation', 1);
             addParameter(p, 'Rejuvenation', []);
@@ -52,34 +54,29 @@
             addParameter(p, 'RIPIterations', 0);
             addParameter(p, 'RankHistogram', [])
             addParameter(p, 'ResamplingThreshold', 0.5);
+            addParameter(p, 'EnsembleGenerator', @(x, y) randn(x, y));
+            
             parse(p, varargin{:});
 
             s = p.Results;
+            
+            [email protected](p.Unmatched);
 
-            obj.Model = s.Model;
-            obj.ModelError = s.ModelError;
             obj.Inflation = s.Inflation;
             obj.Rejuvenation = s.Rejuvenation;
             obj.Localization = s.Localization;
             obj.Parallel = s.Parallel;
             obj.RankHistogram = s.RankHistogram;
             obj.ResamplingThreshold = s.ResamplingThreshold;
             ensN = s.NumEnsemble;
+            
             if ~isempty(obj.RankHistogram)
                 RankValue = zeros(length(obj.RankHistogram), ensN+2);
             end
 
-            kept = p.Unmatched;
-
-            p = inputParser;
-            addParameter(p, 'Observation', datools.observation.Observation(s.Model.NumVars));
-            addParameter(p, 'EnsembleGenerator', @(x) randn(s.Model.NumVars, x));
-            parse(p, kept);
-
-            s = p.Results;
-
             obj.Ensemble = s.EnsembleGenerator(ensN);
-            obj.Observation = s.Observation;
+            obj.Anomalies = (1/sqrt(ensN)) * (obj.Ensemble - mean(obj.Ensemble, 2));
+            %obj.Observation = s.Observation;
             obj.Weights = ones(ensN, 1) / ensN;
 
         end
@@ -93,7 +90,7 @@ function forecast(obj)
             times = zeros(obj.NumEnsemble, 1);
 
             if obj.Parallel
-                rhs = obj.Model.ODEModel.RHS.F;
+                rhs = obj.Model.ODEModel.F;
                 tspan = obj.Model.TimeSpan;
                 solver = obj.Model.Solver;
                 ens = obj.Ensemble;
@@ -112,14 +109,14 @@ function forecast(obj)
                 end
 
                 for ensi = 1:ensN
-                    obj.Ensemble(:, ensi) = obj.ModelError.adderr(obj.Model.TimeSpan(end), ens(:, ensi));
+                    obj.Ensemble(:, ensi) = obj.ModelError.adderr(obj.Model.ODEModel.TimeSpan(end), ens(:, ensi));
                 end
 
             else
                 for ensi = 1:obj.NumEnsemble
                     [time, yend] = obj.Model.solve([], obj.Ensemble(:, ensi));
 
-                    obj.Ensemble(:, ensi) = obj.ModelError.adderr(obj.Model.TimeSpan(end), yend);
+                    obj.Ensemble(:, ensi) = obj.ModelError.adderr(obj.Model.ODEModel.TimeSpan(end), yend);
                     times(ensi) = time;
                 end
             end
@@ -162,7 +159,7 @@ function setMean(obj, xam)
             %   of the object of this class or a derived class to XAM
 
             X = obj.Ensemble;
-            ensN = size(X, 2);
+            ensN = size(X, 2);    % try obj.NumEnsemble
             xm = mean(X, 2);
             A = (X - repmat(xm, 1, ensN));
             X = repmat(xam, 1, ensN) + A;
@@ -188,7 +185,7 @@ function scaleAnomalies(obj, scale)
         function rejuvenate(obj, tau, Xf)
             %REJUVENATE  To reduce particle degeneracy
             %
-            %   REJUVENATE(OBJ, TAU, XF) adds a random noise in the form of
+            %   REJUVENATE(OBJ, TAU, XF) adds a noise in the form of
             %   random combination of background anomalies of the ensembles
             %   of states (using rejuvenation bandwidth TAU) to the
             %   transformation matrix. This matrix is used to rejuvenate