Merge branch 'percentiles'

4. DA Analyses/@kalmanFilter/jointENSRF.m

@@ -1,7 +1,7 @@
-function[output] = jointENSRF( M, D, R, F, w, yloc, meanOnly, fullDevs, reconstruct )
+function[output] = jointENSRF( M, D, R, F, w, yloc, meanOnly, fullDevs, percentiles, reconstruct )
 %% Implements an ensemble square root kalman filter updating observations jointly
 %
-% [output] = dash.jointENSRF( M, D, R, F, w, yloc, meanOnly, reconstruct )
+% [output] = dash.jointENSRF( M, D, R, F, w, yloc, meanOnly, fullDevs, percentiles, reconstruct )
 %
 % ----- Inputs -----
 %
@@ -22,6 +22,9 @@
 %
 % fullDevs: Whether to return full ensemble deviations. Scalar logical
 %
+% percentiles: Which percentiles to return. A vector of values between 0
+%              and 100 (nPerc x 1)
+%
 % reconstruct: Logical vector indicating which state vector elements to
 %              reconstruct. (nState x 1)
 %
@@ -37,6 +40,8 @@
 %
 %   Adev - Updated ensemble deviations (nState x nEns x nTime)
 %
+%   Aperc - Percentiles of the updated ensemble. (nState x nPercentile x nTime)
+%
 %   Ye - Proxy estimates (nObs x nEns)
 %
 %   calibRatio - The calibration ratio. (nObs x nTime)
@@ -49,6 +54,7 @@
 
 % Get sizes
 [nObs, nTime] = size(D);
+nPerc = numel( percentiles );
 nEns = size(M,2);
 
 % Preallocate PSM outputs and calibration ratio
@@ -74,6 +80,7 @@
 elseif ~meanOnly
     Avar = NaN( nState, nTime );
 end
+Aperc = NaN( nState, nPerc, nTime );
 
 % Get (static) Kalman numerator. Clear M for space
 [Mmean, Mdev] = dash.decompose( M );
@@ -92,13 +99,20 @@
     Amean(:,t) = Mmean + K * ( D(obs,t) - Ymean(obs) );    
     calibRatio( obs, t ) = abs( D(obs,t) - Ymean(obs) ).^2 ./ ( diag(Kdenom) );
 
-    % Optionally update the deviations / variance
+    % Optionally update the deviations / variance / percentiles
     if ~meanOnly
         Ka = kalmanFilter.jointKalman( 'Ka', Knum(:,obs), Kdenom, R(obs,t) );
         if fullDevs
             Adev(:,:,t) = Mdev - Ka * Ydev(obs,:);
-        else
+            if nPerc > 0
+                Aperc(:,:,t) = Amean(:,t) + prctile( Adev(:,:,t), percentiles, 2 );
+            end
+        elseif nPerc == 0
             Avar(:,t) = sum(   (Mdev - Ka * Ydev(obs,:)).^2, 2) ./ (nEns-1);
+        else
+            Adev = Mdev - Ka * Ydev(obs,:);
+            Avar(:,t) = sum( Adev.^2, 2 ) ./ (nEns-1);
+            Aperc(:,:,t) = Amean(:,t) + prctile( Adev, percentiles, 2 );
         end
     end
 
@@ -117,6 +131,10 @@
 elseif ~meanOnly
     output.Avar = Avar;
 end
+if nPerc > 0
+    output.settings.percentiles = percentiles(:)';
+    output.Aperc = Aperc;
+end
 output.Ye = Ye;
 output.calibRatio = calibRatio;
 output.R = R;

4. DA Analyses/@kalmanFilter/kalmanFilter.m

@@ -16,6 +16,7 @@
         append;          % Whether to use the appended Ye method
         meanOnly;        % Whether to only calculate the ensemble mean
         fullDevs;        % Whether to return full ensemble deviations
+        percentiles;     % Which percentiles of the ensemble to return
         reconstruct;     % Which state vector elements to reconstruct
         reconH;          % Whether all H indices are reconstructed
     end
@@ -54,10 +55,11 @@
             obj.append = false;
             obj.meanOnly = false;
             obj.fullDevs = false;
+            obj.percentiles = [];
 
             % Block empty constructor, set values
             if isempty(M) || isempty(D) || isempty(R) || isempty(F)
-                error('M, D, R, and F not be empty.');
+                error('M, D, R, and F cannot be empty.');
             end
             obj.setValues( M, D, R, F );
 
@@ -90,10 +92,10 @@
     methods (Static)
 
         % Serial updating scheme
-        [output] = serialENSRF( M, D, R, F, w, fullDevs );
+        [output] = serialENSRF( M, D, R, F, w, fullDevs, percentiles );
 
         % Full inversion
-        [output] = jointENSRF( M, D, R, F, w, yloc, meanOnly, fullDevs, reconstruct );
+        [output] = jointENSRF( M, D, R, F, w, yloc, meanOnly, fullDevs, percentiles, reconstruct );
 
         % Serial kalman gain
         [K, a] = serialKalman( Mdev, Ydev, w, R );

4. DA Analyses/@kalmanFilter/run.m

@@ -15,6 +15,9 @@
 %
 %   Adev - Updated ensemble deviations (nState x nEns x nTime)
 %
+%   Aperc - Percentiles of the updated ensemble (nState x nPercentile x nTime)
+%           (See output.settings for the percentiles calculated) 
+%
 %   Ye - Proxy estimates
 %           Joint Updates:  (nObs x nEns)
 %           Serial Updates: (nObs x nEns x nTime)
@@ -69,7 +72,7 @@
     F = obj.adjustH( F, reconstruct );
 
     % Do the updates
