Added Kalman Filter and PSMs for alpha 4

@dash/dash.m

@@ -1,11 +1,17 @@
 classdef dash
     %% Contains various utility functions
 
-
     methods (Static)
 
-        % Analysis
+        % Maths
+        [Xmean, Xdev] = decompose(X, dim);
+        [C, Ycov] = estimateCovariance(X, Y);
         dist = haversine(latlon1, latlon2);
+        Y = gaspariCohn2D(X, R, scale);
+
+        % Localization                
+        [wloc, yloc] = localizationWeights(type, varargin);
+        [wloc, yloc] = gc2dLocalization(ensCoords, siteCoords, R, scale);
 
         % Misc
         [names, lon, lat, coord, lev, time, run, var]  = dimensionNames;

@dash/decompose.m

@@ -0,0 +1,41 @@
+function[Mmean, Mdev] = decompose(M, dim)
+%% Decomposes an ensemble into mean and deviations
+%
+% Mmean = decompose(M)
+% Returns the ensemble mean.
+%
+% [Mmean, Mdev] = decompose(M)
+% Also returns the ensemble deviations.
+%
+% [...] = decompose(M, dim)
+% Specify which dimension is the ensemble dimension. Default is the second
+% dimension.
+%
+% ----- Inputs -----
+%
+% M: An ensemble. A numeric array. Default size: (nState x nEns)
+%
+% dim: A scalar integer. Specifies which dimension is the ensemble
+%    dimension. By default, the second dimension is used as the ensemble
+%    dimension.
+%
+% ----- Outputs -----
+%
+% Mmean: The ensemble mean. Default size: (nState x 1)
+%
+% Mdev: The ensemble deviations. Default size: (nState x nEns)
+
+% Default ensemble dimension
+if ~exist('dim','var') || isempty(dim)
+    dim = 2;
+end
+
+% Ensemble mean
+Mmean = mean(M, dim);
+
+% Deviations
+if nargout>1
+    Mdev = M - Mmean;
+end
+
+end

@dash/estimateCovariance.m

