Rank Histogram(Observation)

Author: Abhinab Bhattacharjee

src/+datools/+examples/+sandu/l63_experiments.m

@@ -4,11 +4,11 @@
 % Use this to run Lorenz 63 experiments
 modelname = 'Lorenz63';
 % uncomment the filter you want to run
-%filtername = 'EnKF';
+filtername = 'EnKF';
 %filtername = 'ETKF';
 %filtername = 'ETPF';
 %filtername = 'SIR';
-filtername = 'RHF';
+%filtername = 'RHF';
 
 % oservation variance
 variance = 1;
@@ -29,7 +29,7 @@
 
 %% Remaining code%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % rank histogram for the 1st state only (for now)
-histvar = 1:1:1;
+histvar = 1:1:2;
 % decide the type of filter
 switch filtername
     case 'EnKF'
@@ -224,7 +224,7 @@
             y = filter.Observation.observeWithError(model.TimeSpan(1), xt);
 
             % Rank histogram (if needed)
-            datools.utils.stat.RH(filter, xt);
+            datools.utils.stat.RH(filter, xt, y);
 
             % analysis
             % try

src/+datools/+examples/+sandu/runexperiments2.m

@@ -33,6 +33,10 @@ function runexperiments2(user)
 %localization radius
 r = user.localizationradius;
 
+% the states for which the histograms are saved
+histvar = user.histvar;
+RHmeasure = user.RHmeasure;
+
 % plot indices
 rankhistogramplotindex = user.rankhistogramplotindex;
 rmseplotindex = user.rmseplotindex;
@@ -41,7 +45,7 @@ function runexperiments2(user)
 
 %% Remaining code%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % rank histogram for the 1st state only (for now)
-histvar = 1:1:1;
+%histvar = 1:1:1;
 
 % decide the type of filter
 switch filtername
@@ -59,6 +63,8 @@ function runexperiments2(user)
         filtertype = 'Particle';
     case 'SIR'
         filtertype = 'Particle';
+    case 'SIS_EnKF'
+        filtertype = 'Particle';
     case 'RHF'
         filtertype = 'Ensemble';
     case 'EnGMF'
@@ -161,7 +167,7 @@ function runexperiments2(user)
 model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
 nature = datools.Model('Solver', solvernature, 'ODEModel', natureODE);
 
-% Observation Model
+% Observable model
 naturetomodel = datools.observation.Linear(numel(nature0), 'H', ...
     speye(natureODE.NumVars));
 
@@ -226,7 +232,7 @@ function runexperiments2(user)
 
     ns = numsamples;
     sE = zeros(ns, 1);
-    rankvaluesample = zeros(histvar, ensN+1, numsamples);
+    rankvaluesample = zeros(numel(histvar), ensN+1, numsamples);
     rmstempvalsample = nan * ones(steps-spinup, numsamples);
 
     parfor sample = 1:ns
@@ -260,6 +266,8 @@ function runexperiments2(user)
                     localization = [];
                 case 'SIR'
                     localization = [];
+                case 'SIS_EnKF'
+                    localization = [];
                 case 'ETPF'
                     localization = [];
                 case 'ETPF2'
@@ -317,6 +325,16 @@ function runexperiments2(user)
                     'Parallel', false, ...
                     'RankHistogram', histvar, ...
                     'Rejuvenation', rejuvenation);
+            case 'SIS_EnKF'
+                filter = datools.statistical.ensemble.(filtername)(model, ...
+                    'Observation', observation, ...
+                    'NumEnsemble', ensN, ...
+                    'ModelError', modelerror, ...
+                    'EnsembleGenerator', ensembleGenerator, ...
+                    'Parallel', false, ...
+                    'RankHistogram', histvar, ...
+                    'Rejuvenation', rejuvenation, ...
+                    'Inflation', 1.05);
             case 'RHF'
                 filter = datools.statistical.ensemble.(filtername)(model, ...
                     'Observation', observation, ...
@@ -380,17 +398,18 @@ function runexperiments2(user)
                     filter.forecast();
                 end
             catch e
-                dofilter = false
+                dofilter = false;
             end
 
             % observe
-            xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
-            y = filter.Observation.observeWithError(model.TimeSpan(1), xt);
+            xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State); % H(x_true)
+            y = filter.Observation.observeWithError(model.TimeSpan(1), xt); % H(x_true) + noise
 
             % Rank histogram (if needed)
