Skip to content

Commit

Permalink
Update for alpha 5 release
Browse files Browse the repository at this point in the history
  • Loading branch information
JonKing93 committed Apr 13, 2021
1 parent a41aa58 commit 47ebe21
Show file tree
Hide file tree
Showing 79 changed files with 2,877 additions and 544 deletions.
59 changes: 47 additions & 12 deletions @PSM/PSM.m
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@
% estimate - Estimate proxy values from a state vector ensemble
% rename - Renames a PSM

% ----- Written By -----
% Jonathan King, University of Arizona, 2019-2020

properties (SetAccess = private)
name;
rows;
Expand Down Expand Up @@ -102,34 +105,64 @@
end

% Set the rows
function[obj] = useRows(obj, rows)
function[obj] = useRows(obj, rows, nRequired)
%% Specifies the state vector rows used by a PSM
%
% obj = obj.useRows(rows)
%
% obj = obj.useRows(rows, nRows)
% Require a specific number of rows.
%
% ----- Inputs -----
%
% rows: The rows used by the PSM. A vector of positive integers.
% rows: The state vector rows holding the data required to run
% the PSM to estimate observations. If a column vector, uses
% the same rows for every ensemble member. If a matrix, each
% column indicates the rows to use for the corresponding
% ensemble member. Either a logical matrix or matrix of
% linear indices.
%
% nRows: The number of rows required by the PSM. A scalar
% positive integer.
%
% ----- Outputs -----
%
% obj: The updated PSM object

% Default nRows
if ~exist('nRequired','var') || isempty(nRequired)
nRequired = [];
end

% Error check
assert(isvector(rows), 'rows must be a vector');
assert(~isempty(rows), 'rows cannot be empty');
assert(ndims(rows)<=3, 'rows cannot have more than 3 dimensions');
assert(islogical(rows) || isnumeric(rows), 'rows must be numeric or logical');

% Error check numeric indices
if isnumeric(rows)
dash.assertRealDefined(rows, 'rows');
dash.assertPositiveIntegers(rows, 'rows');

% Convert logical to numeric
% Error check logical indices
else
rows = find(rows);
nRows = unique(sum(rows, 1));
assert(isscalar(nRows), 'Each column of rows must select the same number of state vector elements (must have the same number of "true" elements).')
assert(nRows~=0, 'rows is not selecting any state vector rows because it has no "true" elements');

% Convert logical indices to linear
siz = size(rows, 1:3);
[rows, ~, ~] = ind2sub(siz, find(rows));
rows = reshape(rows, nRows, siz(2), siz(3));
end

% Optionally check the number of rows
if ~isempty(nRequired)
assert(size(rows,1)==nRequired, sprintf('You must select %.f rows', nRequired));
end

% Save
obj.rows = rows(:);
obj.rows = rows;
end

% Allow R to be estimated
Expand All @@ -152,16 +185,18 @@
end
end

% Download PSM repositories
methods (Static)
% Estimate Y values for a set of PSMs
[Ye, R] = estimate(X, F)

% Download the code for a PSM
download(psmName, path);

% Get the repository and commit information for a PSM
[repo, commit] = githubLocation(psmName);
end

% Estimate proxy values
methods (Static)
[Ye, R] = estimate(X, F, throwError); % User function call
F = setupEstimate(X, F); % Does error checking before computing estimates
[Ye, R] = computeEstimates(X, F, throwError);
end

% Run individual PSM objects
methods (Abstract)
Expand Down
152 changes: 152 additions & 0 deletions @PSM/computeEstimates.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
function[Ye, R] = computeEstimates(X, F, throwError)
%% Computes proxy estimates for an ensemble and set of PSMs. This is a
% low-level function and does not apply error-checking. Please see
% "PSM.estimate" for the recommended user command.
%
% Ye = PSM.computeEstimates(X, F)
% Estimate observations from the ensemble using the PSMs
%
% [Ye, R] = PSM.computeEstimates(X, F)
% Also estimate error-variances from PSMs when possible
%
% [...] = PSM.computeEstimates(X, F, throwError)
% Specify whether to throw an error when a PSM fails, or issue a warning.
%
% ----- Inputs -----
%
% X: A model ensemble. A numeric array with up to 3
% dimensions. (nState x nEns x nPrior)
%
% F: Cell vector of PSMs
%
% throwError: A scalar logical indicating whether to throw an error when a
% PSM fails (true), or instead issue a warning (false -- Default)
%
% ----- Outputs -----
%
% Ye: Proxy estimates for the prior(s). (nSite x nEns x nPrior)
%
% R: Error-variances for the priors. (nSite x nEns x nPrior)