-    output = obj.serialENSRF( M, obj.D, obj.R, F, w, obj.fullDevs );
+    output = obj.serialENSRF( M, obj.D, obj.R, F, w, obj.fullDevs, obj.percentiles );
 
     % Unappend if necessary
     output.Append = false;
@@ -92,7 +95,7 @@
     end
 
     % Do the updates
-    output = obj.jointENSRF( M, obj.D, obj.R, obj.F, w, yloc, obj.meanOnly, obj.fullDevs, reconstruct );
+    output = obj.jointENSRF( M, obj.D, obj.R, obj.F, w, yloc, obj.meanOnly, obj.fullDevs, obj.percentiles, reconstruct );
 end
 
 end

4. DA Analyses/@kalmanFilter/serialENSRF.m

@@ -1,4 +1,4 @@
-function[output] = serialENSRF( M, D, R, F, w, fullDevs )
+function[output] = serialENSRF( M, D, R, F, w, fullDevs, percentiles )
 %% Implements an ensemble square root kalman filter with serial updates.
 %
 % [output] = dash.serialENSRF( M, D, R, F, w )
@@ -16,6 +16,12 @@
 %
 % w: Covariance localization weights. (nState x nObs)
 %
+% percentiles: A vector of values between 0 and 100 specifying which
+%              percentiles to return. (nPerc x 1)
+%
+% fullDevs: Scalar logical indicating whether to return full ensemble
+%           deviations.
+%
 % ----- Outputs -----
 %
 % output: A structure with the following fields
@@ -24,8 +30,12 @@
 %
 %   Amean - The updated ensemble mean (nState x nTime)
 %
+%   Adev - Updated ensemble deviations (nState x nEns x nTime)
+%
 %   Avar - Updated ensemble variance (nState x nTime)
 %
+%   Aperc - Percentiles of the updated ensemble. (nState x nPerc x nTime)
+%
 %   Ye - Proxy estimates (nObs x nEns x nTime)
 %
 %   calibRatio - The calibration ratio. (nObs x nTime)
@@ -39,6 +49,7 @@
 % Get some sizes
 [nObs, nTime] = size(D);
 [nState, nEns] = size(M);
+nPerc = numel( percentiles );
 
 % Decompose the initial ensemble. Clear the ensemble to free memory.
 [Mmean, Mdev] = dash.decompose(M);
@@ -51,6 +62,7 @@
 else
     Avar = NaN( nState, nTime );
 end
+Aperc = NaN( nState, nPerc, nTime );
 Ye = NaN( nObs, nEns, nTime );
 sites = false( nObs, nTime );
 calibRatio = NaN( nObs, nTime );
@@ -90,6 +102,9 @@
     else
         Avar(:,t) = sum( Ad.^2, 2 ) ./ (nEns - 1);
     end
+    if nPerc > 0
+        Aperc(:,:,t) = Am + prctile( Ad, percentiles, 2 );
+    end
     progressbar(t/nTime);
 end
 
@@ -105,6 +120,10 @@
 else
     output.Avar = Avar;
 end
+if nPerc > 0
+    output.settings.percentiles = percentiles(:)';
+    output.Aperc = Aperc;
+end
 output.Ye = Ye;
 output.calibRatio = calibRatio;
 output.R = R;

4. DA Analyses/@kalmanFilter/settings.m

@@ -24,6 +24,9 @@
 % Specify whether to return full ensembles deviations, or just the
 % variance. Default is just the variance.
 %
+% obj.settings( ..., 'percentiles', percentiles )
+% Specify which percentiles of the ensemble to return
+%
 % obj.settings( ..., 'reconstruct', reconstruct )
 % Specify which state vector elements to reconstruct. Not recommended. See
 % "kalmanFilter.reconstructVars" instead.
@@ -50,14 +53,17 @@
 % fullDevs: A scalar logical indicating whether to return full ensemble
 %           deviations. Default is false.
 %
+% percentiles: A vector of values between 0 and 100, specifying which
+%              ensemble percentiles to return.
+%
 % reconstruct: A logical vector specifying which state vector elements to
 %              reconstruct. (nState x 1)
 
 % Parse inputs
-[type, weights, inflate, append, meanOnly, fullDevs, recon] = parseInputs( varargin, ...
-    {'type','localize','inflate','append','meanOnly','returnDevs','reconstruct'}, ...
-    {obj.type, obj.localize, obj.inflate, obj.append, obj.meanOnly, obj.fullDevs, obj.reconstruct}, ...
-    {[],[],[],[],[],[],[]} );
+[type, weights, inflate, append, meanOnly, fullDevs, percentiles, recon] = parseInputs( varargin, ...
+    {'type','localize','inflate','append','meanOnly','returnDevs','percentiles','reconstruct'}, ...
+    {obj.type, obj.localize, obj.inflate, obj.append, obj.meanOnly, obj.fullDevs, obj.percentiles, obj.reconstruct}, ...
+    {[],[],[],[],[],[],[],[]} );
 
 % Error checking
 if ~isstrflag(type)
@@ -128,6 +134,15 @@
     end
 end
 
+% Percentiles
+if ~isempty( percentiles )
+    if meanOnly
+        error('Cannot calculate ensemble percentiles when only updating the ensemble mean.');
+    elseif ~isnumeric(percentiles) || ~isvector(percentiles) || any(percentiles<0) || any(percentiles>100) || any(isnan(percentiles)) 
+        error('percentiles must be a vector of numeric values between 0 and 100 that do not contain NaN values.');
+    end
+end
+
 % Save values
 obj.type = type;
 obj.localize = weights;
@@ -136,6 +151,7 @@
 obj.meanOnly = meanOnly;
 obj.fullDevs = fullDevs;
 obj.reconstruct = reconstruct;
+obj.percentiles = percentiles;
 obj.reconH = reconH;
 
 end