-            datools.utils.stat.RH(filter, xt);
-            rankvaluesample(:, :, sample) = filter.RankValue(1, 1:end-1);
+%             datools.utils.stat.RH(filter, xt, y, Hxa, RHmeasure);
+%             rankvaluesample(:, :, sample) = filter.RankValue(:, 1:end-1);
 
+            
             % analysis
             try
                 if dofilter
@@ -401,6 +420,15 @@ function runexperiments2(user)
             end
 
             xa = filter.BestEstimate;
+            xaensembles = filter.Ensemble;
+            
+            % observable
+            hxa = naturetomodel.observeWithoutError(model.TimeSpan(1), xaensembles); % H(x_analysis)
+            Hxa = filter.Observation.observeWithError(model.TimeSpan(1), hxa); % H(x_analysis) + noise
+            
+            % Rank histogram (if needed)
+            datools.utils.stat.RH(filter, xt, y, Hxa, RHmeasure);
+            rankvaluesample(:, :, sample) = filter.RankValue(:, 1:end-1);
 
             err = xt - xa;
 
@@ -446,13 +474,22 @@ function runexperiments2(user)
     if isnan(resE)
         resE = 1000;
     end
+    
+%     totalklval = 0;
+%     for klcounter = 1:numel(histvar)
+%         klvalue = datools.utils.stat.KLDivergence(rankvalue(klcounter,:), (1 / (ensN+1))*ones(1, ensN+1));
+%         totalklval = totalklval + klvalue;
+%     end
+    
+    rankvalue = sum(rankvalue,1);
 
     switch filtertype
         case 'Ensemble'
             rmses(ensNi, infi) = resE;
 
             [xs, pval, rhplotval(ensNi, infi)] = datools.utils.stat.KLDiv(rankvalue, ...
-                (1 / ensN)*ones(1, ensN+1));
+                (1 / (ensN+1))*ones(1, ensN+1));
+            %rhplotval(ensNi, infi) = totalklval;
         case 'Particle'
             rmses(ensNi, reji) = resE;
 

src/+datools/+examples/+sandu/runnerscript.m

@@ -24,39 +24,46 @@
 % filtername = 'LETKF';
 
 %% User inputs
-options.modelname = 'Lorenz63';
+options.modelname = 'QG';
 
-options.filtername = 'EnKF';
+options.filtername = 'ETKF';
+
+%options.ensNs = [8, 12, 16, 20, 24, 28, 32];%[25, 50, 75, 100];
+options.ensNs = [16, 32, 48, 64, 80, 96, 112];
+%options.ensNs = [32, 48, 64, 80, 96, 112, 128];
+%options.ensNs = [32, 64, 128, 256, 512, 1024, 2048];
 
-options.ensNs = [4, 8, 12, 16, 20, 24, 28];%[25, 50, 75, 100];
 
 options.infs = round(linspace(1, 1.10, 7), 2);    % [1, 1.025, 1.05, 1.075, 1.10];
 
 options.variance = 1; % observation variance
 
-options.observeindicies = 1:1:3;  % observation indices
+options.observeindicies = linspace(0,8001,150);  % observation indices
 
 options.rejs = round(logspace(-1.5, -0.25, 7), 2); % round(2*logspace(-1.5, -0.5, 6), 2)
 
-options.spinups = 500;
+options.spinups = 50;
 
 options.steps = 11 * options.spinups; % change as deemed fit
 
-options.Dt = 0.12; % 0.12(L63), 0.05(L96), 0.0109(QG)
+options.Dt = 0.0109; % 0.12(L63), 0.05(L96), 0.0109(QG)
 
 options.odesolver = 'ode45'; % ode45 , RK4
 
 options.localize = false; % set to true if localization is needed
 