% Default and error check error throws
if ~exist('throwError','var') || isempty(throwError)
throwError = false;
end
dash.assertScalarType(throwError, 'throwError', 'logical', 'logical');

% Sizes
nSite = numel(F);
[nState, nEns, nPrior] = size(X);

% Preallocate
Ye = NaN(nSite, nEns, nPrior);
if nargout>1
R = NaN(nSite, nEns, nPrior);
end

% Get the estimates for each PSM.
for s = 1:nSite
psm = F{s};

% Propagate rows over ensemble members and priors
[nRows, nCols, nDim3] = size(psm.rows);
if nCols==1
psm.rows = repmat(psm.rows, [1 nEns, 1]);
end
if nDim3==1
psm.rows = repmat(psm.rows, [1 1 nPrior]);
end

% Extract the data required to run the PSM
siz = [nRows, nEns, nPrior];
[~, c, p] = ind2sub(siz, (1:prod(siz))' );
index = sub2ind([nState, nEns, nPrior], psm.rows(:), c, p);
Xpsm = reshape( X(index), siz );

% Notify user if PSMs fail
ranPSM = false(1, nPrior);
psmFailed = false(1, nPrior);
Rfailed = false(1, nPrior);

% Run the PSM for each prior
for p = 1:nPrior
Xrun = Xpsm(:,:,p);

% Run the PSM with R estimation
if nargout>1 && psm.estimatesR
try
[Yrun, Rrun] = psm.runPSM(Xrun);
ranPSM(p) = true;
Rname = sprintf('R estimates for %s for prior %.f', psm.messageName(s), p);
dash.assertVectorTypeN(Rrun, 'numeric', nEns, Rname);
assert(isrow(Rrun), sprintf('%s must be a row vector', Rname));
catch ME
Rfailed(p) = true;

if throwError
rethrow(ME);
end
end
end

% Run the PSM without R estimation.
if ~ranPSM(p)
try
Yrun = psm.runPSM(Xrun);
ranPSM(p) = true;
catch ME
psmFailed(p) = true;
Rfailed(p) = true;

if throwError
rethrow(ME);
end
end
end

% Error check the Y output
if ranPSM(p)
try
Yname = sprintf("Y estimates for %s for prior %.f", psm.messageName(s), p);
dash.assertVectorTypeN(Yrun, 'numeric', nEns, Yname);
assert(isrow(Yrun), sprintf('%s must be a row vector', Yname));
catch ME
psmFailed(p) = true;
Rfailed(p) = true;

if throwError
rethrow(ME);
end
end
end

% Only record values if the PSM was successful
if ~psmFailed(p)
Ye(s,:,p) = Yrun;
if nargout>1 && psm.estimatesR && ~Rfailed(p)
R(s,:,p) = Rrun;
end
end
end

% Notify user of failed PSMs
if any(psmFailed)
whichPrior = '';
if ~all(psmFailed)
whichPrior = sprintf(' for prior(s) %s', dash.messageList(find(psmFailed)));
end
warning('%s failed to run%s', psm.messageName(s), whichPrior);
end

% Notify of failed R estimation
Rfailed = Rfailed & ~psmFailed;
if nargout>1 && psm.estimatesR && any(Rfailed)
whichPrior = '';
if ~all(Rfailed)
whichPrior = sprintf(' for prior(s) %s', dash.messageList(find(Rfailed)));
end
warning('%s failed to estimate R values%s', psm.messageName(s), whichPrior);
end
end

end
9 changes: 6 additions & 3 deletions @PSM/download.m
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,12 @@
% ----- Inputs -----
%
% psmName: The name of a PSM to download. A string. Options are
% 'bayspar': The BAYSPAR TEX86 PSM
% 'bayspline': The BAYSPLINE UK37 PSM
% 'baymag': The BAYMAG Mg/Ca PSM
% 'bayfox': The BayFOX Bayesian model for d18Oc of planktic foraminifera
% 'baymag': The BayMAG Bayesian model for Mg/Ca of planktic foraminifera
% 'bayspar': The BaySPAR Bayesian model for TEX86
% 'bayspline': The BaySPLINE Bayesian model for UK'37
% 'prysm': The PRySM python package of proxy system modelling tools
% 'vslite': The Vaganov-Shaskin Lite model of tree-ring width
%
% path: Indicates the folder where the PSM code should be downloaded. A
% string.
Expand Down
Loading

0 comments on commit 47ebe21

Please sign in to comment.