From 47ebe21445257aab3dc9cad741135a79802ea4be Mon Sep 17 00:00:00 2001 From: JonKing93 Date: Tue, 13 Apr 2021 16:03:34 -0700 Subject: [PATCH] Update for alpha 5 release --- @PSM/PSM.m | 59 ++++++-- @PSM/computeEstimates.m | 152 +++++++++++++++++++ @PSM/download.m | 9 +- @PSM/estimate.m | 135 ++++------------- @PSM/githubLocation.m | 7 +- @PSM/setupEstimate.m | 52 +++++++ @dash/assertScalarType.m | 2 +- @dash/dash.m | 1 + @dash/version.m | 12 ++ @ensemble/ensemble.m | 6 +- @ensemble/loadSettings.m | 9 +- @ensembleFilter/Rcovariance.m | 45 ++++++ @ensembleFilter/checkInput.m | 48 ++++++ @ensembleFilter/checkMissingR.m | 31 ++++ @ensembleFilter/ensembleFilter.m | 54 +++++++ @ensembleFilter/estimates.m | 66 +++++++++ @ensembleFilter/finalize.m | 21 +++ @ensembleFilter/finalizeWhich.m | 10 ++ @ensembleFilter/observations.m | 37 +++++ @ensembleFilter/parseWhich.m | 54 +++++++ @ensembleFilter/prior.m | 69 +++++++++ @ensembleFilter/rename.m | 16 ++ @ensembleFilter/uncertainties.m | 101 +++++++++++++ @ensembleMetadata/buildFromPrimitives.m | 2 +- @ensembleMetadata/closestLatLon.m | 51 ++++++- @ensembleMetadata/ensembleMetadata.m | 2 +- @ensembleMetadata/sizes.m | 6 +- @evolvingEnsemble/evolvingEnsemble.m | 114 +++++++++++++++ @kalmanFilter/assertEditableCovariance.m | 4 +- @kalmanFilter/checkCovariance.m | 2 +- @kalmanFilter/covarianceEstimate.m | 12 +- @kalmanFilter/estimateCovariance.m | 8 +- @kalmanFilter/finalize.m | 23 ++- @kalmanFilter/index.m | 20 ++- @kalmanFilter/inflate.m | 2 +- @kalmanFilter/kalmanFilter.m | 43 ++---- @kalmanFilter/observations.m | 73 ++-------- @kalmanFilter/prior.m | 57 +++----- @kalmanFilter/run.m | 99 ++++++++++--- @optimalSensor/checkR.m | 9 ++ @optimalSensor/computeMetric.m | 13 ++ @optimalSensor/estimates.m | 38 +++++ @optimalSensor/meanMetric.m | 10 ++ @optimalSensor/metric.m | 37 +++++ @optimalSensor/optimalSensor.m | 53 +++++++ @optimalSensor/prior.m | 38 +++++ @optimalSensor/psms.m | 33 +++++ @optimalSensor/rename.m | 3 + @optimalSensor/run.m | 98 +++++++++++++ @optimalSensor/saveMeanArgs.m | 20 +++ @particleFilter/bayesWeights.m | 20 +++ @particleFilter/bestWeights.m | 25 ++++ @particleFilter/estimates.m | 30 ++++ @particleFilter/particleFilter.m | 69 +++++++++ @particleFilter/prior.m | 45 ++++++ @particleFilter/run.m | 66 +++++++++ @particleFilter/weighting.m | 42 ++++++ @particleFilter/weights.m | 22 +++ @stateVectorVariable/stateVectorVariable.m | 3 + README.md | 28 ++++ bayfoxPSM.m | 123 ++++++++++++++++ baymagPSM.m | 139 ++++++++++-------- baysparPSM.m | 146 +++++++++---------- baysplinePSM.m | 117 +++++++-------- dataSource.m | 3 + linearPSM.m | 14 +- matSource.m | 3 + ncSource.m | 3 + opendapSource.m | 3 + posteriorCalculation.m | 3 + posteriorIndex.m | 26 +++- posteriorPercentiles.m | 3 + posteriorVariance.m | 3 + progressbar.m | 3 + prysmCellulose.m | 139 ++++++++++++++++++ prysmCoral.m | 161 +++++++++++++++++++++ prysmIcecore.m | 105 ++++++++++++++ prysmSpeleothem.m | 131 +++++++++++++++++ vslitePSM.m | 80 +++++++++- 79 files changed, 2877 insertions(+), 544 deletions(-) create mode 100644 @PSM/computeEstimates.m create mode 100644 @PSM/setupEstimate.m create mode 100644 @dash/version.m create mode 100644 @ensembleFilter/Rcovariance.m create mode 100644 @ensembleFilter/checkInput.m create mode 100644 @ensembleFilter/checkMissingR.m create mode 100644 @ensembleFilter/ensembleFilter.m create mode 100644 @ensembleFilter/estimates.m create mode 100644 @ensembleFilter/finalize.m create mode 100644 @ensembleFilter/finalizeWhich.m create mode 100644 @ensembleFilter/observations.m create mode 100644 @ensembleFilter/parseWhich.m create mode 100644 @ensembleFilter/prior.m create mode 100644 @ensembleFilter/rename.m create mode 100644 @ensembleFilter/uncertainties.m create mode 100644 @evolvingEnsemble/evolvingEnsemble.m create mode 100644 @optimalSensor/checkR.m create mode 100644 @optimalSensor/computeMetric.m create mode 100644 @optimalSensor/estimates.m create mode 100644 @optimalSensor/meanMetric.m create mode 100644 @optimalSensor/metric.m create mode 100644 @optimalSensor/optimalSensor.m create mode 100644 @optimalSensor/prior.m create mode 100644 @optimalSensor/psms.m create mode 100644 @optimalSensor/rename.m create mode 100644 @optimalSensor/run.m create mode 100644 @optimalSensor/saveMeanArgs.m create mode 100644 @particleFilter/bayesWeights.m create mode 100644 @particleFilter/bestWeights.m create mode 100644 @particleFilter/estimates.m create mode 100644 @particleFilter/particleFilter.m create mode 100644 @particleFilter/prior.m create mode 100644 @particleFilter/run.m create mode 100644 @particleFilter/weighting.m create mode 100644 @particleFilter/weights.m create mode 100644 README.md create mode 100644 bayfoxPSM.m create mode 100644 prysmCellulose.m create mode 100644 prysmCoral.m create mode 100644 prysmIcecore.m create mode 100644 prysmSpeleothem.m diff --git a/@PSM/PSM.m b/@PSM/PSM.m index 942f7086..5c57c80a 100644 --- a/@PSM/PSM.m +++ b/@PSM/PSM.m @@ -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; @@ -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 @@ -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) diff --git a/@PSM/computeEstimates.m b/@PSM/computeEstimates.m new file mode 100644 index 00000000..01688e3a --- /dev/null +++ b/@PSM/computeEstimates.m @@ -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 \ No newline at end of file diff --git a/@PSM/download.m b/@PSM/download.m index 45281fa0..ed130807 100644 --- a/@PSM/download.m +++ b/@PSM/download.m @@ -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. diff --git a/@PSM/estimate.m b/@PSM/estimate.m index 58de1c67..99bf3d9f 100644 --- a/@PSM/estimate.m +++ b/@PSM/estimate.m @@ -1,124 +1,45 @@ -function[Ye, R] = estimate(X, F) -%% Estimates proxies values from a state vector ensemble +function[Ye, R] = estimate(X, F, throwError) +%% Estimates proxy values from a state vector ensemble given a set of PSMs % -% Ye = PSM.estimate(X, psms) -% Estimates proxies values for a state vector ensemble given a -% set of PSMs. +% Ye = PSM.estimate(X, F) +% Estimate observations from the ensemble using the PSMs % -% [Ye, R] = PSM.estimates(X, psms) -% Also estimates proxy uncertainties (R). +% [Ye, R] = PSM.estimate(X, F) +% Also estimate error-variances from PSMs when possible +% +% [...] = PSM.estimate(X, F, throwError) +% Throws an error when a PSM fails, rather than issuing a warning. Useful +% when trying to debug new PSMs. % % ----- Inputs ----- % -% X: A state vector ensemble. A numeric array with the -% following dimensions (State vector x ensemble members x priors) +% X: A state vector ensemble. A numeric array with up to three dimensions. +% Size is (nState x nEns x nPrior). +% +% F: A cell vector of PSM objects. % -% psms: A set of PSMs. Either a scalar PSM object or a cell -% vector whose elements are PSM objects. +% 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. A numeric array with the dimensions -% (Proxy sites x ensemble members x priors) +% Ye: Observation estimates for the state vector ensemble and PSMs. % -% R: Proxy uncertainty estimates. A numeric array with -% dimensions (proxy sites x ensemble members x priors) +% R: Estimate error-variances as estimated by the PSMs. -% Parse and error check the ensemble. Get size -if isa(X, 'ensemble') - isens = true; - assert(isscalar(X), 'ens must be a scalar ensemble object'); - [nState, nEns] = X.metadata.sizes; - nPriors = 1; - m = matfile(X.file); -else - assert(isnumeric(X), 'X must be numeric'); - assert(ndims(X)<=3, 'X cannot have more than 3 dimensions'); - isens = false; - [nState, nEns, nPriors] = size(X); +% Defaults +if ~exist('throwError','var') || isempty(throwError) + throwError = []; end -% Error check the ensemble, get sizes -dash.assertRealDefined(X, 'X'); - -% Parse the PSM vector -nSite = numel(F); -[F, wasCell] = dash.parseInputCell(F, nSite, 'F'); -name = "F"; +% Error check +F = PSM.setupEstimate(X, F); -% Error check the individual PSMs -for s = 1:nSite - if wasCell - name = sprintf('Element %.f of F', s); - end - dash.assertScalarType(F{s}, name, 'PSM', 'PSM'); - - % Check the rows of the PSM do not exceed the number of rows - if max(F{s}.rows) > nState - error('The ensemble has %.f rows, but %s uses rows that are larger (%.f)', ... - nState, F{s}.messageName(s), max(F{s}.rows)); - end -end - -% Preallocate -Ye = NaN(nSite, nEns, nPriors); -if nargout>1 - R = NaN(nSite, nEns, nPriors); +% Generate the estimates +if nargout==1 + Ye = PSM.computeEstimates(X, F, throwError); +else + [Ye, R] = PSM.computeEstimates(X, F, throwError); end -% Get the values needed to run each PSM -for s = 1:nSite - if ~isens - Xpsm = X(F{s}.rows,:,:); - - % If using an ensemble object, first attempt to read all rows at once - else - try - rows = dash.equallySpacedIndices(F{s}.rows); - Xpsm = m.X(rows,:); - [~, keep] = ismember(F{s}.rows, rows); - Xpsm = Xpsm(keep,:); - - % If unsuccessful, load values iteratively - catch - nRows = numel(F{s}.rows); - Xpsm = NaN(nRows, nEns); - for r = 1:nRows - Xpsm(r,:) = m.X(F{s}.rows(r),:); - end - end - end - - % Get the values for each prior and run the PSM - for p = 1:nPriors - Xrun = Xpsm(:,:,p); - if nargout>1 - [Yrun, Rrun] = F{s}.runPSM(Xrun); - else - Yrun = F{s}.runPSM(Xrun); - end - - % Error check the R output - if nargout>1 - name = sprintf('R values for %s for prior %.f', F{s}.messageName(s), p); - dash.assertVectorTypeN(Rrun, 'numeric', nEns, name); - if ~isrow(Rrun) - error('%s must be a row vector', name); - end - end - - % Error check the Y output - name = sprintf('Y values for %s for prior %.f', F{s}.messageName(s), p); - dash.assertVectorTypeN(Yrun, 'numeric', nEns, name); - if ~isrow(Yrun) - error('%s must be a row vector', name); - end - - % Save - Ye(s,:,p) = Yrun; - if nargout>1 - R(s,:,p) = Rrun; - end - end -end end \ No newline at end of file diff --git a/@PSM/githubLocation.m b/@PSM/githubLocation.m index af897f10..357f822f 100644 --- a/@PSM/githubLocation.m +++ b/@PSM/githubLocation.m @@ -13,12 +13,15 @@ % % commit: The commit hash for the PSM code. A string. -% The Github location associated with each PSM +% The Github location associated with each PSM. Order is 1. key (name for +% DASH), 2. Github Repository, 3. Commit info = [ "bayspar", "https://github.com/jesstierney/BAYSPAR", "13446fb098445d9a8899d9f9a4ea9d81cd916ac9"; "bayspline", "https://github.com/jesstierney/BAYSPLINE", "1e6f9673bcc55b483422c6d6e1b1f63148636094"; "baymag", "https://github.com/jesstierney/BAYMAG", "358de1545d47cbde328fa543c66ab50a20680b00"; - "vslite", "https://github.com/suztolwinskiward/VSLite", "f86cc33ee0eb9a2b9818994542d3c4179e618631" + "bayfox", "https://github.com/jesstierney/bayfoxm", "cb98a5259c7340c3b19669c45a599799b9a491c2"; + "vslite", "https://github.com/suztolwinskiward/VSLite", "f86cc33ee0eb9a2b9818994542d3c4179e618631"; + "prysm", "https://github.com/sylvia-dee/PRYSM", "13dc4fbc1a4493e86a4568d2d83d8495f6f40fe1" ]; % Error check the psm name diff --git a/@PSM/setupEstimate.m b/@PSM/setupEstimate.m new file mode 100644 index 00000000..edecc5ac --- /dev/null +++ b/@PSM/setupEstimate.m @@ -0,0 +1,52 @@ +function[F] = setupEstimate(X, F) +%% Error checks an ensemble and PSMs before estimating proxies. +% +% PSM.checkEstimate(X, F) +% +% ----- Inputs ----- +% +% X: A model ensemble +% +% F: Either a scalar PSM object or a cell vector of PSM objects +% +% ----- Outputs ----- +% +% F: PSM inputs as a cell vector. + +% Error check the ensemble and get sizes +assert(isnumeric(X), 'X must be numeric'); +dash.assertRealDefined(X, 'X'); +assert(ndims(X)<=3, 'X cannot have more than 3 dimensions'); +[nState, nEns, nPrior] = size(X); + +% Parse the PSM vector +nSite = numel(F); +[F, wasCell] = dash.parseInputCell(F, nSite, 'F'); +name = "F"; + +% Error check the individual PSMs +for s = 1:nSite + psm = F{s}; + if wasCell + name = sprintf('Element %.f of F', s); + end + dash.assertScalarType(psm, name, 'PSM', 'PSM'); + + % Check the PSM rows match ensemble sizes + name = psm.messageName(s); + maxrow = max(psm.rows, [], 'all'); + [~, nEns2, nPrior2] = size(psm.rows); + + if maxrow > nState + error('The ensemble has %.f rows, but %s uses rows that are larger (%.f)', ... + nState, psm.messageName(s), maxrow); + elseif nEns2~=1 && nEns2~=nEns + error('The ensemble has %.f members (columns), but %s specifies rows for %.f ensemble members', ... + nEns, psm.messageName(s), nEns2); + elseif nPrior2~=1 && nPrior2~=nPrior + error('The ensemble has %.f priors (elements along dimension 3), but %s specifies rows for %.f priors', ... + nPrior, psm.messageName(s), nPrior2); + end +end + +end \ No newline at end of file diff --git a/@dash/assertScalarType.m b/@dash/assertScalarType.m index 393cfbfc..ba0f5c89 100644 --- a/@dash/assertScalarType.m +++ b/@dash/assertScalarType.m @@ -16,7 +16,7 @@ % string. if ~isscalar(input) || ~isa(input, type) - error('%s must be a scalar %s.', name, typeName); + error('%s must be a %s scalar.', name, typeName); end end \ No newline at end of file diff --git a/@dash/dash.m b/@dash/dash.m index 29279268..ab6c6dcc 100644 --- a/@dash/dash.m +++ b/@dash/dash.m @@ -20,6 +20,7 @@ [X, order] = permuteDimensions(X, index, iscomplete, nDims); tf = bothNaN(A, B); s = loadMatfileFields(file, fields, extName); + str = version; % Input parsing varargout = parseInputs(inArgs, flags, defaults, nPrev); diff --git a/@dash/version.m b/@dash/version.m new file mode 100644 index 00000000..ab36a433 --- /dev/null +++ b/@dash/version.m @@ -0,0 +1,12 @@ +function[ver] = version +%% Returns the version of the current DASH installation +% +% ver = dash.version +% +% ----- Outputs --- +% +% ver: A version string. + +ver = "4.0.0-alpha-5.0.0"; + +end \ No newline at end of file diff --git a/@ensemble/ensemble.m b/@ensemble/ensemble.m index 9ccff2e2..82eff93f 100644 --- a/@ensemble/ensemble.m +++ b/@ensemble/ensemble.m @@ -57,9 +57,9 @@ obj.file = dash.checkFileExists(filename); obj = obj.update; - % Members and variables are unspecified. - obj.members = []; - obj.variables = []; + % Use all members and variables by default + obj.members = 1:obj.metadata.nEns; + obj.variables = obj.metadata.variableNames; % Update name (error checking via ensembleMetadata) if exist('name','var') diff --git a/@ensemble/loadSettings.m b/@ensemble/loadSettings.m index 69d1a6ed..e18f79ea 100644 --- a/@ensemble/loadSettings.m +++ b/@ensemble/loadSettings.m @@ -11,16 +11,9 @@ % % v: Variable indices to load. -% Defaults if unspecified +% Get the object settings members = obj.members; -if isempty(members) - members = 1:obj.metadata.nEns; -end - variables = obj.variables; -if isempty(variables) - variables = obj.metadata.variableNames; -end % Check the variables are still consistent with the .ens file filevars = obj.metadata.variableNames; diff --git a/@ensembleFilter/Rcovariance.m b/@ensembleFilter/Rcovariance.m new file mode 100644 index 00000000..6765f082 --- /dev/null +++ b/@ensembleFilter/Rcovariance.m @@ -0,0 +1,45 @@ +function[Rcov] = Rcovariance(obj, t, s) +%% Returns the R covariance used by a filter in a particular time step. +% +% Rcov = obj.Rcovariance(t) +% Return the R covariance matrix for a queried time step. +% +% Rcov = obj.Rcovariance(t, s) +% Return the R covariance for specific sites. +% +% ----- Inputs ----- +% +% t: The index of an assimilated time step. A scalar positive integer. +% +% s: The indices of observation sites. Either a vector of positive +% integers, or a logical vector with one element per site. +% +% ----- Outputs ----- +% +% Rcov: The R covariance matrix + +% Error check +assert(~isempty(obj.R), 'You have not yet specified R uncertainties'); +dash.assertScalarType(t, 't', 'numeric', 'numeric'); +t = dash.checkIndices(t, 't', obj.nTime, 'number of time steps'); + +% Default and error check sites +if ~exist('s','var') || isempty(s) + s = ~isnan(obj.Y(:,t)); +end +s = dash.checkIndices(s, 's', obj.nSite, 'the number of observation sites'); + +% Empty whichR +if isempty(obj.whichR) + obj.whichR = ones(obj.nTime, 1); +end + +% Get the covariance matrix +c = obj.whichR(t); +if ~obj.Rcov + Rcov = diag(obj.R(s, c)); +else + Rcov = obj.R(s, s, c); +end + +end \ No newline at end of file diff --git a/@ensembleFilter/checkInput.m b/@ensembleFilter/checkInput.m new file mode 100644 index 00000000..f10b28fc --- /dev/null +++ b/@ensembleFilter/checkInput.m @@ -0,0 +1,48 @@ +function[nDim1, nDim2, nDim3] = checkInput(X, name, allowNaN, requireMatrix) +%% Error checks an input as numeric, real, defined, and not empty. +% Optionally checks for a matrix and optionally allows NaN. Returns input +% size. +% +% [nDim1, nDim2, nDim3] = ensembleFilter.checkInput(X, name, allowNaN, requireMatrix) +% +% ----- Inputs ----- +% +% X: The input being checked. Either M, D, R, or Y +% +% name: The name of the input being checked. A string +% +% allowNaN: Scalar logical indicating whether to allow NaN values +% +% requireMatrix: Scalar logical indicating whether to require the input to +% be a matrix. +% +% ----- Outputs ----- +% +% nDim1, nDim2, nDim3: The size of dimensions 1, 2 and 3 for the input. + +% Defaults +if ~exist('allowNaN','var') || isempty(allowNaN) + allowNaN = false; +end +if ~exist('requireMatrix','var') || isempty(requireMatrix) + requireMatrix = false; +end + +% Numeric +assert(isnumeric(X), sprintf('%s must be numeric', name)); + +% Optionally check for matrix +if requireMatrix + assert(ismatrix(X), sprintf('%s must be a matrix', name)); +end + +% No Inf or complex, optionally allow NaN +dash.assertRealDefined(X, name, allowNaN); + +% No empty arrays (they break size checks) +assert(~isempty(X), sprintf('%s cannot be empty', name)); + +% Return sizes +[nDim1, nDim2, nDim3] = size(X); + +end \ No newline at end of file diff --git a/@ensembleFilter/checkMissingR.m b/@ensembleFilter/checkMissingR.m new file mode 100644 index 00000000..0eb3bdaf --- /dev/null +++ b/@ensembleFilter/checkMissingR.m @@ -0,0 +1,31 @@ +function[] = checkMissingR(obj) +%% Checks that each observation has an R variance or R covariances in each +% required time step. +% +% obj.checkMissingR + +% If observations are set, get the time steps for each set of R values +if ~isempty(obj.Y) + obj = obj.finalizeWhich; + for c = 1:obj.nR + times = find(obj.whichR==c); + + % R variances + if ~obj.Rcov + hasobs = any(~isnan(obj.Y(:,times)), 2); + missing = hasobs & isnan(obj.R(:,c)); + assert(~any(missing), sprintf('R variance set %.f must have a value for site %.f', c, find(missing,1))); + + % R covariances need both row and column + else + [site1, site2] = find(isnan(obj.R(:,:,c))); + for k = 1:numel(times) + t = times(k); + missing = find( ~isnan(obj.Y(site1,t)) & ~isnan(obj.Y(site2,t)), 1); + assert(isempty(missing), sprintf('R covariance set %.f must have a covariance betwee site %.f and site %.f', c, site1(missing), site2(missing))); + end + end + end +end + +end \ No newline at end of file diff --git a/@ensembleFilter/ensembleFilter.m b/@ensembleFilter/ensembleFilter.m new file mode 100644 index 00000000..668eb375 --- /dev/null +++ b/@ensembleFilter/ensembleFilter.m @@ -0,0 +1,54 @@ +classdef (Abstract) ensembleFilter + %% Implements common utilities for ensemble-based data assimilation + % filters (i.e. Kalman filters and particle filters) + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + + properties (SetAccess = protected) + name; + + % Essential inputs + Y; + Ye; + X; + whichPrior; + R; + whichR; + Rcov; + + % Sizes + nState = 0; + nEns = 0; + nPrior = 0; + nSite = 0; + nTime = 0; + nR = 0; + end + + % Basic inputs + methods + obj = rename(obj, name); + obj = observations(obj, Y); + obj = uncertainties(obj, R, whichR, isCov); + obj = prior(obj, X, whichPrior); + obj = estimates(obj, Ye, whichPrior); + end + + % User query + methods + Rcov = Rcovariance(obj, t, s); + end + + % Utilities + methods + whichArg = parseWhich(obj, whichArg, name, nIndex, indexName, resetTime); + obj = finalize(obj, actionName); + obj = finalizeWhich(obj); + checkMissingR(obj); + end + methods (Static) + [nDim1, nDim2, nDim3] = checkInput(X, name, allowNaN, requireMatrix); + end + +end \ No newline at end of file diff --git a/@ensembleFilter/estimates.m b/@ensembleFilter/estimates.m new file mode 100644 index 00000000..ebbdee57 --- /dev/null +++ b/@ensembleFilter/estimates.m @@ -0,0 +1,66 @@ +function[obj] = estimates(obj, Ye, whichPrior) +%% Specify the model estimates of the observations/proxies for a filter +% +% obj = obj.estimates(Ye) +% Specify estimates for static or transient priors to use in each time +% step. +% +% obj = obj.estimates(Ye, whichPrior) +% Specifies which prior / associated set of estimates to use in each time +% step. +% +% ----- Inputs ----- +% +% Ye: The model estimates of observations/proxies. A numeric array that +% cannot contain NaN or Inf. (nSite x nEns x nPrior) +% +% whichPrior: A vector with one element per time step. Each element +% is the index of the prior to use for that time step. The indices refer +% to the index along the third dimension of M. (nTime x 1) +% +% ----- Outputs ----- +% +% obj: The updated ensembleFilter object + +% Check data type and get sizes +[nSite, nEns, nPrior] = obj.checkInput(Ye, 'Ye'); + +% Default and parse whichPrior +if ~exist('whichPrior','var') || isempty(whichPrior) + whichPrior = []; +end +resetTime = isempty(obj.Y) && isempty(obj.whichR); +whichPrior = obj.parseWhich(whichPrior, 'whichPrior', nPrior, 'priors', resetTime); +nTime = numel(whichPrior); + +% Size checks +if nSite~=obj.nSite + assert(isempty(obj.Y), sprintf('You previously specified observations for %.f sites, but Ye has %.f sites (rows)', obj.nSite, nSite)); + assert(isempty(obj.R), sprintf('You previously specified R uncertainties for %.f sites, but Ye has %.f sites (rows)', obj.nSite, nSite)); +end +if ~isempty(obj.X) + assert( nEns==obj.nEns, sprintf('You previously specified a prior with %.f ensemble members, but Ye has %.f ensemble members (columns)', obj.nEns, nEns)); + assert( nPrior==obj.nPrior, sprintf('You previously specified a transient prior with %.f priors, but Ye has %.f priors (elements along dimenison 3)', obj.nPrior, nPrior)); +end +if nTime~=obj.nTime && ~isempty(whichPrior) + assert(isempty(obj.Y), sprintf('You previously specified observations for %.f time steps, but here you specify estimates for %.f time steps', obj.nTime, nTime)); + assert(isempty(obj.whichR), sprintf('You previously specified R uncertainties for %.f time steps, but here you specify estimates for %.f time steps', obj.nTime, nTime)); +end + +% If there is a whichPrior for the prior, require the same time steps +if ~isempty(whichPrior) && ~isempty(obj.whichPrior) && ~isempty(obj.X) + assert(isequal(whichPrior, obj.whichPrior), 'The time steps for these estimates do not match the time steps for the associated transient priors. (check the whichPrior input)'); +end + +% Save +obj.Ye = Ye; +obj.nSite = nSite; +obj.nEns = nEns; +obj.nPrior = nPrior; + +if ~isempty(whichPrior) + obj.whichPrior = whichPrior; + obj.nTime = nTime; +end + +end \ No newline at end of file diff --git a/@ensembleFilter/finalize.m b/@ensembleFilter/finalize.m new file mode 100644 index 00000000..37f8e737 --- /dev/null +++ b/@ensembleFilter/finalize.m @@ -0,0 +1,21 @@ +function[obj] = finalize(obj, actionName) +%% Checks that a filter has all essential fields and fills empty whichPrior +% +% obj = obj.finalize(actionName) +% +% ----- Inputs ----- +% +% actionName: The action that the filter is being finalized for. A string +% +% ----- Outputs ----- +% +% obj: The updated ensembleFilter object + +% Check for essential inputs +assert(~isempty(obj.Y), sprintf('You must provide observations before you %s', actionName)); +assert(~isempty(obj.Ye), sprintf('You must provide proxy estimates before you %s.', actionName)); + +% Fill empty whichPrior and whichR +obj = obj.finalizeWhich; + +end diff --git a/@ensembleFilter/finalizeWhich.m b/@ensembleFilter/finalizeWhich.m new file mode 100644 index 00000000..c9720e20 --- /dev/null +++ b/@ensembleFilter/finalizeWhich.m @@ -0,0 +1,10 @@ +function[obj] = finalizeWhich(obj) + +if isempty(obj.whichPrior) + obj.whichPrior = ones(obj.nTime, 1); +end +if isempty(obj.whichR) + obj.whichR = ones(obj.nTime, 1); +end + +end \ No newline at end of file diff --git a/@ensembleFilter/observations.m b/@ensembleFilter/observations.m new file mode 100644 index 00000000..1ac3a257 --- /dev/null +++ b/@ensembleFilter/observations.m @@ -0,0 +1,37 @@ +function[obj] = observations(obj, Y) +%% Specify the observations for a filter +% +% obj = obj.observations(Y) +% +% ----- Inputs ----- +% +% Y: The observations/proxy values. A numeric matrix. Each row has the +% observations for one site over all time steps. Use NaN in time steps +% with no observations. (nSite x nTime) +% +% ----- Outputs ----- +% +% obj: The updated ensembleFilter object + +% Error check +[nSite, nTime] = obj.checkInput(Y, 'Y', true, true); + +% Check that sizes don't conflict +if nSite~=obj.nSite + assert(isempty(obj.Ye), sprintf('You previously specified estimates for %.f sites, but Y has %.f sites (rows)', obj.nSite, nSite)); + assert(isempty(obj.R), sprintf('You previously specified observation uncertainties for %.f sites, but Y has %.f sites (rows)', obj.nSite, nSite)); +end +if nTime~=obj.nTime + assert(isempty(obj.whichPrior), sprintf('You previously specified a transient prior for %.f time steps, but Y has %.f time steps (columns)', obj.nTime, nTime)); + assert(isempty(obj.whichR), sprintf('You previously specified R uncertainties for %.f time steps, but Y has %.f time steps (columns)', obj.nTime, nTime)); +end + +% Set values +obj.nSite = nSite; +obj.nTime = nTime; +obj.Y = Y; + +% Check for missing R +obj.checkMissingR; + +end \ No newline at end of file diff --git a/@ensembleFilter/parseWhich.m b/@ensembleFilter/parseWhich.m new file mode 100644 index 00000000..da7bdec9 --- /dev/null +++ b/@ensembleFilter/parseWhich.m @@ -0,0 +1,54 @@ +function[whichArg] = parseWhich(obj, whichArg, name, nIndex, indexName, resetTime) +%% Parses which* inputs for evolving settings +% +% whichArg = obj.parseWhich(whichArg, name, nIndex, indexName) +% +% ----- Inputs ----- +% +% whichArg: The which* input +% +% name: The name of the which* input in the calling function +% +% nIndex: The indices that can be referenced by whichArg +% +% indexName: The name of the quantity being indexed +% +% resetTime: Scalar logical indicating that whichArg is the only input +% controlling the number of time steps +% +% ----- Outputs ----- +% +% whichArg: The which* input adjusted for any default settings + +% Reset time +if resetTime + obj.nTime = 0; +end + +% If unset, nTime is the number of whichArg elements +empty = isempty(whichArg); +if obj.nTime==0 && ~empty + obj.nTime = numel(whichArg); +elseif obj.nTime==0 + obj.nTime = nIndex; +end + +% If empty, nIndex must be 1 or nTime. Leave empty if 1, but set to indices +% if nTime +if empty + if nIndex==obj.nTime && nIndex~=1 + whichArg = (1:obj.nTime)'; + elseif nIndex~=1 + error(['The number of %s (%.f) does not match the number of time steps ',... + '(%.f), so you must use the "%s" input to specify which ',... + '%s to use in each time step.'], indexName, nIndex, obj.nTime, name, indexName); + end + +% Otherwise, error check the indices +else + dash.assertVectorTypeN(whichArg, 'numeric', obj.nTime, name); + dash.checkIndices(whichArg, name, nIndex, sprintf('the number of %s', indexName)); +end +whichArg = whichArg(:); + +end \ No newline at end of file diff --git a/@ensembleFilter/prior.m b/@ensembleFilter/prior.m new file mode 100644 index 00000000..5314bf51 --- /dev/null +++ b/@ensembleFilter/prior.m @@ -0,0 +1,69 @@ +function[obj] = prior(obj, X, whichPrior) +%% Specify the prior(s) for a filter +% +% obj = obj.prior(X) +% Specify a static or evolving offline prior to use in each time step. +% +% obj = obj.prior(X, whichPrior) +% Specifies evolving offline priors and indicates which prior to use in +% each time step. +% +% ----- Inputs ----- +% +% X: An offline prior. A numeric array that cannot contain Inf or complex values. +% +% Static: X is a matrix (nState x nEns) and will be used as the +% prior in each time step. +% +% Transient: X is an array (nState x nEns x nPrior). If you do +% not provide a second input, must have one prior per time step. +% +% whichPrior: A vector with one element per time step. Each element +% is the index of the prior to use for that time step. The indices refer +% to the index along the third dimension of M. (nTime x 1) +% +% ----- Outputs ----- +% +% obj: The updated ensembleFilter object + +% Error check data type and get sizes +[nState, nEns, nPrior] = obj.checkInput(X, 'X', true); + +% Default and parse whichPrior +if ~exist('whichPrior','var') || isempty(whichPrior) + whichPrior = []; + if ~isempty(obj.whichPrior) && nPrior>1 + whichPrior = obj.whichPrior; + end +end +resetTime = isempty(obj.Y) && isempty(obj.whichR); +whichPrior = obj.parseWhich(whichPrior, 'whichPrior', nPrior, 'priors', resetTime); +nTime = numel(whichPrior); + +% Size checks +if ~isempty(obj.Ye) + assert( nEns==obj.nEns, sprintf('You previously specified estimates for %.f ensemble members, but X has %.f ensemble members (columns)', obj.nEns, nEns)); + assert( nPrior==obj.nPrior, sprintf('You previously specified estimates for %.f priors, but X has %.f priors (elements along dimension 3)', obj.nPrior, nPrior)); +end +if nTime~=obj.nTime && ~isempty(whichPrior) + assert(isempty(obj.Y), sprintf('You previously specified observations for %.f time steps, but here you specify transient priors for %.f time steps', obj.nTime, nTime)); + assert(isempty(obj.whichR), sprintf('You previously specified R uncertainties for %.f time steps, but here you specify transient priors for %.f time steps', obj.nTime, nTime)); +end + +% If there is a whichPrior for estimates, require the same time steps +if ~isempty(whichPrior) && ~isempty(obj.whichPrior) && ~isempty(obj.Ye) + assert(isequal(whichPrior, obj.whichPrior), 'The time steps for the transient priors do not match the time steps specified for the associated estimates. (check the whichPrior input)'); +end + +% Save +obj.X = X; +obj.nState = nState; +obj.nEns = nEns; +obj.nPrior = nPrior; + +if ~isempty(whichPrior) + obj.whichPrior = whichPrior; + obj.nTime = nTime; +end + +end \ No newline at end of file diff --git a/@ensembleFilter/rename.m b/@ensembleFilter/rename.m new file mode 100644 index 00000000..bf830d2f --- /dev/null +++ b/@ensembleFilter/rename.m @@ -0,0 +1,16 @@ +function[obj] = rename(obj, name) +%% Renames an ensembleFilter object +% +% obj = obj.rename( name ) +% +% ----- Inputs ----- +% +% name: The new name. A string. +% +% ----- Outputs ----- +% +% obj: The updated ensembleFilter object + +obj.name = dash.assertStrFlag(name, 'name'); + +end \ No newline at end of file diff --git a/@ensembleFilter/uncertainties.m b/@ensembleFilter/uncertainties.m new file mode 100644 index 00000000..cdae84f0 --- /dev/null +++ b/@ensembleFilter/uncertainties.m @@ -0,0 +1,101 @@ +function[obj] = uncertainties(obj, R, whichR, isCov) +%% Specify observation uncertainties / error-covariance for a filter +% +% obj.uncertainties(Rvar) +% Specify R uncertainties, the error-variance of observations. +% +% obj.uncertainties(Rvar, whichR) +% Specify which R variance to use in each assimilated time step. +% +% obj.uncertainties(Rcov, whichR, isCov) +% Specify error-covariances for the observations. +% +% ----- Inputs ----- +% +% Rvar: Error-variances for the observations. A matrix with one row per +% observations. If a column vector, uses the same error-variance in each +% assimilated time step. If the number of columns matches the number of +% time steps, uses the relevant variances for each time step. For any +% other number of columns, you must specify which R variance to use in +% each time step using the second input. All R variances must be +% positive. +% +% Rcov: Error-covariances for the observations. An array with dimensions +% (nSite x nSite x nR). If you do not specify the second input, must +% either have one element along the third dimension, or one element per +% time step. Each error-covariance matrix must be symmetric with positive +% diagonal elements. +% +% whichR: A vector with one element per time step. Each element is the +% index of the set of R variances/covariances to use for that time step. +% The indices refer to the column of Rvar, or the index along the third +% dimension of Rcov, as appropriate. +% +% isCov: A scalar logical that indicates whether the uncertainties are R +% error-variances (false -- Default) or error-covariances (true) +% +% ----- Outputs ----- +% +% obj: The updated filter object + +% Default and error check covariance toggle +if ~exist('isCov','var') || isempty(isCov) + isCov = false; +end +dash.assertScalarType(isCov, 'isCov', 'logical', 'logical'); + +% Get sizes +if ~isCov + [nSite, nR] = obj.checkInput(R, 'Rvar', true, true); +else + [nSite, nSite2, nR] = obj.checkInput(R, 'Rcov', true, false); + assert(nSite==nSite2, sprintf('The number of rows in Rcov (%.f) does not match the number of columns (%.f)', nSite, nSite2)); +end + +% Default and parse whichR +if ~exist('whichR','var') || isempty(whichR) + whichR = []; +end +resetTime = isempty(obj.Y) && isempty(obj.whichPrior); +whichR = obj.parseWhich(whichR, 'whichR', nR, 'R uncertainties', resetTime); +nTime = numel(whichR); + +% Size checks +if nSite~=nSite + assert(isempty(obj.Y), sprintf('You previously specified observations for %.f sites, but R has %.f sites (rows)', obj.nSite, nSite)); + assert(isempty(obj.Ye), sprintf('You previously specified estimates for %.f sites, but R has %.f sites (rows)', obj.nSite, nSite)); +end +if nTime~=obj.nTime && ~isempty(whichR) + assert(isempty(obj.Y), sprintf('You previously specified observations for %.f time steps, but these R values are for %.f time steps', obj.nTime, nTime)); + assert(isempty(obj.whichPrior), sprintf('You previously specified a transient prior for %.f time steps, but these R values are for %.f time steps', obj.nTime, nTime)); +end + +% Check values are variances or covariances +if ~isCov + assert(~any(R(:)<=0), 'R variances must be greater than 0.'); + +else + for c = 1:nR + Rc = R(:,:,c); + Rnan = isnan(Rc); + Rc(Rnan) = 1; + assert( issymmetric(Rnan) && issymmetric(Rc), sprintf('R covariance %.f is not a symmetric matrix', c)); + assert( ~any(diag(Rc)<=0), sprintf('The diagonal elements of R covariance %.f must be greater than 0', c)); + end +end + +% Save values +obj.R = R; +obj.Rcov = isCov; +obj.nSite = nSite; +obj.nR = nR; + +if ~isempty(whichR) + obj.whichR = whichR; + obj.nTime = nTime; +end + +% Check for missing R values +obj.checkMissingR; + +end \ No newline at end of file diff --git a/@ensembleMetadata/buildFromPrimitives.m b/@ensembleMetadata/buildFromPrimitives.m index b4aa2b26..a8a2bcb9 100644 --- a/@ensembleMetadata/buildFromPrimitives.m +++ b/@ensembleMetadata/buildFromPrimitives.m @@ -76,4 +76,4 @@ end -end +end \ No newline at end of file diff --git a/@ensembleMetadata/closestLatLon.m b/@ensembleMetadata/closestLatLon.m index 0fcf3177..98e28c7e 100644 --- a/@ensembleMetadata/closestLatLon.m +++ b/@ensembleMetadata/closestLatLon.m @@ -1,4 +1,4 @@ -function[rows] = closestLatLon(obj, latlon, varName, verbose) +function[rows] = closestLatLon(obj, latlon, varName, exclude, verbose) %% Returns the state vector rows closest to a set of lat-lon coordinates % % rows = obj.closestLatLon( coordinate ) @@ -8,7 +8,11 @@ % Returns the state vector rows that are closest to the given coordinate % for a particular variable. % -% rows = obj.closestLatLon( coordinate, varName, verbose ) +% rows = obj.closestLatLon( coordinate, varName, exclude ) +% Specify rows that should not be selected as the closest +% latitude-longitude point. +% +% rows = obj.closestLatLon( coordinate, varName, exclude, verbose ) % Optionally disable console messages. % % ----- Inputs ----- @@ -19,6 +23,11 @@ % % varName: The name of a variable in the state vector ensemble. A string. % +% exclude: A logical matrix that indicates which state vector elements to +% exclude from consideration. Must have one row per state vector +% element or one row per state vector element for the specified +% variable, as appropriate. +% % verbose: A scalar logical that indicates whether to return console % messages when determining coordinates (true -- default) or not (false) % @@ -43,12 +52,42 @@ % Get the state vector lat-lon coordinates and get the distance to each svLatLon = obj.latlon(varName, verbose); +nState = size(svLatLon, 1); dist = dash.haversine(svLatLon, latlon); -% Find the state vector rows associated with the minimum distances -rows = find(dist==min(dist)); -if ~isempty(varName) - rows = obj.findRows(varName, rows); +% Default and error check exclude +if ~exist('exclude','var') || isempty(exclude) + exclude = false(nState, 1); +end +assert(islogical(exclude) && ismatrix(exclude), 'exclude must be a logical matrix'); +assert(size(exclude,1)==nState, sprintf('exclude must have %.f rows (one per state vector element', nState)); +nCols = size(exclude, 2); + +% Find the state vector rows associated with the minimum distances for each +% set of excluded indices +for m = 1:nCols + currentDist = dist; + currentDist(exclude(:,m)) = NaN; + currentRows = find( currentDist == min(currentDist) ); + nRows = numel(currentRows); + + % Preallocate the rows for each ensemble member + if m==1 + rows = NaN(nRows, nCols); + end + + % If there are more rows than fit in rows, preallocate more NaN + % rows for the other sets of excluded indices + nAdd = max( nRows-size(rows,1), 0 ); + if nAdd > 0 + rows = [rows; NaN(nAdd, nCols)]; %#ok + end + + % Convert variable rows to state vector rows and save + if ~isempty(varName) + currentRows = obj.findRows(varName, currentRows); + end + rows(1:nRows, m) = currentRows; end end \ No newline at end of file diff --git a/@ensembleMetadata/ensembleMetadata.m b/@ensembleMetadata/ensembleMetadata.m index 695f4f8b..683caf7d 100644 --- a/@ensembleMetadata/ensembleMetadata.m +++ b/@ensembleMetadata/ensembleMetadata.m @@ -179,7 +179,7 @@ meta = columns(obj, cols, varNames, dims, alwaysStruct); [latlon] = latlon(obj, varName, verbose); - rows = closestLatLon(obj, latlon, varName, verbose); + rows = closestLatLon(obj, latlon, varName, exclude, verbose); rows = findRows(obj, varName, varRows); [nState, nEns] = sizes(obj, vars); diff --git a/@ensembleMetadata/sizes.m b/@ensembleMetadata/sizes.m index a6b9f1d6..fcda9607 100644 --- a/@ensembleMetadata/sizes.m +++ b/@ensembleMetadata/sizes.m @@ -2,11 +2,11 @@ %% Returns the size of a state vector ensemble, or variables in the state % vector ensemble. % -% [nState, nEns] = obj.size +% [nState, nEns] = obj.sizes % Returns the size of the complete state vector ensemble. % -% [nEls, nEns] = obj.size(varNames) -% [nEls, nEns] = obj.size(v) +% [nEls, nEns] = obj.sizes(varNames) +% [nEls, nEns] = obj.sizes(v) % Returns the size of variables in the state vector ensemble. % % ------ Inputs ----- diff --git a/@evolvingEnsemble/evolvingEnsemble.m b/@evolvingEnsemble/evolvingEnsemble.m new file mode 100644 index 00000000..05fa4959 --- /dev/null +++ b/@evolvingEnsemble/evolvingEnsemble.m @@ -0,0 +1,114 @@ +classdef evolvingEnsemble + %% Implements evolving ensembles that are too large to fit in active memory + + properties + ensembles; + nState; + nEns; + nPrior = 0; + end + + % User methods + methods + function[obj] = evolvingEnsemble(varargin) + %% Concatenate ensemble objects into an evolving ensemble + % + % obj = evolvingEnsemble(ens1, ens2, ..., ensN) + % + % ----- Inputs ----- + % + % ensN: An ensemble object. All ensemble objects must have the + % same number of LOADED ensemble members and state vector + % elements. + % + % ----- Outputs ----- + % + % obj: The evolvingEnsemble object + + % Require at least one ensemble. Get sizes + assert(~isempty(varargin), 'You must provide at least one ensemble object to create an evolvingEnsemble.'); + dash.assertScalarType(varargin{1}, 'The first input', 'ensemble', 'ensemble'); + [obj.nState, obj.nEns] = varargin{1}.loadedMetadata.sizes; + + % Add to the evolving ensemble + obj.add(varargin); + end + + function[obj] = add(obj, varargin) + % Adds additional ensemble objects to an evolving ensemble + % + % obj = obj.add(ens1, ens2, ..., ensN) + % + % ------ Inputs ----- + % + % ensN: An ensemble object. + % + % ----- Outputs ----- + % + % obj: The updated evolving ensemble + + % Check each input is a scalar ensemble object + nPriors = numel(varargin); + for k = 1:nPriors + dash.assertScalarType(varargin{k}, sprintf('Input %.f', k), 'ensemble', 'ensemble'); + end + + % Ensure all ensemble objects are the same size + for k = 1:nPriors + [S, E] = varargin{k}.loadedMetadata.sizes; + assert(S==obj.nState, sprintf(['Each ensemble object must load %.f ',... + 'state vector elements, but input ensemble %.f loads %.f elements.'], obj.nState, k, S)); + assert(E==obj.nEns, sprintf(['Each ensemble must load %.f ensemble members, ',... + 'but input ensemble %.f loads %.f members.'], obj.nEns, k, E)); + end + + % Add to the evolving ensemble + obj.ensembles = [obj.ensembles; varargin{:}]; + obj.nPrior = obj.nPrior + nPriors; + end + + function[obj] = remove(obj, p) + % Remove ensembles from the evolving ensemble + % + % obj = obj.remove(p) + % + % ----- Inputs ----- + % + % p: The index of the priors to remove from the evolving ensemble + % + % ----- Outputs ----- + % + % obj: The updated evolvingEnsemble object + + % Check the indices are allowed + p = dash.checkIndices(p, 'p', obj.nPrior, 'the number of priors in the evolving ensemble'); + p = unique(p); + assert(numel(p)pf.nEns + error(['You previously specified a weighting scheme for the best %.f particles, ',... + 'but the new prior only has %.f particles (ensemble members, columns)'], pf.weightArgs, pf.nEns); +end + +end \ No newline at end of file diff --git a/@particleFilter/particleFilter.m b/@particleFilter/particleFilter.m new file mode 100644 index 00000000..845d111b --- /dev/null +++ b/@particleFilter/particleFilter.m @@ -0,0 +1,69 @@ +classdef particleFilter < ensembleFilter + %% Implements a particle filter + % + % particleFilter Methods: + % particleFilter - Creates a new particleFilter object + % + % observations - Specify the proxy observations and uncertainties + % prior - Specify the prior ensemble + % estimates - Specify the proxy estimates + % + % weighting - Select a weighting scheme for the particle filter + % run - Run the particle filter + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + + properties + weightType; % A switch for the weighting type: 1. Bayes, 2. Best N average + weightArgs; % Arguments used to calculate weights + end + + % Constructor + methods + function[pf] = particleFilter(name) + %% Creates a new particleFilter object + % + % pf = particleFilter + % Create a new particle filter + % + % pf = particleFilter(name) + % Optionally give the particle filter an identifying name + % + % ----- Inputs ----- + % + % name: An optional name for the particle filter. A string. + % + % ----- Outputs ----- + % + % pf: The new particle filter object + + % Name + if ~exist('name','var') || isempty(name) + name = ""; + end + pf = pf.rename(name); + + % Use the bayes weighting scheme by default + pf = pf.weighting('bayes'); + end + end + + % User methods + methods + pf = prior(pf, X, whichPrior); + pf = estimates(pf, Y, whichPrior); + out = run(pf); + pf = weighting(pf, type, N); + end + + % Weights utilities + methods + w = weights(pf, sse); + end + methods (Static) + w = bayesWeights(sse); + w = bestWeights(sse, N); + end + +end \ No newline at end of file diff --git a/@particleFilter/prior.m b/@particleFilter/prior.m new file mode 100644 index 00000000..74088e6d --- /dev/null +++ b/@particleFilter/prior.m @@ -0,0 +1,45 @@ +function[pf] = prior(pf, X, whichPrior) +%% Specify the prior(s) for a particle filter +% +% pf = pf.prior(X) +% Specify a static or evolving offline prior to use in each time step. +% +% pf = pf.prior(X, whichPrior) +% Specifies evolving offline priors and indicates which prior to use in +% each time step. +% +% ----- Inputs ----- +% +% M: A static offline or evolving offline prior. A numeric array that +% cannot contain Inf or complex values. +% +% Static: M is a matrix (nState x nEns) and will be used as the +% prior in each time step. +% +% Evolving: M is an array (nState x nEns x nPrior). If you do +% not provide a second input, must have one prior per time step. +% +% whichPrior: A vector with one element per time step. Each element +% is the index of the prior to use for that time step. The indices refer +% to the index along the third dimension of M. (nTime x 1) +% +% ----- Outputs ----- +% +% pf: The updated particleFilter object + +% Default whichPrior +if ~exist('whichPrior','var') || isempty(whichPrior) + whichPrior = []; +end + +% Record current nEns and apply standard error checking +nEns = pf.nEns; +pf = prior@ensembleFilter(pf, X, whichPrior); + +% Check for size conflicts with the weighting scheme +assert(nEns==pf.nEns || pf.weightType~=1 || pf.weightArgs>pf.nEns, sprintf(... + ['You previously specified a weighting scheme for the best %.f particles, ',... + 'but the new prior only has %.f particles (ensemble members, columns)'], ... + pf.weightArgs, pf.nEns)); + +end \ No newline at end of file diff --git a/@particleFilter/run.m b/@particleFilter/run.m new file mode 100644 index 00000000..22a4e5ed --- /dev/null +++ b/@particleFilter/run.m @@ -0,0 +1,66 @@ +function[out] = run(pf) +%% Runs a particleFilter object + +% Check for essential inputs. Set defaults +pf = pf.finalize('run a particle filter'); +update = ~isempty(pf.X); + +% Preallocate +out = struct(); +sse = NaN(pf.nEns, pf.nTime); +if update + out.A = NaN(pf.nState, pf.nTime); +end + +% If using R variances, get the innovation for each prior +if ~pf.Rcov + pf.Y = permute(pf.Y, [1 3 2]); + for p = 1:pf.nPrior + t = find(pf.whichPrior==p); + innov = pf.Y(:,:,t) - pf.Ye(:,:,p); + + % Vectorize the sum of squared errors + R = pf.R(:,pf.whichR(t)); + R = permute(R, [1 3 2]); + sse(:,t) = squeeze( sum( (1./R).*(innov).^2, 1, 'omitnan') ); + end + +% For covariances, start by finding the unique R covariances +else + sites = ~isnan(pf.Y)'; + Rcovs = [sites, pf.whichR]; + [Rcovs, ~, whichR] = unique(Rcovs, 'rows'); + + % Invert each covariance + nCovs = size(Rcovs, 1); + for c = 1:nCovs + times = find(whichR==c); + s = sites(times(1),:); + Rinv = pf.Rcovariance(times(1), s)^-1; + + % Get the innovations for each time step + for k = 1:numel(times) + t = times(k); + p = pf.whichPrior(t); + innov = pf.Y(s,t) - pf.Ye(s,:,p); + + % Get the SSE for each ensemble member + for m = 1:pf.nEns + sse(m,t) = innov(:,m)' * Rinv * innov(:,m); + end + end + end +end + +% Compute the particle filter weights +out.weights = pf.weights(sse); + +% Optionally update the prior +if update + for p = 1:pf.nPrior + t = find(pf.whichPrior==p); + out.A(:,t) = pf.X(:,:,p) * out.weights(:, t); + end +end + +end \ No newline at end of file diff --git a/@particleFilter/weighting.m b/@particleFilter/weighting.m new file mode 100644 index 00000000..ed4d4f40 --- /dev/null +++ b/@particleFilter/weighting.m @@ -0,0 +1,42 @@ +function[pf] = weighting(pf, type, N) +%% Select the weighting scheme for a particle filter +% +% pf = pf.weighting('bayes') +% The posterior is computed as the weighted mean of all particles. +% Particles weights are calculated using Bayes' theorem. This is the +% default weighting scheme. +% +% pf = pf.weighting('best', N) +% Computes the posterior using the average of the best N particles in each +% time step. +% +% ----- Inputs ----- +% +% N: The number of particles used to compute the posterior +% +% ----- Outputs ----- +% +% pf: The updated particleFilter object + +% Error check and save the weighting type +type = dash.assertStrFlag(type, 'type'); +assert(any(strcmpi(type, ["best","bayes"])), 'The weighting scheme must either be "best" or "bayes"'); + +% Error check and set up "bayes" scheme +if strcmpi(type, "bayes") + assert(nargin==2, 'The "bayes" weighting scheme does not use additional input arguments'); + pf.weightType = 0; + pf.weightArgs = []; + +% Error check "best" scheme number. Require a scalar integer within the +% number of ensemble members +elseif strcmpi(type, "best") + assert(isscalar(N), 'N must be a scalar'); + dash.assertPositiveIntegers(N, 'N'); + assert(N<=pf.nEns, sprintf('N cannot be larger than the number of ensemble members (%.f)', pf.nEns)); + + pf.weightType = 1; + pf.weightArgs = N; +end + +end \ No newline at end of file diff --git a/@particleFilter/weights.m b/@particleFilter/weights.m new file mode 100644 index 00000000..cf6a1f17 --- /dev/null +++ b/@particleFilter/weights.m @@ -0,0 +1,22 @@ +function[w] = weights(pf, sse) +%% Returns weights for a particle filter. +% +% w = pf.weights(sse) +% +% ----- Inputs ----- +% +% sse: Sum of squared errors from a particle filter +% +% ----- Outputs ----- +% +% w: The weights + +% Switch to the appropriate calculator +if pf.weightType==0 + w = pf.bayesWeights(sse); +elseif pf.weightType==1 + N = pf.weightArgs; + w = pf.bestWeights(sse, N); +end + +end \ No newline at end of file diff --git a/@stateVectorVariable/stateVectorVariable.m b/@stateVectorVariable/stateVectorVariable.m index 768c0110..278009d5 100644 --- a/@stateVectorVariable/stateVectorVariable.m +++ b/@stateVectorVariable/stateVectorVariable.m @@ -2,6 +2,9 @@ % This class implements a custom structure to hold design parameters % for a variable in a state vector. + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = private) % Set by constructor name; % The identifying name of the variable diff --git a/README.md b/README.md new file mode 100644 index 00000000..5c7916ba --- /dev/null +++ b/README.md @@ -0,0 +1,28 @@ +# DASH +A package for paleoclimate data assimilation workflow. +[Website](https://jonking93.github.io/DASH) + +DASH is a Matlab package to help implement paleoclimate data assimilation. It includes modules to help +1. Catalogue and organize climate data +2. Design state vector ensembles +3. Implement proxy system modles +4. Kalman Filters +5. Particle Filters, and +6. Optimal Sensors +and automate other workflow tasks. + +You can find comprehensive documentation for DASH at: [https://jonking93.github.io/DASH](https://jonking93.github.io/DASH) + + +### Branches +The following is an overview of the branches of DASH + +main: This branch holds the most recent release of DASH. Any code on the main branch (following the v4.0.0 release) should run succesfully. + +dev: This is the branch for active development of the DASH code. This branch may have the most cutting edge updates to DASH, but may change and/or break without warning. + +gh-pages: This branch holds the source code for the DASH website + +gh-pages-dev: Use this branch for website development before pushing to production. + + diff --git a/bayfoxPSM.m b/bayfoxPSM.m new file mode 100644 index 00000000..b670cad9 --- /dev/null +++ b/bayfoxPSM.m @@ -0,0 +1,123 @@ +classdef bayfoxPSM < PSM + % Implements the BayFOX PSM, a Bayesian model for d18Oc of planktic + % foraminifera by Jess Tierney. + % + % Find it on Github at: https://github.com/jesstierney/bayfoxm + % + % Or read the paper: + % Malevich, S. B., Vetter, L., & Tierney, J. E. (2019). Global Core Top + % Calibration of δ18O in Planktic Foraminifera to Sea Surface + % Temperature. Paleoceanography and Paleoclimatology, 34(8), 1292-1315. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2020 + + properties + species; + end + + % Run directly + methods (Static) + function[Y, R] = run(t, d18Osw, species) + %% Runs the BayFOX planktic foraminifera forward model + % + % [d18Oc, R] = bayfoxPSM.run(SST, d18Osw, species) + % Applies the BayFOX PSM to a set of sea surface temperatues + % and d18O (seawater) values to estimate d18O of calcite. + % + % ----- Inputs ----- + % + % Please see the BayFOX documentation of the function + % "bayfox_forward" for details on inputs. + % + % ----- Outputs ----- + % + % d18Oc: Estimates of d18O of foram calcite. A vector. + % + % R: Error-variance uncertainties estimated from the posterior. A vector. + + assert(isnumeric(t) && isvector(t), 'T must be a numeric vector'); + assert(isnumeric(d18Osw) && isvector(d18Osw), 'd18Osw must be a numeric vector'); + + % Run the PSM, collect Y and R from posterior + d18Oc = bayfox_forward(t, d18Osw, species); + Y = mean(d18Oc, 2); + R = var(d18Oc, [], 2); + + % Shape to the SST vector + if isrow(t) + Y = Y'; + R = R'; + end + end + end + + methods + % Constructor + function[obj] = bayfoxPSM(rows, species, name) + %% Creates a new bayfoxPSM object + % + % obj = bayfoxPSM(rows, species) + % Creates a BayFOX PSM for a set of SST and d18O(sea-water) + % values for a planktic foraminifera species. + % + % obj = bayfoxPSM(rows, species, name) + % Optionally names the PSM. + % + % ----- Inputs ----- + % + % rows: The state vector rows with the SST values and the + % d18Osw values needed to run the PSM. The first row should + % indicate SST values, the second row should indicate + % d18Osw. + % + % species: A string indicating the target species. Options are: + % 'bulloides' = G. bulloides + % 'incompta' = N. incompta + % 'pachy' = N. pachyderma + % 'ruber' = G. ruber + % 'sacculifer' = T. sacculifer + % 'all' = use the pooled annual (non-species specific) model + % 'all_sea' = use the pooled seasonal (non-species specific) model + % + % name: An optional name for the PSM. A string + % + % ----- Outputs ----- + % + % obj: The new baysplinePSM object + + % Set name, estimatesR, rows + if ~exist('name','var') + name = ""; + end + obj@PSM(name, true); + obj = obj.useRows(rows, 2); + + % Set the species parameter + obj.species = species; + end + + function[Y, R] = runPSM(obj, X) + % Runs a BAYFOX PSM given data from a state vector ensemble + % + % [Y, R] = obj.run(X) + % + % ----- Inputs ----- + % + % X: The data used to run the PSM. A numeric matrix with two + % rows. First row is SST, second row is d18Osw + % + % ----- Outputs ----- + % + % Y: Estimates of d18O of calcite + % + % R: Proxy uncertainties estimated from the posterior + + % Split out the climate variables and run the model + T = X(1,:); + d18Osw = X(2,:); + [Y, R] = bayfoxPSM.run(T, d18Osw, obj.species); + end + end + +end \ No newline at end of file diff --git a/baymagPSM.m b/baymagPSM.m index d05742d7..6dae932a 100644 --- a/baymagPSM.m +++ b/baymagPSM.m @@ -1,7 +1,13 @@ classdef baymagPSM < PSM - % Implements the BAYMAG Mg/Ca PSM. - % Requires the Curve Fitting Toolbox + % Implements the BAYMAG PSM, A Bayesian model for Mg/Ca of planktic + % foraminiera by Jess Tierney. % + % Prerequisites: Requires the Curve Fitting Toolbox + % + % Find it on Github at: https://github.com/jesstierney/bayfoxm + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 properties age; @@ -13,71 +19,103 @@ options; end + % Run directly + methods (Static) + function[Y, R] = run(age, t, omega, salinity, pH, clean, species, options) + %% Runs the BayMAG planktic foraminifera forward model + % + % [mgca, R] = baymagPSM.run(age, SST, omega, salinity, pH, clean, species, options) + % Applies the BayMAG PSM to a set of sea surface temperatures + % to estimate Mg/Ca of planktic foraminifera. + % + % ----- Inputs ----- + % + % Please see the BayMAG documentation of the function + % "baymag_forward_ln" for details on the inputs. + % + % options: A cell vector holding the "varargin" inputs for + % "baymax_forward_ln". + % + % ----- Outputs ----- + % + % mgca: Mg/Ca estimates of planktic foraminifera. A vector. + % + % R: Error-variance uncertainties estimated from the + % posterior. A vector. + + % Error check the SSTs + assert(isnumeric(t), 't must be numeric'); + assert(isvector(t), 't must be a vector'); + + % Default for unspecified options. + if ~exist('options','var') || isempty(options) + options = {}; + end + + % Run the PSM + mgca = baymag_forward_ln(age, t, omega, salinity, pH, clean, species, options{:}); + Y = mean(mgca, 2); + R = var(mgca, [], 2); + + % Shape to the sst vector + if isrow(t) + Y = Y'; + R = R'; + end + end + end + methods % Constructor function[obj] = baymagPSM(row, age, omega, salinity, pH, clean, species, options, name) - %% Creates a new BAYMAG PSM + %% Creates a new baymagPSM object % % obj = baymagPSM(row, age, omega, salinity, pH, clean, species) - % Creates a BAYMAG PSM for a site. + % Creates a BAYMAG PSM for a set of SST values for a proxy site. % % obj = baymagPSM(..., options) % Specify optional parameters. (Please see the documentation - % for the baymag_forward function in the BAYMAG package) + % for the "baymag_forward_ln" function in the BAYMAG package) % % obj = baymagPSM(..., options, name) % Optionally name the PSM % % ----- Inputs ----- % - % Please see the documentation of the baymag_forward function - % in the BAYMAG package for more detailed descriptions of the - % PSM parameters. - % - % row: The state vector row with the SST values required to run - % the PSM. A positive integer. - % - % age: The ages of the temperatures. A numeric scalar. - % - % omega: A scalar indicating bottom water saturation. - % - % salinity: A scalar of salinity (psu) - % - % pH: Scalar of pH (total scale) - % - % clean: A flag used to describe the cleaning technique. - % - % species: A string indicating the target species. + % row: The state vector row with the SST values needed to run + % the PSM. % % options: A cell vector with up to 4 elements. Elements are - % the optional parameters detailed in the BAYMAG documentation. + % the optional parameters detailed in the BAYMAG + % documentation of the "baymag_forward_ln" function. % % name: An optional name for the PSM. A string. % + % Please see the documentation of the "baymag_forward_ln" + % function in the BAYMAG package for details on the remaining + % inputs. + % % ----- Outputs ----- % % obj: The new baysparPSM object - % Set the name and estimatesR + % Set the name, rows, and estimatesR if ~exist('name','var') name = ""; end obj@PSM(name, true); - - % Check and set row - assert(isscalar(row), 'row must be a scalar'); - obj = obj.useRows(row); + obj = obj.useRows(row, 1); % Error check the optional argument cell if ~exist('options','var') || isempty(options) options = {}; + else + dash.assertVectorTypeN(options, 'cell', [], 'options'); + assert(numel(options)<=4, 'options cannot have more than 4 elements.'); end - dash.assertVectorTypeN(options, 'cell', [], 'options'); - assert(numel(options)<=4, 'options cannot have more than 4 elements.'); % Set the inputs obj.age = age; - obj.t = t; obj.omega = omega; obj.salinity = salinity; obj.pH = pH; @@ -88,13 +126,14 @@ % Run the PSM function[Y, R] = runPSM(obj, T) - %% Runs a BAYMAG PSM + %% Runs a BAYMAG PSM given data from a state vector ensemble % % Y = obj.run(T) % % ----- Inputs ----- % - % T: The temperatures use to run the PSM. A numeric row vector + % T: The sea surface temperatures use to run the PSM. A + % numeric row vector. % % ----- Outputs ----- % @@ -102,37 +141,9 @@ % % R: Proxy uncertainties estimated from the posterior - % Run the forward model, convert to row - mgca = baymag_forward_ln(obj.age, T, obj.omega, obj.salinity, obj.pH, obj.clean, obj.species, obj.options{:}); - Y = mean(mgca, 2)'; - R = var(mgca, [], 2)'; + % Run the forward model + [Y, R] = baymagPSM.run(obj.age, T, obj.omega, obj.salinity, obj.pH, obj.clean, obj.species, obj.options{:}); end end - - % Run directly - methods (Static) - function[Y, R] = run(age, t, omega, salinity, pH, clean, species, options) - - % Error check the SSTs - assert(isnumeric(t), 't must be numeric'); - assert(isvector(t), 't must be a vector'); - - % Default for unspecified options. - if ~exist('options','var') || isempty(options) - options = {}; - end - - % Run the PSM - mgca = baymag_forward_ln(age, t, omega, salinity, pH, clean, species, options{:}); - Y = mean(mgca, 2); - R = var(mgca, [], 2); - - % Shape to the sst vector - if isrow(ssts) - Y = Y'; - R = R'; - end - end - end end \ No newline at end of file diff --git a/baysparPSM.m b/baysparPSM.m index f372d39f..dfc84531 100644 --- a/baysparPSM.m +++ b/baysparPSM.m @@ -1,7 +1,17 @@ classdef baysparPSM < PSM - % Implements the BAYSPAR TEX86 PSM. - % Requires the Curve Fitting Toolbox + % Implements the BAYSPAR PSM, a Bayesian model for TEX86 by Jess Tierney % + % Prerequisites: Requires the Curve Fitting Toolbox + % + % Find it on Github at: https://github.com/jesstierney/BAYSPAR + % + % Or read the paper: + % Tierney, J.E. & Tingley, M.P. (2014) A Bayesian, spatially-varying + % calibration model for the TEX86 proxy. Geochimica et Cosmochimica + % Acta, 127, 83-106. https://doi.org/10.1016/j.gca.2013.11.026. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 properties lat; @@ -9,6 +19,52 @@ options; end + methods (Static) + % Run directly + function[Y, R] = run(lat, lon, ssts, options) + %% Runs the BaySPAR TEX86 forward model. + % + % [tex86, R] = baysparPSM.run(lat, lon, ssts, options) + % Applies the BaySPAR PSM to a set of sea surface temperatures + % for a proxy site. + % + % ----- Inputs ----- + % + % Please see the BayFOX documentation on the function + % "TEX_forward" for details on the inputs. + % + % options: A cell vector holding the "varargin" inputs for + % "TEX_forward". + % + % ----- Outputs ----- + % + % tex86: TEX86 estimates. A vector + % + % R: Error-variance uncertainties estimated from the posterior. + % A vector. + + % Error check the SSTs + assert(isnumeric(ssts), 'ssts must be numeric'); + assert(isvector(ssts), 'ssts must be a vector'); + + % Default options. Error check PSM parameters + if ~exist('options','var') || isempty(options) + options = {}; + end + + % Run the PSM + tex = TEX_forward(lat, lon, ssts, options{:}); + Y = mean(tex, 2); + R = var(tex, [], 2); + + % Shape to the sst vector + if isrow(ssts) + Y = Y'; + R = R'; + end + end + end + methods % Constructor function[obj] = baysparPSM(row, lat, lon, options, name) @@ -46,21 +102,25 @@ % % obj: The new baysparPSM object - % Set the name and estimatesR + % Set the name, rows, and estimatesR if ~exist('name','var') name = ""; end obj@PSM(name, true); - - % Check and set row - assert(isscalar(row), 'row must be a scalar'); - obj = obj.useRows(row); + obj = obj.useRows(row,1); % Error check the PSM inputs. if ~exist('options','var') || isempty(options) options = {}; + else + dash.assertVectorTypeN(options, 'cell', [], 'options'); + assert(numel(options)<=3, 'options cannot have more than 3 elements'); end - [obj.lat, obj.lon, obj.options] = baysparPSM.checkInputs(lat, lon, options); + + % Save the inputs + obj.lat = lat; + obj.lon = lon; + obj.options = options; end % Run the PSM @@ -71,7 +131,7 @@ % % ----- Inputs ----- % - % SSTs: The SSTs use to run the PSM. A numeric row vector + % SSTs: The SSTs use to run the PSM. A numeric row vector. % % ----- Outputs ----- % @@ -80,72 +140,8 @@ % R: Proxy uncertainties estimated from the posterior % Run the forward model, convert to row - tex = TEX_forward(obj.lat, obj.lon, SSTs, obj.options{:}); - Y = mean(tex,2)'; - R = var(tex,[],2)'; - - end - end - - methods (Static) - % Run directly - function[Y, R] = run(lat, lon, ssts, options) + [Y, R] = baysparPSM.run(obj.lat, obj.lon, SSTs, obj.options); - % Error check the SSTs - assert(isnumeric(ssts), 'ssts must be numeric'); - assert(isvector(ssts), 'ssts must be a vector'); - - % Default options. Error check PSM parameters - if ~exist('options','var') || isempty(options) - options = {}; - end - baysparPSM.checkInputs(lat, lon, options); - - % Run the PSM - tex = TEX_forward(lat, lon, ssts, options{:}); - Y = mean(tex, 2); - R = var(tex, [], 2); - - % Shape to the sst vector - if isrow(ssts) - Y = Y'; - R = R'; - end end - - % Error check the PSM inputs - function[lat, lon, options] = checkInputs(lat, lon, options) - %% Error checks the PSM parameters - % - % [lat, lon, options] = baysparPSM.checkInputs(lat, lon, options) - % - % ----- Inputs ----- - % - % lat: Site latitude. A numeric scalar. - % - % lon: Site longitude. A numeric scalar. - % - % options: Calibration options. A cell vector with up to 3 - % elements. - - % Error check lat and lon - dash.assertScalarType(lat, 'lat', 'numeric', 'numeric'); - dash.assertScalarType(lon, 'lon', 'numeric', 'numeric'); - dash.assertRealDefined(lat, 'lat'); - dash.assertRealDefined(lon, 'lon'); - - % Error check calibration options - dash.assertVectorTypeN(options, 'cell', [], 'options'); - if numel(options)>3 - error('options cannot contain more than 3 elements.'); - end - end end - -end - - - - - - \ No newline at end of file +end \ No newline at end of file diff --git a/baysplinePSM.m b/baysplinePSM.m index d3180dc2..9e7d6b96 100644 --- a/baysplinePSM.m +++ b/baysplinePSM.m @@ -1,15 +1,59 @@ classdef baysplinePSM < PSM - % Implements the BAYSPLINE UK37 PSM - % Requires the Curve Fitting Toolbox + % Implements the BAYSPLINE PSM, a Bayesian model for UK37 by Jess + % Tierney. % - % baysplinePSM Methods: - % baysplinePSM - Creates a new baysplinePSM object - % run - Runs BAYSPLINE directly + % Prerequisites: Requires the Curve Fitting Toolbox + % + % Find it on Github at: https://github.com/jesstierney/BAYSPLINE + % + % Or read the paper: + % Tierney, J.E. & Tingley, M.P. (2018) BAYSPLINE: A New Calibration for + % the Alkenone Paleothermometer. Paleoceanography and + % Paleoclimatology 33, 281-301, [http://doi.org/10.1002/2017PA003201]. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + + % Run directly + methods (Static) + function[Y, R] = run(ssts) + %% Runs the BAYSPLINE UK37 forward model. + % + % [UK37, R] = baysplinePSM.run(ssts) + % Applies BAYSPLINE to a set of SSTs to estimate UK37 values + % and associated uncertainty. + % + % ----- Inputs ----- + % + % ssts: A vector of sea surface temperatures (in Celsius). + % + % ----- Outputs ----- + % + % Y: A vector of UK37 estimates. + % + % R: UK37 uncertainties estimated from the posterior. A vector. + + % Error check ssts + assert(isnumeric(ssts), 'ssts must be numeric'); + assert(isvector(ssts), 'ssts must be a vector'); + + % Run the PSM. Get the posterior mean and estimate R + uk = UK_forward(ssts); + Y = mean(uk, 2); + R = var(uk, [], 2); + + % Shape to the sst vector + if isrow(ssts) + Y = Y'; + R = R'; + end + end + end methods % Constructor function[obj] = baysplinePSM(row, name) - %% Creates a new BAYSPLINE PSM + %% Creates a new baysplinePSM object % % obj = baysplinePSM(row) % Creates a BAYSPLINE PSM for a set of SST values. @@ -20,7 +64,7 @@ % ----- Inputs ----- % % row: The state vector row with the SST values required to run - % the PSM. A positive integer. + % the PSM. % % name: An optional name for the PSM. A string % @@ -28,15 +72,12 @@ % % obj: The new baysplinePSM object - % Set name, estimatesR + % Set name, estimatesR, row if ~exist('name','var') name = ""; end obj@PSM(name, true); - - % Check and set rows - assert(isscalar(row), 'row must be a scalar'); - obj = obj.useRows(row); + obj = obj.useRows(row, 1); end @@ -44,60 +85,22 @@ function[UK, R] = runPSM(~, SSTs) %% Runs a baysplinePSM object % - % [UK, R] = obj.run(SSTs) + % [UK37, R] = obj.run(SSTs) % % ----- Inputs ----- % - % SSTs: The SSTs for the BAYSPLINE model. A numeric row vector. + % SSTs: The SSTs used to run the PSM. A numeric row vector. % % ----- Outputs ----- % - % UK: A row vector of UK37 estimates. + % UK37: A row vector of UK37 estimates. A vector. % - % R: Uncertainty estimates from the posterior. - + % R: Error-variance uncertainty estimated from the posterior. + % A vector. + % Currently, there are no additional parameters, so can just % call the static method directly. [UK, R] = baysplinePSM.run(SSTs); end end - - % Run directly - methods (Static) - function[Y, R] = run(ssts) - %% Runs the BAYSPLINE PSM directly - % - % [Y, R] = baysplinePSM.run(ssts) - % Applies BAYSPLINE to a set of SSTs to estimate UK37 values - % and associated uncertainty. - % - % ----- Inputs ----- - % - % ssts: A vector of sea surface temperatures (in Celsius). - % Please see the BAYSPLINE documentation for more details. - % - % ----- Outputs ----- - % - % Y: A vector of UK37 estimates. - % - % R: UK37 uncertainties estimated from the posterior. A vector. - - % Error check ssts - assert(isnumeric(ssts), 'ssts must be numeric'); - assert(isvector(ssts), 'ssts must be a vector'); - - % Run the PSM. Get the posterior mean and estimate R - uk = UK_forward(ssts); - Y = mean(uk, 2); - R = var(uk, [], 2); - - % Shape to the sst vector - if isrow(ssts) - Y = Y'; - R = R'; - end - end - end - -end - \ No newline at end of file +end \ No newline at end of file diff --git a/dataSource.m b/dataSource.m index c7c8a9b4..aad0ca5f 100644 --- a/dataSource.m +++ b/dataSource.m @@ -4,6 +4,9 @@ % implement functionality for different types of data sourcess. (For % example, netCDF and .mat files and opendap files). + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = protected) source; % The data source. A filename or opendap url var; % (For hdf data sources) The name of the variable. diff --git a/linearPSM.m b/linearPSM.m index 616ef1bd..37fde23f 100644 --- a/linearPSM.m +++ b/linearPSM.m @@ -1,11 +1,7 @@ classdef linearPSM < PSM - % Implements a linear PSM of the following equation + % Implements a general linear PSM of the following form: % % Y = a1 X1 + a2 X2 + ... + an Xn + b - % - % linearPSM Methods: - % linearPSM - Creates a new linearPSM object - % run - Runs a linear PSM directly % ----- Written By ----- % Jonathan King, University of Arizona, 2019-2020 @@ -30,7 +26,7 @@ % obj = linearPSM(rows, slope, ...) % Uses the same slope for all variables. % - % obj = linaerPSM(rows, slopes) + % obj = linearPSM(rows, slopes) % Uses a default intercept of 0. % % ----- Inputs ----- @@ -62,7 +58,7 @@ if ~exist('intercept','var') intercept = []; end - [obj.slopes, obj.intercept] = linearPSM.checkInputs(slopes, intercept, numel(obj.rows)); + [obj.slopes, obj.intercept] = linearPSM.checkInputs(slopes, intercept, size(obj.rows,1)); end % Run the PSM @@ -104,7 +100,8 @@ % Each row is a different X value. % % slopes: The slopes for each variable. A numeric vector in the - % order as the listed rows. + % order as the listed rows. If a scalar, uses the same slope + % for all variables. % % intercept: The intercept for the linear equation. A numeric % scalar. @@ -160,6 +157,7 @@ elseif nSlopes~=nRows error('The number of slopes (%.f) does not match the number of rows (%.f)', nSlopes, nRows); end + slopes = slopes(:); % Default and error check intercept if ~exist('intercept','var') || isempty(intercept) diff --git a/matSource.m b/matSource.m index 9f23db20..8a3d3c8d 100644 --- a/matSource.m +++ b/matSource.m @@ -1,6 +1,9 @@ classdef matSource < dataSource %% Reads data from a .mat file data source + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = private) m; % A matfile object end diff --git a/ncSource.m b/ncSource.m index e987694f..f5322c86 100644 --- a/ncSource.m +++ b/ncSource.m @@ -2,6 +2,9 @@ %% Used to read data from source based on a NetCDF format. Includes % local NetCDF files and OPeNDAP requests. + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = protected) nDims; % The number of dimensions recorded in the NetCDF end diff --git a/opendapSource.m b/opendapSource.m index 1426910b..90217d73 100644 --- a/opendapSource.m +++ b/opendapSource.m @@ -2,6 +2,9 @@ %% Reads data accessed via an OPeNDAP url. Attempts to load and save % the entire variable when using repeated loads to increase speed. + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = private) X; % The loaded and saved dataset attemptFullLoad; % Whether to attempt to load the entire dataset diff --git a/posteriorCalculation.m b/posteriorCalculation.m index 62cf4e03..9ebc474d 100644 --- a/posteriorCalculation.m +++ b/posteriorCalculation.m @@ -2,6 +2,9 @@ %% Implements calculations that require the posterior deviations % (and optionally the posterior mean) from a Kalman Filter + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (Abstract, SetAccess = immutable) timeDim; % The time dimension in the output quantity end diff --git a/posteriorIndex.m b/posteriorIndex.m index ea1b3821..d045657c 100644 --- a/posteriorIndex.m +++ b/posteriorIndex.m @@ -1,5 +1,8 @@ classdef posteriorIndex < posteriorCalculation + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = immutable) name; timeDim = 2; @@ -7,26 +10,39 @@ properties (SetAccess = private) weights; - denom; rows; + nanflag; end methods function[siz] = outputSize(~, ~, nTime, nEns) siz = [nEns, nTime]; end - function[obj] = posteriorIndex(name, weights, rows) + function[obj] = posteriorIndex(name, weights, rows, nanflag) obj.name = name; obj.weights = weights; - obj.denom = sum(weights,1); obj.rows = rows; + obj.nanflag = nanflag; end function[index] = calculate(obj, Adev, Amean) + + % Extract the relevant data Amean = Amean(obj.rows,:); Adev = Adev(obj.rows,:); - devsum = sum(obj.weights .* Adev, 1) ./ obj.denom; - meansum = sum(obj.weights .* Amean, 1) ./ obj.denom; + % Exclude NaN values + nans = isnan(obj.weights); + Amean(nans,:) = NaN; + Adev(nans, :) = NaN; + + % Get the denominator after accounting for NaN + w = obj.weights; + w(isnan(Amean(:,1))) = NaN; + denom = sum(w, 'omitnan'); + + % Get the index + devsum = sum(obj.weights.*Adev, 1, obj.nanflag) ./ denom; + meansum = sum(obj.weights.*Amean, 1, obj.nanflag) ./ denom; index = devsum' + meansum; end function[name] = outputName(obj) diff --git a/posteriorPercentiles.m b/posteriorPercentiles.m index fa703c4c..0fb806cd 100644 --- a/posteriorPercentiles.m +++ b/posteriorPercentiles.m @@ -1,6 +1,9 @@ classdef posteriorPercentiles < posteriorCalculation %% Calculates percentiles of a posterior ensemble. + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = immutable) timeDim = 3; end diff --git a/posteriorVariance.m b/posteriorVariance.m index 9f19a335..1977785e 100644 --- a/posteriorVariance.m +++ b/posteriorVariance.m @@ -1,6 +1,9 @@ classdef posteriorVariance < posteriorCalculation %% Calculates the variance of a posterior ensemble. + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties (SetAccess = immutable) timeDim = 2; end diff --git a/progressbar.m b/progressbar.m index ac3868b1..0ddc8d83 100644 --- a/progressbar.m +++ b/progressbar.m @@ -1,6 +1,9 @@ classdef progressbar < handle % Implements a progress bar + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 + properties show; handle; diff --git a/prysmCellulose.m b/prysmCellulose.m new file mode 100644 index 00000000..6d7e9ddc --- /dev/null +++ b/prysmCellulose.m @@ -0,0 +1,139 @@ +classdef prysmCellulose < PSM + %% Implements the cellulose module of the PRYSM package + % + % Prerequisites: Python 3.4, numpy, scipy, and rpy2. See the PRYSM + % documentation for more details. + % + % Find PRYSM on Github at: https://github.com/sylvia-dee/PRYSM + % + % Or read the paper: + % Dee, S. G., Russell, J. M., Morrill, C., Chen, Z., & Neary, A. + % (2018). PRYSM v2. 0: A proxy system model for lacustrine archives. + % Paleoceanography and Paleoclimatology, 33(11), 1250-1269. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2020 + + properties + d18Os; + d18Op; + d18Ov; + model = "Evans"; + useIsotopes = true; + end + + methods (Static) + function[d18O] = run(T, P, RH, d18Os, d18Op, d18Ov, model, useIsotopes) + %% Runs the PRYSM Cellulose module + % + % d18O = prysmCellulose.run(T, P, RH, d18Os, d18Op, d18Ov, model, useIsotopes) + % Runs the Prysm cellulose model given temperature, + % precipitation, and relative humidity for a site. + % + % ----- Inputs ----- + % + % Please see the documentation of PRYSM module + % "psm.cellulose.sensor.cellulose_sensor" for details on the + % inputs. + % + % ----- Outputs ----- + % + % d18O: Estimates of d18O of cellulose. + + % Placeholder for time input + time = 1:size(T,2); + + % Convert numpy arrays + time = py.numpy.array(time); + T = py.numpy.array(T); + P = py.numpy.array(P); + RH = py.numpy.array(RH); + + % Get the model switch + if strcmp(model, 'Roden') + flag = 0; + elseif strcmp(model, 'Evans') + flag = 1; + end + + % Run the model, convert to Matlab numeric + d18O = py.psm.cellulose.sensor.cellulose_sensor(time, T, P, RH, d18Os, d18Op, d18Ov, flag, useIsotopes); + d18O = numeric(d18O); + end + end + + methods + function[obj] = prysmCellulose(rows, d18Os, d18Op, d18Ov, model, useIsotopes, name) + %% Creates a new prysmCellulose object + % + % obj = prysmCellulose(rows, d18Os, d18Op, d18Ov, model, useIsotopes) + % Creates a PRYSM cellulose PSM for temperature, precipitation, + % and relative humidity for a proxy site. + % + % obj = prysmCellulose(..., name) + % Optionally name the PSM + % + % ----- Inputs ----- + % + % rows: The state vector rows with the temperature, + % precipitation, and relative humidity data for the proxy site. + % First row should be temperature, second is precipitation, + % third should be relative humidity. + % + % model: A string. Either "Roden" or "Evans" + % + % name: An optional name for the PSM. A string. + % + % Please see the documentation of PRYSM module + % "psm.cellulose.sensor.cellulose_sensor" for details on the + % remaining inputs. + + % Name, R estimation, rows + if ~exist('name','var') + name = ""; + end + obj@PSM(name, false); + obj = obj.useRows(rows, 3); + + % Error check model inputs + dash.assertScalarType(d18Os, 'd18Os', 'numeric', 'numeric'); + dash.assertScalarType(d18Op, 'd18Op', 'numeric', 'numeric'); + dash.assertScalarType(d18Ov, 'd18Ov', 'numeric', 'numeric'); + obj.d18Os = d18Os; + obj.d18Op = d18Op; + obj.d18Ov = d18Ov; + + dash.assertStrFlag(model); + if exist('model','var') && ~isempty(model) + dash.checkStrsInList(model, ["Roden","Evans"], 'model', 'model name'); + obj.model = model; + end + if exist('useIsotopes','var') && ~isempty(useIsotopes) + dash.assertScalarType(useIsotopes, 'useIsotopes', 'logical', 'logical'); + obj.useIsotopes = useIsotopes; + end + end + + function[d18O] = runPSM(obj, X) + %% Runs a PRYSM cellulose PSM given data from a state vector ensemble + % + % d18O = obj.run(X) + % + % ----- Inputs ----- + % + % X: The data used to run the PSM. A numeric matrix with three + % rows. First row is temperature, second row is precipitation, + % third row is relative humidity. + % + % ----- Outputs ----- + % + % d18O: Estimates of d18O of cellulose. + + T = X(1,:); + P = X(2,:); + RH = X(3,:); + d18O = obj.run(T, P, RH, obj.d18Os, obj.d18Op, obj.d18Ov, obj.model, obj.useIsotopes); + end + end + +end \ No newline at end of file diff --git a/prysmCoral.m b/prysmCoral.m new file mode 100644 index 00000000..347bdb9b --- /dev/null +++ b/prysmCoral.m @@ -0,0 +1,161 @@ +classdef prysmCoral < PSM + %% Implements the coral module of the PRYSM package + % + % Prerequisites: Python 3.4, numpy, scipy, and rpy2. See the PRYSM + % documentation for more details. + % + % Find PRYSM on Github at: https://github.com/sylvia-dee/PRYSM + % + % Or read the paper: + % Dee, S. G., Russell, J. M., Morrill, C., Chen, Z., & Neary, A. + % (2018). PRYSM v2. 0: A proxy system model for lacustrine archives. + % Paleoceanography and Paleoclimatology, 33(11), 1250-1269. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2020 + + properties + lat; + lon; + useSSS; + species = "Default"; + bcoeffs = prysmCoral.b_default; + end + properties (Constant, Hidden) + b_default = [0.3007062, 0.2619054, 0.436509, 0.1552032, .15]; + end + + methods (Static) + function[coral] = run(lat, lon, SST, SSS, d18O, species, bcoeffs) + %% Runs the PRYSM coral module + % + % coral = prysmCoral.run(SST, SSS, d18O, lat, lon, species, bcoeffs) + % Runs the Prysm coral model given sea surface temperatures and + % either sea surface salinities or d18O (sea-water) for a site. + % + % ----- Inputs ----- + % + % Please see the documentation of PRYSM module + % "psm.coral.sensor.pseudocoral" for details on the inputs. + % + % ----- Outputs ----- + % + % coral: Coral estimates. + + % d18O placeholder if using SSS + if isempty(d18O) + d18O = -1; + end + + % B coefficients + b = prysmCoral.b_default; + if exist('bcoeffs','var') && ~isempty(bcoeffs) + dash.assertVectorTypeN(bcoeffs, 'numeric', [], 'bcoeffs'); + assert(numel(bcoeffs)<=5, 'bcoeffs cannot have more than 5 elements'); + for k = 1:numel(bcoeffs) + if ~isnan(bcoeffs(k)) + b(k) = bcoeffs(k); + end + end + end + + % Numpy arrays + SST = py.numpy.array(SST); + SSS = py.numpy.array(SSS); + d18O = py.numpy.array(d18O); + + % Run PSM, convert output to Matlab numeric + coral = py.psm.coral.sensor.pseudocoral(lat, lon, SST, SSS, d18O, species, b(1), b(2), b(3), b(4), b(5)); + coral = numeric(coral); + end + end + + methods + function[obj] = prysmCoral(rows, useSSS, lat, lon, species, bcoeffs, name) + %% Creates a new prysmCoral object + % + % obj = prysmCoral(rows, useSSS, lat, lon, species, bcoeffs, name) + % Creates a PRYSM coral PSM for sea surface temperature and + % either sea surface salinity or d18O (sea-water) for a site. + % + % obj = prysmCoral(..., name) + % Optionally name the PSM + % + % ----- Inputs ----- + % + % rows: The state vector rows with the temperature, + % precipitation, and relative humidity data for the proxy site. + % First row should be sea surface temperature, second should + % either be sea surface salinity or d18O (sea-water) + % + % useSSS: A scalar logical indicating whether the second set of + % state vector rows are for sea surface salinity (true) or + % d18O of sea-water (false) + % + % name: An optional name for the PSM. A string. + % + % Please see the documentation of PRYSM module + % "psm.coral.sensor.pseudocoral" for details on the + % remaining inputs. + % + % ----- Outputs ----- + % + % obj: The new prysmCoral object + + % Name, R estimation, rows + if ~exist('name', 'var') + name = ""; + end + obj@PSM(name, false); + obj = obj.useRows(rows, 2); + + % Model inputs + dash.assertScalarType(useSSS, 'useSSS', 'logical', 'logical'); + obj.useSSS = useSSS; + dash.assertScalarType(lat,'lat','numeric','numeric'); + obj.lat = lat; + dash.assertScalarType(lon,'lon','numeric','numeric'); + obj.lon = lon; + dash.assertStrFlag(species, 'species'); + allowed = ["Default","Porites_sp","Porites_lob","Porites_lut","Porites_aus","Montast","Diploas"]; + dash.checkStrsInList(species, allowed, 'species', 'species name'); + obj.species = species; + + % b coefficients + if exist('bcoeffs','var') && ~isempty(bcoeffs) + dash.assertVectorTypeN(bcoeffs, 'numeric', [], 'bcoeffs'); + assert(numel(bcoeffs)<=5, 'bcoeffs cannot have more than 5 elements'); + for k = 1:numel(bcoeffs) + if ~isnan(bcoeffs(k)) + obj.bcoeffs(k) = bcoeffs(k); + end + end + end + end + + function[d18O] = runPSM(obj, X) + %% Runs a PRYSM coral PSM given data from a state vector ensemble + % + % d18O = obj.run(X) + % + % ----- Inputs ----- + % + % X: The data used to run the PSM. A numeric matrix with two + % rows. First row is sea-surface temperatures, second row is + % either sea-surface salinity or d18O (sea-water). + % + % ----- Outputs ----- + % + % d18O: Estimates of d18O of cellulose. + + SST = X(1,:); + SSS = X(2,:); + d18O = []; + if ~obj.useSSS + d18O = SSS; + end + d18O = obj.run(obj.lat, obj.lon, SST, SSS, d18O, obj.species, obj.bcoeffs); + end + end + +end \ No newline at end of file diff --git a/prysmIcecore.m b/prysmIcecore.m new file mode 100644 index 00000000..1555f939 --- /dev/null +++ b/prysmIcecore.m @@ -0,0 +1,105 @@ +classdef prysmIcecore < PSM + %% Implements the icecore module of the PRYSM package + % + % Prerequisites: Python 3.4, numpy, scipy, and rpy2. See the PRYSM + % documentation for more details. + % + % Find PRYSM on Github at: https://github.com/sylvia-dee/PRYSM + % + % Or read the paper: + % Dee, S. G., Russell, J. M., Morrill, C., Chen, Z., & Neary, A. + % (2018). PRYSM v2. 0: A proxy system model for lacustrine archives. + % Paleoceanography and Paleoclimatology, 33(11), 1250-1269. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2020 + + properties + altitudeDifference = 0; + end + + methods (Static) + function[d18O] = run(d18O, altitudeDifference) + %% Runs the PRYSM icecore module + % + % d18O = prysmIcecore.run(d18O, altitudeDifference) + % Runs the Prysm ice core model given d18O (precipitation) for + % a site. + % + % ----- Inputs ----- + % + % Please see the documentation of PRYSM module + % "psm.icecore.sensor.icecore_sensor" for details on the + % inputs. + % + % ----- Outputs ----- + % + % d18O: Estimates of icecore d18O + + % Placeholder time input + time = 1:size(d18O,2); + + % Numpy arrays + d18O = py.numpy.array(d18O); + time = py.numpy.array(time); + + % Run the PSM, convert back to Matlab numeric type + d18O = py.psm.icecore.sensor.icecore_sensor(time, d18O, altitudeDifference); + d18O = numeric(d18O); + end + end + + methods + function[obj] = prysmIcecore(row, altitudeDifference, name) + %% Creates a new prysmIcecore object + % + % obj = prysmIcecore(row, altitudeDifference) + % Creates a PRYSM icecore PSM for d18O (precipitation) at a + % proxy site. + % + % obj = prysmIcecore(..., name) + % Optionally name the PSM + % + % ----- Inputs ----- + % + % row: The state vector row with d18O (precipitation) needed to + % run the PSM. + % + % name: An optional name for the PSM (a string) + % + % Please see the documentation of PRYSM module + % "psm.icecore.sensor.icecore_sensor" for details on the + % remaining inputs. + + % Name, R estimation, rows + if ~exist('name','var') + name = ""; + end + obj@PSM(name, false); + obj = obj.useRows(row, 1); + + % Error check, set model inputs + if exist('altitudeDifference','var') || isempty(altitudeDifference) + dash.assertScalarType(altitudeDifference, altitudeDifference, 'numeric', 'numeric'); + obj.altitudeDifference = altitudeDifference; + end + end + + function[d18O] = runPSM(obj, d18O) + %% Runs a PRYSM ice core PSM given data from a state vector ensemble + % + % d18O = obj.run(X) + % + % ----- Inputs ----- + % + % X: The d18O (precipitation) data used to run the PSM. A + % numeric row vector. + % + % ----- Outputs ----- + % + % d18O: Estimates of ice core d18O. + + d18O = obj.run(d18O, obj.altitudeDifference); + end + end +end \ No newline at end of file diff --git a/prysmSpeleothem.m b/prysmSpeleothem.m new file mode 100644 index 00000000..23153add --- /dev/null +++ b/prysmSpeleothem.m @@ -0,0 +1,131 @@ +classdef prysmSpeleothem < PSM + %% Implements the speleothem module of the PRYSM package + % + % Prerequisites: Python 3.4, numpy, scipy, and rpy2. See the PRYSM + % documentation for more details. + % + % Find PRYSM on Github at: https://github.com/sylvia-dee/PRYSM + % + % Or read the paper: + % Dee, S. G., Russell, J. M., Morrill, C., Chen, Z., & Neary, A. + % (2018). PRYSM v2. 0: A proxy system model for lacustrine archives. + % Paleoceanography and Paleoclimatology, 33(11), 1250-1269. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2020 + + + properties + timeStep = 1/12; + model = 'Adv-Disp'; + tau0 = 1; + Pe = 1; + end + + methods (Static) + function[d18O] = run(d18O, T, timeStep, model, tau0, Pe) + %% Runs the PRYSM speleothem module + % + % d18O = prysmSpeleothem.run(d18O, T, timeStep, model, tau0, Pe) + % Runs the Prysm cellulose model given d18O (dripwater) and + % temperature at a speleothem site. + % + % ----- Inputs ----- + % + % Please see the documentation of PRYSM module + % "psm.speleo.sensor.speleo_sensor" for details on the + % inputs. + % + % ----- Outputs ----- + % + % d18O: Estimates of speleothem d18O (calcite). + + % Placeholder for time input + time = 1:size(T,2); + + % Convert to numpy arrays + d18O = py.numpy.array(d18O); + T = py.numpy.array(T); + time = py.numpy.array(time); + + % Run the PSM, convert back to Matlab numeric type + d18O = py.psm.speleo.sensor.speleo_sensor(time, d18O, T, timeStep, model, tau0, Pe); + d18O = numeric(d18O); + end + end + + methods + function[obj] = prysmSpeleothem(rows, timeStep, model, tau0, Pe, name) + %% Creates a new prysmSpeleothem object + % + % obj = prysmSpeleothem(rows, timeStep, model, tau0, Pe) + % Creates a PRYSM speleothem PSM given dripwater d18O and + % temperature for a speleothem site. + % + % obj = prysmSpeleothem(..., name) + % Optionally name the PSM + % + % ----- Inputs ----- + % + % rows: The state vector rows with the d18O and temperature + % data. First row should be d18O, second row is temperature. + % + % name: An optional name for the PSM. A string. + % + % Please see the documentation of PRYSM module + % "psm.speleo.sensor.speleo_sensor" for details on the + % remaining inputs. + % + % ----- Outputs ----- + % + % obj: The new prysmSpeleothem object + + % Name, R estimation, rows + if ~exist('name','var') + name = ""; + end + obj@PSM(name, false); + obj = obj.useRows(rows, 2); + + % Error check and set model inputs + assert(isscalar(timeStep) && timeStep>0, 'timeStep must be a positive scalar'); + obj.timeStep = timeStep; + if exist('model','var') && ~isempty(model) + dash.assertStrFlag(model); + assert(any(strcmp(model, ["Well-Mixed","Adv-Disp"])), 'model can either be "Well-Mixed" or "Adv-Disp"'); + obj.model = model; + end + if exist('tau0', 'var') && ~isempty(tau0) + dash.assertScalarType(tau0, 'tau0', 'numeric', 'numeric'); + obj.tau0 = tau0; + end + if exist('Pe', 'var') && ~isempty(Pe) + dash.assertScalarType(Pe, 'Pe', 'numeric', 'numeric'); + obj.Pe = Pe; + end + if exist('timeStep','var') && ~isempty(timeStep) + assert(timeStep==1/12, 'timeStep must be 1/12 (monthly)'); + end + end + + function[d18O] = runPSM(obj, X) + %% Runs a PRYSM speleothem PSM given data from a state vector ensemble + % + % d18O = obj.run(X) + % + % ----- Inputs ----- + % + % X: The data used to run the PSM. A numeric matrix with three + % rows. First row is d18O (dripwater), second row is + % temperature. + % + % ----- Outputs ----- + % + % d18O: Estimates of speleothem d18O. + + d18O = X(1,:); + T = X(2,:); + d18O = obj.run(d18O, T, obj.timeStep, obj.model, obj.tau0, obj.Pe); + end + end +end \ No newline at end of file diff --git a/vslitePSM.m b/vslitePSM.m index 3fbaee5b..e8677b70 100644 --- a/vslitePSM.m +++ b/vslitePSM.m @@ -1,4 +1,16 @@ classdef vslitePSM < PSM + %% Implements the VS-Lite tree-ring PSM + % + % Find it on Github at: https://github.com/suztolwinskiward/VSLite + % + % Or read the paper: + % Tolwinski-Ward, S. E., Evans, M. N., Hughes, M. K., + % and Anchukaitis, K. J.: An efficient forward model of the climate + % controls on interannual variation in tree-ring width, Clim. Dynam., + % 36, 2419–2439, doi:10.1007/s00382-010-0945-5, 2011. + + % ----- Written By ----- + % Jonathan King, University of Arizona, 2019-2020 properties phi; @@ -11,6 +23,26 @@ methods (Static) function[Y] = run(phi, T1, T2, M1, M2, T, P, varargin) + %% Runs the VS-Lite forward model + % + % trw = vslitePSM.run(phi, T1, T2, M1, M2, T, P, varargin) + % Runs the VS-Lite PSM given monthly temperature and + % precipitation data at a site. + % + % ----- Inputs ----- + % + % T: A matrix of temperature data with 12 rows (one per month). + % First row is January + % + % P: A matrix of precipitation data with 12 rows (one per moth). + % First row is January. + % + % Please see the documentation on VS-Lite function + % "VSLite_v2_3" for details on the other inputs. + % + % ----- Outputs ----- + % + % trw: Estimates of tree ring width. % Error check vars = {phi, T1, T2, M1, M2}; @@ -30,15 +62,41 @@ end methods - % Constructor function[obj] = vslitePSM(rows, phi, T1, T2, M1, M2, options, name) + %% Creates a new vslitePSM object + % + % obj = vslitePSM(rows, phi, T1, T2, M1, M2, options) + % Creates a VS-Lite PSM for a set of temperature and + % precipitation values + % + % obj = vslitePSM(..., name) + % Optionally name the PSM + % + % ----- Inputs ----- + % + % rows: The state vector rows with the monthly temperature and + % precipitation data for the site. The first 12 rows are + % January-December temperature. Rows 13-24 are + % January-December precipitation. + % + % options: A cell vector with the "varargin" inputs for the + % VS-Lite function. + % + % name: An optional name for the PSM. A string. + % + % ----- Outputs ----- + % + % obj: A new vslitePSM object + + % Name, rows, etc if ~exist('name','var') || isempty(name) name = ""; end obj@PSM(name, false); - obj = obj.useRows(rows); + obj = obj.useRows(rows, 24); + % Save the parameters obj.phi = phi; obj.T1 = T1; obj.T2 = T2; @@ -47,12 +105,24 @@ obj.options = options; end - % function[Y] = runPSM(obj, X) + %% Runs the VS-Lite PSM given data from a state vector ensemble + % + % trw = obj.runPSM(X) + % + % ----- Inputs ----- + % + % X: State vector data used to run the PSM. A numeric matrix + % with 24 rows. Rows 1-12 are monthly temperatures from + % January to December. Rows 13-24 are monthly precipitations + % from January to December. + % + % ----- Outputs ----- + % + % trw: Estimates of tree-ring width. T = X(1:12, :); - P = X(1:12, :); - + P = X(13:24, :); Y = obj.run(obj.phi, obj.T1, obj.T2, obj.M1, obj.M2, T, P, obj.options{:}); end end