-options.localizationradius = 4; % localization rdius
+options.localizationradius = 4; % localization radius
+
+options.ns = 10;
 
-options.ns = 20;
+options.histvar = 1:1:1;
+options.RHmeasure = 'Truth';
 
 %plotting parameters
 options.rankhistogramplotindex = 1:2:numel(options.ensNs);
 options.rmseplotindex = 1:2:numel(options.ensNs);
-options.rmseheatmapplotindex = 1:1:numel(options.ensNs);
-options.kldivergenceplotindex = 1:1:numel(options.ensNs);
+options.rmseheatmapplotindex = 1:2:numel(options.ensNs);
+options.kldivergenceplotindex = 1:2:numel(options.ensNs);
 
 %%
 % call the experiment file

src/+datools/+utils/+stat/KLDiv.m

@@ -10,8 +10,7 @@
 function [xs, pval, KLVal] = KLDiv(pdf1, pdf2)
 
 if numel(pdf1) ~= numel(pdf2)
-    fprintf('Both vectors must be of same length');
-    return;
+    error('Both vectors must be of same length');
 end
 
 % check if the distributions are normalized

src/+datools/+utils/+stat/RH.m

@@ -1,7 +1,7 @@
-function RH(model, xt)
+function RH(filter, xt, y, Hxa, observationvariable)
 %function to plot the rank histogram
 
-flag = isa(model, 'datools.statistical.ensemble.EnF');
+flag = isa(filter, 'datools.statistical.ensemble.EnF');
 if ~flag
     frpintf('Model should be ensemble type');
     exit;
@@ -10,25 +10,46 @@ function RH(model, xt)
     fprintf('need to pass the truth at this time step');
     exit;
 end
-var = length(model.RankHistogram);
+var = length(filter.RankHistogram);
+varobservtion = filter.Observation.Indices;
+xa = filter.Ensemble;
+[~, ensN] = size(xa);
 
-xf = model.Ensemble;
-[~, ensN] = size(xf);
-
-for i = 1:var
-    tempVar = model.RankHistogram(i);
-    model.RankValue(i, ensN+2) = tempVar;
-    ensFor = sort(xf(tempVar, :));
-
-    if xt(tempVar) < ensFor(1)
-        model.RankValue(i, 1) = model.RankValue(i, 1) + 1;
-    elseif xt(tempVar) > ensFor(length(ensFor))
-        model.RankValue(i, ensN+1) = model.RankValue(i, ensN+1) + 1;
-    else
-        index = find(ensFor > xt(tempVar), 1);
-        model.RankValue(i, index) = model.RankValue(i, index) + 1;
-    end
 
+switch observationvariable
+    case 'Truth'
+        for i = 1:var
+            tempVar = filter.RankHistogram(i);
+            filter.RankValue(i, ensN+2) = tempVar;
+            ensFor = sort(xa(tempVar, :));
+            
+            if xt(tempVar) < ensFor(1)
+                filter.RankValue(i, 1) = filter.RankValue(i, 1) + 1;
+            elseif xt(tempVar) > ensFor(length(ensFor))
+                filter.RankValue(i, ensN+1) = filter.RankValue(i, ensN+1) + 1;
+            else
+                index = find(ensFor > xt(tempVar), 1);
+                filter.RankValue(i, index) = filter.RankValue(i, index) + 1;
+            end
+            
+        end
+    case 'Observation'
+        for i = 1:numel(varobservtion)
+            tempVar = filter.RankHistogram(varobservtion(i));
+            filter.RankValue(i, ensN+2) = tempVar;
+            ensFor = sort(Hxa(i, :));
+            
+            if y(i) < ensFor(1)
+                filter.RankValue(i, 1) = filter.RankValue(i, 1) + 1;
+            elseif y(i) > ensFor(length(ensFor))
+                filter.RankValue(i, ensN+1) = filter.RankValue(i, ensN+1) + 1;
+            else
+                index = find(ensFor > y(i), 1);
+                filter.RankValue(i, index) = filter.RankValue(i, index) + 1;
+            end
+            
+        end
 end
 
+
 end