@@ -0,0 +1,37 @@
+function[C, Ycov] = estimateCovariance(X, Y)
+%% Estimates covariance for a state vector ensemble and proxy estimates
+%
+% [C, Ycov] = dash.estimateCovariance(X, Y)
+%
+% ----- Inputs -----
+%
+% X: A state vector ensemble. A numeric matrix with dimensions
+%    (State vector x Ensemble members). May not contain NaN or Inf.
+%
+% Y: Model estimates of proxies/observations. A numeric matrix of size 
+%    (Proxy sites x Ensemble members). Must have the same number of columns
+%     as X. May not contain NaN or Inf.
+%
+% ----- Outputs -----
+%
+% C: The covariance of the state vector elements with the proxy sites. A
+%    numeric matrix of size (State vector x Proxy sites)
+%
+% Ycov: The covariance of the proxy sites with one another. A symmetric
+%    matrix of size (Proxy sites x Proxy sites)
+
+% Error check the inputs
+[~, nEns] = kalmanFilter.checkInput(X, 'X', false, true);
+[~, nEns2] = kalmanFilter.checkInput(Y, 'Y', false, true);
+assert(nEns==nEns2, 'Y and X must have the same number of columns');
+
+% Decompose
+[~, Xdev] = dash.decompose(X, 2);
+[~, Ydev] = dash.decompsoe(Y, 2);
+
+% Estimate covariances
+unbias = 1/(nEns-1);
+C = unbias .* (Xdev * Ydev');
+Ycov = unbias .* (Ydev * Ydev');
+
+end

@dash/gaspariCohn2D.m

@@ -0,0 +1,72 @@
+function[Y] = gaspariCohn2D(X, R, scale)
+%% Applies a Gaspari-Cohn 5th order polynomial in 2 dimensions for a
+% specified cutoff radius and length scale.
+%
+% weights = gaspariCohn2D(dist, R)
+% Determines localization weights based on a cutoff radius. Uses a length
+% scale of 0.5
+%
+% weights = gaspariCohn2D(dist, R, scale)
+% Sets the length scale for the Gaspari-Cohn polynomial.
+%
+% weights = gaspariCohn2D(dist, R, 'optimal')
+% Uses a length scale of sqrt(10/3) based on Lorenc (2003).
+%
+% ----- Inputs -----
+%
+% X: An numeric array. Cannot contain negative values.
+%
+% R: The cutoff radius. A positive scalar. Should use the same units as X.
+%
+% scale: A scalar on the interval (0 0.5]. Used to determine the cutoff
+%    radius. If unspecified, scale is set to 0.5 and the localization
+%    radius is equal to R. If less than 0.5, the cutoff radius will be
+%    smaller than R.
+%
+% ----- Outputs -----
+%
+% weights: A vector of localization weights
+
+% Set defaults
+if ~exist('scale','var') || isempty(scale)
+    scale = 0.5;
+elseif dash.isstrflag(scale) && strcmp(scale, 'optimal')
+    scale = sqrt(10/3);
+end
+
+% Error check
+dash.assertScalarType(R, 'R', 'numeric', 'numeric');
+dash.assertRealDefined(R, 'R', false, true);
+assert(R>0, 'R must be positive');
+
+dash.assertScalarType(scale, 'scale', 'numeric', 'numeric');
+dash.assertRealDefined(scale, 'scale');
+assert(scale>=0&scale<=0.5, 'scale must be on the interval [0 0.5]');
+
+assert(isnumeric(X), 'dist must be numeric');
+dash.assertRealDefined(X, 'dist', true, true);
+assert(~any(X(:)<0), 'dist cannot have negative elements');
+
+% Get the length scale and localization radius
+c = scale * R;
+Rloc = 2*c; % Note that Rloc <= R, they are not strictly equal
+
+% Get points inside/outside the radius.
+outRloc = (X > Rloc);
+inScale = (X <= c);
+outScale = (X > c) & (X <= Rloc);
+
+% Preallocate weights
+Y = ones( size(X) );
+
+% Apply the polynomial to the distances
+x = X / c;
+Y(inScale) = polyval([-.25,.5,.625,-5/3,0,1], x(inScale));
+Y(outScale) = polyval([1/12,-.5,.625,5/3,-5,4], x(outScale)) - 2./(3*x(outScale));
+Y(outRloc) = 0;
+
+% Weights should never be negative. Remove near-zero negative weights
+% resulting from round-off errors.
+Y( Y<0 ) = 0;
+
+end

@dash/gc2dLocalization.m

@@ -0,0 +1,58 @@
+function[wloc, yloc] = gc2dLocalization(ensCoords, siteCoords, R, scale)
+%% Determines localization weights using a Gaspari-Cohn 5th order polynomial
+% in 2 dimensions.
+%
+% [wloc, yloc] = gc2dLocalization(ensCoords, siteCoords, R)
+% Determines localization weights for a specified cutoff radius.
+%
+% [wloc, yloc] = gc2dLocalization(ensCoords, siteCoords, R, scale)
+% Also sets the length scale for the polynomial.
+%
+% ----- Inputs -----
+%
+% ensCoords: Latitude-longitude coordinates for an ensemble. A numeric
+%    matrix with two columns. The first column contains latitude
+%    coordinates and the second column is longitude.
+%
+% siteCoords: Latitude-longitude coordinates for proxy/observation sites. A
+%    numeric matrix with two columns. The first column hold latitude
+%    coordinates and the second column is longitude.
+%
+% R: The localization radius in kilometers. A positive scalar.
+%
+% scale: The length scale for the Gaspari-Cohn polynomial. A scalar on the
+%    interval (0, 0.5]. By default, the length scale is set to 0.5, which
+%    sets the covariance cutoff radius equal to R.
+%
+% ----- Outputs -----
+%
+% wloc: Localization weights between each state vector element (rows) and
+%    the proxy sites (columns). A numeric matrix.
+%
+% yloc: Localization weights between each proxy site and the other proxy
+%    sites. A symmetric numeric matrix.
+
+% Default for unset scale
+if ~exist('scale','var')
+    scale = [];
+end
+
+% Error check the coordinates
+assert(isnumeric(ensCoords), 'ensCoords must be numeric');
+assert(isnumeric(siteCoords), 'siteCoords must be numeric');
+dash.assertRealDefined(ensCoords, 'ensCoords', true, true);
+dash.assertRealDefined(siteCoords, 'siteCoords', true, true);
+assert(ismatrix(ensCoords), 'ensCoords must be a matrix');
+assert(ismatrix(siteCoords), 'siteCoords must be a matrix');
+assert(size(ensCoords,2)==2, 'ensCoords must have 2 columns');
+assert(size(siteCoords,2)==2, 'siteCoords must have 2 columns');
+
+% Get the distances
+wdist = dash.haversine(ensCoords, siteCoords);
+ydist = dash.haversine(siteCoords);
+
+% Apply the Gaspari-Cohn polynomial
+wloc = dash.gaspariCohn2D(wdist, R, scale);
+yloc = dash.gaspariCohn2D(ydist, R, scale);
+
+end

@dash/haversine.m

@@ -18,8 +18,8 @@
 %
 % dist: The distance between the coordinates. 
 %
-%    If you provided a single input, dist is a symmetric matrix with each
-%    rows and columns corresponding to the coordinates listed in latlon.
+%    If you provided a single input, dist is a symmetric matrix. The Nth 
+%    row/column corresponds to the Nth row of latlon.
 %
 %    If you provided two inputs, the rows of dist correspond to latlon1 and
 %    the columns correspond to latlon2.

@dash/localizationWeights.m

@@ -0,0 +1,47 @@
+function[wloc, yloc] = localizationWeights(type, varargin)
+%% Returns localization weights
+%
+% [wloc, yloc] = localizationWeights('gc2d', ensCoords, siteCoords, R)
+% [...] = localizationWeights('gc2d', ensCoords, siteCoords, R, scale)
+% Calculates localization weights using a Gaspari-Cohn 5th order polynomial
+% over 2 dimensions. Calculates weights using horizontal (latitude - 
+% longitude) distances.
+%
+% ----- Inputs -----
+%
+% ensCoords: Latitude-longitude coordinates for an ensemble. A numeric
+%    matrix with two columns. The first column contains latitude
+%    coordinates and the second column is longitude.
+%
+% siteCoords: Latitude-longitude coordinates for proxy/observation sites. A
+%    numeric matrix with two columns. The first column hold latitude
+%    coordinates and the second column is longitude.
+%
+% R: The localization radius in kilometers. A positive scalar.
+%
+% scale: The length scale for the Gaspari-Cohn polynomial. A scalar on the
+%    interval (0, 0.5]. By default, the length scale is set to 0.5, which
+%    sets the covariance cutoff radius equal to R.
+%
+% ----- Outputs -----
+%
+% wloc: Localization weights between each state vector element (rows) and
+%    the proxy sites (columns). A numeric matrix.
+%
+% yloc: Localization weights between each proxy site and the other proxy
+%    sites. A symmetric numeric matrix.
+
+% Error check the type
+dash.assertStrFlag(type, 'The first input');
+allowedTypes = "gc2d";
+assert( any(strcmpi(type, allowedTypes)), ...
+    sprintf(['The first input must be a recognized localization scheme. ',...
+    'Recognized types are: %s.'], dash.messageList(allowedTypes)) );
+
+% Switch to the appropriate localization scheme
+if strcmpi(type, 'gc2d')
+    [wloc, yloc] = dash.gc2dLocalization(varargin{:});
+end
+
+end
+

@ensembleMetadata/remove.m

@@ -13,21 +13,23 @@
 % obj: The updated ensembleMetadata object
 
 % Error check
-v = dash.checkStrsInList(varNames, obj.variableNames, 'varNames', 'variable in the state vector');
+if ~isempty(varNames)
+    v = dash.checkStrsInList(varNames, obj.variableNames, 'varNames', 'variable in the state vector');
 
-% Remove from arrays
-obj.variableNames(v) = [];
-obj.varLimit(v,:) = [];
-obj.nEls(v) = [];
-obj.dims(v) = [];
-obj.stateSize(v) = [];
-obj.isState(v) = [];
-obj.meanSize(v) = [];
-obj.metadata = rmfield(obj.metadata, varNames);
+    % Remove from arrays
+    obj.variableNames(v) = [];
+    obj.varLimit(v,:) = [];
+    obj.nEls(v) = [];
+    obj.dims(v) = [];
+    obj.stateSize(v) = [];
+    obj.isState(v) = [];
+    obj.meanSize(v) = [];
+    obj.metadata = rmfield(obj.metadata, varNames);
 
-% Update variable limits
-last = cumsum(obj.nEls);
-obj.varLimit = [last-obj.nEls+1, last];
+    % Update variable limits
+    last = cumsum(obj.nEls);
+    obj.varLimit = [last-obj.nEls+1, last];
+end
 
 end
 

@kalmanFilter/assertEditableCovariance.m

@@ -0,0 +1,23 @@
+function[] = assertEditableCovariance(kf, type)
+%% Checks that a prior and observations have been specified before
+% setting covariance options.
+%
+% kf.assertEditableCovariance(type)
+%
+% ----- Inputs -----
+%
+% type: A string listing the type of covariance modification being made.
+%    Use for error messages.
+
+% Require a prior
+if isempty(kf.M)
+    error(['You must specify a prior (using the "prior" command) before ',...
+        'setting %s.'], type);
+
+% Require observations
+elseif isempty(kf.D)
+    error(['You must specify the observations/proxies (using the ',...
+        '"observations" command) before setting %s.'], type);
+end
+
+end