classdef PRISM_SRD < handle
%
%   Code originally authored by Robert Spiers, Idaho State University
%   Last Updated: June 21, 2022
%
% SUMMARY:
%
% MATLAB Code to quantitatively characterize similarity between a source
% set and target sample. Termed the Physicochemically Responsive Integrated
% Similarity Measure (PRISM). Protocol developed by Robert Spiers and Dr.
% John Kalivas. Publication can be found online.
%
%==========================================================================
%   USAGE:
%==========================================================================
%   
%       obj = PRISM_SRD(Xsrc, Ysrc, Xtrg, duplicateSamples, LVs, PCs);
%       obj.execute();
%
%       Note: First operation only instantiates the PRISM object. It is not
%       until the 'execute' operation is performed that the similarity
%       statistics are calculated.
%
%==========================================================================
%   INPUTS:
%==========================================================================
%
%       -> Xsrc (m x n): Matrix of source spectra, where m is the number of
%                       samples and n is the number of wavelengths.
%
%       -> Ysrc (m x 1): Vector of source analyte values corresponding to
%                       the spectra in Xsrc.
%
%       -> Xtrg (k x n): Matrix of target spectra where k is the number of
%                       targets. Each target spectrum is processed
%                       individually apart from the others. k=1 is valid.
%
%       -> duplicateSamples (scalar): Numeric flag indicating whether the
%                       target spectra are duplicate samples of the source
%                       spectra (e.g., same samples measured on different
%                       instruments or at different temperatures).
%               :: 0 indicates no duplicates
%               :: 1 indicates the full target set is duplicate of the full
%                       source set. They NEED to be in the same order also.
%
%       -> LVs (optional input): Specifies either which latent variables to
%               use or how the algorithm should determine how many latent 
%               variables to use
%               :: 'CVBJ' (default): Use a RMSEC, RMSECV, 2-norm,
%                       Jaggedness U-curve to find the latent variable at
%                       the base of the curve.
%               :: integer input: Use number of latent variables
%                       corresponding with the integer amount
%               :: Vector of integers: Use all the specified latent
%                       variables and apply data fusion to combine the
%                       latent measurements.
%
%       -> PCs (optional input): Specifies either which principal
%               components to use or how the algorithm should determine how
%               many principal components to use
%               :: Decimal less than 1: Number of principal components
%                       required to span the percent variance specified by 
%                       this input
%               :: 0.99 (default): Use 99% variance rule.
%               :: Vector of integers: Calculate similarity for all
%                       specified principal components, then fuse.
%
%==========================================================================
%   OUTPUTS
%==========================================================================
%
%       -> PRISMscore (k x 1): Vector containing the final PRISM scores for
%                       each of the input target samples.
%
%       -> PRISMx (k x 1): Vector containing the PRISM scores for only the
%                       spectral similarity merits.
%
%       -> PRISMy (k x 1): Vector containing the PRISM scores for only the
%                       model-informed similarity merits.
%
%       -> srcSimMatX (pX x m): Matrix containing the fused (but not 
%                       normalized) similarity measurements for the source
%                       samples compared to themselves. pX is the number of
%                       spectral similarity measures.
%
%       -> srcSimMatY (pY x m): Matrix containing the fused (but not 
%                       normalized) similarity measurements for the source
%                       samples compared to themselves. pY is the number of
%                       model-informed similarity measures.
%
%       -> trgSimMatX (pX x k): Matrix containing the fused (but not 
%                       normalized) similarity measurements for the target
%                       samples compared to the source set.
%
%       -> trgSimMatY (pY x k): Matrix containing the fused (but not 
%                       normalized) similarity measurements for the target
%                       samples compared to the source set.
%
%       -> Any property with 'Full' in the title corresponds to the
%       aforementioned 'SimMat' matrices but has the data before fusion.
%       They are structured as such:
%           :: Outer cell array ((m or k) x 1): Each sample left out
%           :: Inner matrices (p x 1 x (m-1 or m)): Sample left out 
%                       compared to each sample left in source (m-1) at
%                       each similarity measure (p) before unsigned and
%                       signed measures are evaluated.
%           :: Innermost cells: Similarity measurements for all the
%                       principal components or latent variables. Will be a
%                       scalar if measure does not require a component
%                       selection or there is only one PC/LV.
%
%==========================================================================
%   OTHER NOTES
%==========================================================================
%
%       -> Fusion via SRD is not the only option. Altering this code to
%       facilitate mean or median fusion is simple. Simply make alterations
%       to the .calculatePRISM() instance method.
%
%       -> The PRISM is implemented using object-oriented programming in
%       MATLAB. Although MATLAB is generally not used with object-oriented
%       programming, this is the easiest way to organize the static and
%       instance methods. If you are interested in altering the code to
%       suit your needs, you may wish to study the implementation of object
%       oriented programming in MATLAB first.
%
%       -> The PRISM should NOT be completely re-run for every new target
%       sample. .getSourceSimilarity() normalizes the measurements and only
%       requires source data. For a particular source dataset, this
%       operation only needs to happen once; this will drastically decrease
%       computational time. Use the following operations to do new target
%       samples when the obj.getSourceSimilarity() has already been run.
%               :: obj.resetSecondary(Xtrg);
%               :: obj.getTargetSimilarity();
%               :: obj.calculatePRISM();
%
%       -> PRISM is a subclass of the 'handle' superclass. This means that 
%       any instance of the PRISM is only a pointer to the data and is not
%       the data itself. This means that the data cannot be directly copied
%       over into a new variable: this would quickly overwhelm RAM storage.
%       Instead, one needs to use the obj.clone() method if they want to
%       copy over an instance without changing the original.

    properties
        Xsrc                    % Source spectra
        Ysrc                    % Source analyte values
        Xtrg                    % Target spectra
        duplicateSamples        % Flag indicating if there are duplicates
        LVs                     % Which latent variables to analyze
        PCs                     % Which PCs to analyze
        srcSimMatX              % See 'OUTPUTS' section for all of these
        srcSimMatY
        trgSimMatX
        trgSimMatY
        srcFullSimX
        srcFullSimY
        trgFullSimX
        trgFullSimY
        PRISMx
        PRISMy
        PRISMscore              % Final PRISM score indicating similarity
    end

    methods
        % Instantiation function for the PRISM
        function obj = PRISM_SRD(Xsrc, Ysrc, Xtrg, duplicateSamples, LVs, PCs)
            if nargin < 6
                PCs = 0.99;         % Default to 99% rule
            end
            if nargin < 5
                LVs = 'CVBJ';       % Default to CVBJ U-curve
            end
            if nargin < 4
                warning('Unspecified whether the target samples are duplicates of the source samples. Algorithm forced to assume.');
                if size(Xsrc, 1) == size(Xtrg, 1)
                    duplicateSamples = 1;
                    warning('Algorithm assumes the samples are duplicates.');
                else
                    duplicateSamples = 0;
                    warning('Algorithm assumes there are not duplicate samples.');
                end
            end
            if nargin < 3
                error('Insufficient input arguments. Please input source spectra, source analyte values, and target spectra');
            end

            obj.Xsrc = Xsrc;
            obj.Ysrc = Ysrc;
            obj.Xtrg = Xtrg;
            obj.duplicateSamples = duplicateSamples;
            
            % Parsing optional inputs for LVs and PCs
            if ischar(LVs)
                LVs = PRISM.getLVs(Xsrc, Ysrc, LVs);
            end
            obj.LVs = LVs;

            if PCs < 1
                PCs = 1:PRISM.getPCs(Xsrc, PCs);
            end
            obj.PCs = PCs;
        end

        function execute(obj)
            % Executes similarity analysis for the PRISM
            obj.getSourceSimilarity();
            obj.getTargetSimilarity();
            obj.calculatePRISM();
        end
        
        function getSourceSimilarity(obj)
            % Calculate the source to source similarity measurements
            [nXsim, nYsim, ~, ~] = PRISM.getNumSimilarity();
            simsX = NaN(nXsim, size(obj.Xsrc, 1));
            simsY = NaN(nYsim, size(obj.Ysrc, 1));
            fullSimsX = cell(size(obj.Xsrc, 1), 1);
            fullSimsY = cell(size(obj.Ysrc, 1), 1);

            Xs = obj.Xsrc;
            Ys = obj.Ysrc;
            PC = obj.PCs;
            LV = obj.LVs;
            
            parfor m = 1:size(obj.Xsrc, 1)
%                 fprintf('    Primary sample #%i \n', m);

                % Remove a source sample from the source set
                Xsrctmp = Xs;
                Ysrctmp = Ys;
                Xsrctmp(m, :) = [];
                Ysrctmp(m, :) = [];
                
                % Calculate similarity measures to rest of source set
                simX = PRISM.getXSimilarity(Xsrctmp, Xs(m, :), PC);
                simY = PRISM.getYSimilarity(Xsrctmp, Ysrctmp, Xs(m, :), LV);
                
                % Apply the fusion rules
                simsX(:, m) = PRISM.applyFusion(simX);
                simsY(:, m) = PRISM.applyFusion(simY);
                
                fullSimsX(m) = {simX};
                fullSimsY(m) = {simY};
            end
            
            obj.srcSimMatX = simsX;
            obj.srcSimMatY = simsY;
            obj.srcFullSimX = fullSimsX;
            obj.srcFullSimY = fullSimsY;
        end

        function getTargetSimilarity(obj)
            % Calculate the target to source similarity measurements
            [nXsim, nYsim, nXpre, nYpre] = PRISM.getNumSimilarity;
            
            simsX = NaN(nXsim, size(obj.Xtrg, 1));
            simsY = NaN(nYsim, size(obj.Xtrg, 1));
            fullSimsX = cell(size(obj.Xtrg, 1), 1);
            fullSimsY = cell(size(obj.Xtrg, 1), 1);
            
            % Different protocol whether the target sample is the same as
            % one of the source samples
            if obj.duplicateSamples == 0
                % If not duplicates, then all target samples can execute
                % simultaneously (efficient)
                simX = PRISM.getXSimilarity(obj.Xsrc, obj.Xtrg, obj.PCs);
                simY = PRISM.getYSimilarity(obj.Xsrc, obj.Ysrc, obj.Xtrg, obj.LVs);
            else
                % Otherwise, needs to execute iteratively since the source
                % set must change every iteration
                simX = cell(nXpre, size(obj.Xtrg, 1), size(obj.Xsrc, 1) - 1);
                simY = cell(nYpre, size(obj.Xtrg, 1), size(obj.Xsrc, 1) - 1);

                Xs = obj.Xsrc;
                Ys = obj.Ysrc;
                PC = obj.PCs;
                LV = obj.LVs;
                Xt = obj.Xtrg;
                parfor m = 1:size(obj.Xtrg, 1)
%                     fprintf('    Secondary sample #%i \n', m);
                    Xstmp = Xs;
                    Ystmp = Ys;
                    Xstmp(m, :) = [];      % Remove the sample from the source distribution
                    Ystmp(m, :) = [];
                    
                    simX(:, m, :) = PRISM.getXSimilarity(Xstmp, Xt(m, :), PC);
                    simY(:, m, :) = PRISM.getYSimilarity(Xstmp, Ystmp, Xt(m, :), LV);
                end
            end
            
            % Apply fusion rules to each target sample
            for m = 1:size(obj.Xtrg, 1)
                simsX(:, m) = PRISM.applyFusion(squeeze(simX(:, m, :)));
                simsY(:, m) = PRISM.applyFusion(squeeze(simY(:, m, :)));
                
                fullSimsX(m) = {simX(:, m, :)};
                fullSimsY(m) = {simY(:, m, :)};
            end
            
            obj.trgSimMatX = simsX;
            obj.trgSimMatY = simsY;
            obj.trgFullSimX = fullSimsX;
            obj.trgFullSimY = fullSimsY;
        end

        function calculatePRISM(obj)
            % Combine all similarity measurements into the final PRISM score

            % Calculate the PRISMx

            obj.PRISMx = NaN(size(obj.Xtrg, 1), 1);
            for m = 1:size(obj.Xtrg, 1)
                if obj.duplicateSamples == 0
                    simMat = [obj.srcSimMatX, obj.trgSimMatX(:, m)];        % Append the source fused data with the target fused data
                else
                    simSrc = PRISM.purgeSourceMatrix(obj.srcFullSimX, m);   % If there are duplicate samples, need to remove sample-in measurements from the source set
                    simMat = [simSrc, obj.trgSimMatX(:, m)];
                end

                for n = 1:size(simMat, 1)
                    simMat(n, :) = simMat(n, :) / norm(simMat(n, :));       % Row normalization to get all similarity measurents on same scale
                end

                SRDout = SRDlightweight(simMat, 'max');                     % SRD fusion rule. Could replace this with sum or median fusion if desired.
                ranks = SRDout.normSRD;

                % ranks = -1 * mean(simMat);                                % Un-comment this code if mean fusion is desired

                % obj.PRISMx(m) = (ranks(end) - mean(ranks(1:end-1))) / std(ranks(1:end-1));
                obj.PRISMx(m) = (ranks(end) - mean(ranks)) / std(ranks);    % z-score calculated including the target sample, otherwise z-values highly inflated.
            end

            % Calculate the PRISMy using identical protocol

            obj.PRISMy = NaN(size(obj.Xtrg, 1), 1);
            for m = 1:size(obj.Xtrg, 1)
                if obj.duplicateSamples == 0
                    simMat = [obj.srcSimMatY, obj.trgSimMatY(:, m)];
                else
                    simSrc = PRISM.purgeSourceMatrix(obj.srcFullSimY, m);
                    simMat = [simSrc, obj.trgSimMatY(:, m)];
                end

                for n = 1:size(simMat, 1)
                    simMat(n, :) = simMat(n, :) / norm(simMat(n, :));
                end

                SRDout = SRDlightweight(simMat, 'max');
                ranks = SRDout.normSRD;

                % ranks = -1 * mean(simMat);                                % Un-comment this code if mean fusion is desired

%                 obj.PRISMy(m) = (ranks(end) - mean(ranks(1:end-1))) / std(ranks(1:end-1));
                obj.PRISMy(m) = (ranks(end) - mean(ranks)) / std(ranks);
            end

            % Use Stouffer's z-score combination rule to fuse the z-scores.

            obj.PRISMscore = (obj.PRISMx + obj.PRISMy) / sqrt(2);
        end

        function newPRISM = clone(obj)
            % Clone method to produce a new PRISM instance from the existing
            % one. This is necessary since PRISM is a subclass of the handle
            % superclass. In other words, newPRISM = obj; will NOT work because
            % then the data for obj will be changed when newPRISM is changed.
            newPRISM = PRISM(obj.Xsrc, obj.Ysrc, obj.Xtrg, obj.duplicateSamples, obj.LVs, obj.PCs);     % Create new PRISM instance
            props = properties(obj);
            for m = 1:length(props)
                newPRISM.(props{m}) = obj.(props{m});           % Copy all properties
            end
        end

        function PE = getPredictionError(obj, Yt)
            % Convenience function to calculate the prediction error of the
            % target samples when using the source model. Should only be used
            % for evaluation processes
            
            % Yt is the target analyte values
            Xs = obj.Xsrc;
            Ys = obj.Ysrc;
            Xt = obj.Xtrg;

            LV = obj.getLVs(Xs, Ys, 'CVBJ');

            if obj.duplicateSamples
                PE = obj.PEduplicate(Xs, Ys, Xt, Yt, LV);
                return;
            end

            b = PLS1ErikB(Xs - mean(Xs), Ys - mean(Ys), LV);
            yhatt = (Xt - mean(Xs)) * b(:, LV) + mean(Ys);
            PE = abs(yhatt - Yt);
        end

        function resetSecondary(obj, Xtrg)
            obj.Xtrg = Xtrg;
        end

    end

    methods (Static)        % General Static Utility Methods
        function [nX, nY, nXb4, nYb4] = getNumSimilarity
            % getNumSimilarity stores the number of similarity measures for
            % each situation. To be used by other methods.
            % Stores total number of X and Y measures, as well as the number of
            % X and Y measures before (b4) signed and unsigned averages.

            nX = 26;
            nY = 18;
            nXb4 = 13;
            nYb4 = 10;
        end

        function PE = PEduplicate(Xs, Ys, Xt, Yt, LV)
            % Prediction error needs to be assessed differently if there are
            % duplicate samples in the source and target sets. This method
            % removes the source duplicate sample before modeling.

            yhatt = NaN(size(Yt));
            for m = 1:size(Xt, 1)
                Xstmp = Xs;
                Ystmp = Ys;
                Xstmp(m, :) = [];
                Ystmp(m, :) = [];

                b = PLS1ErikB(Xstmp - mean(Xstmp), Ystmp - mean(Ystmp), LV);
                yhatt(m) = (Xt(m, :) - mean(Xstmp)) * b(:, LV) + mean(Ystmp);
            end
            PE = abs(yhatt - Yt);
        end

        function PCs = getPCs(X, percent)
            % Use the percent variance rule to calculate the number of
            % principal components to use.
            % Notably, the data is first mean-centered, which means more
            % PCs must be spanned before percent is fulfilled.
            [~, S, ~] = svd(X - mean(X));
            PCs = find(cumsum(diag(S)) > (percent * sum(diag(S))), 1);
        end

        function LVs = getLVs(X, Y, selMeth)
            % Use a U-curve on the source set to determine the optimal
            % number of latent variables
            maxLV = PRISM.convergencePLS(X, Y);     % Determine how many latent variables are required before convergence in the 2-norms is observed

            splitNum = floor(0.80 * size(X, 1));        % 80/20 data split into calibration set and cross validation set    
            numSplits = 100;                            % 100 random data splits

            RMSEC = NaN(numSplits, maxLV - 1);
            RMSECV = NaN(numSplits, maxLV - 1);
            bnorms = NaN(numSplits, maxLV - 1);
            Js = NaN(numSplits, maxLV - 1);

            warning('off', 'all');

            for n = 1:numSplits
                rng = randperm(size(X, 1));
                Xcal = X(rng(1:splitNum), :);
                ycal = Y(rng(1:splitNum), :);
                Xcv = X(rng((splitNum + 1):end), :);
                ycv = Y(rng((splitNum + 1):end), :);

                B = PLS1ErikB(Xcal - mean(Xcal), ycal - mean(ycal), maxLV - 1);

                for k = 1:size(B, 2)
                    b = B(:, k);
                    yhatc = (Xcal - mean(Xcal)) * b + mean(ycal);
                    yhatcv = (Xcv - mean(Xcal)) * b + mean(ycal);

                    RMSEC(n, k) = sqrt(sum((yhatc - ycal).^2) / length(ycal));
                    RMSECV(n, k) = sqrt(sum((yhatcv - ycv).^2) / length(ycv));
                    bnorms(n, k) = norm(b);
                    Js(n, k) = sqrt(mean((diff(b)).^2));
                end
            end

            warning('on', 'all');

            RangeScale = @(v) (v - min(v)) / (max(v) - min(v));

            RMSEC = RangeScale(mean(RMSEC));
            RMSECV = RangeScale(mean(RMSECV));
            bnorms = RangeScale(mean(bnorms));
            Js = RangeScale(mean(Js));

            Ucurve = zeros(1, length(RMSEC));

            for c = 1:length(selMeth)                       % Put the relevant terms into the U-curve
                switch selMeth(c)
                    case 'C'
                        Ucurve = Ucurve + RMSEC;
                    case 'V'
                        Ucurve = Ucurve + RMSECV;
                    case 'B'
                        Ucurve = Ucurve + bnorms;
                    case 'J'
                        Ucurve = Ucurve + Js;
                end
            end

            [~, LVs] = min(Ucurve, [], 'omitnan');          % Find the minimum of the U-curve

        end

        function maxLV = convergencePLS(X, Y)
            % Find the number of latent variables required before 2-norm
            % convergence. This algorithm is described and justified in the
            % supporting information section of doi.org/10.1021/acs.jcim.0c01493
            bhats = PLS1ErikB(X - mean(X), Y - mean(Y), size(X, 1) - 1);
            difnorm = diff(vecnorm(bhats));
            difnorm = (difnorm - min(difnorm)) / (max(difnorm) - min(difnorm));

            maxLV = find(difnorm >= 0.01, 1, 'last');
        end
    end

    methods (Static)        % PRISM-specific static methods (compute similarity, fusion, etc.)
        function fused = applyFusion(sim, whichSigned)
            % Apply fusion rules to data structure containing similarity
            % measurements.
            % Namely: average over PCs and LVs, then average over source
            % samples in, but do these averages signed and unsigned
            % depending on which merit is analyzed.
            [~, ~, ~, nYbasic] = PRISM.getNumSimilarity;
            nSim = size(sim, 1);
            
            if nargin < 2
                if nSim == nYbasic
                    whichSigned = 1:8;          % For the y-measures, the last two merits are not to be signed
                else
                    whichSigned = 1:nSim;
                end
            end

            
            meaned = NaN(size(sim, 1), size(sim, 3));
            for j = 1:nSim
                for m = 1:size(sim, 3)
                    meaned(j, m) = mean(sim{j, 1, m});          % Average over the PCs or LVs
                end
            end

            fused = NaN(size(sim, 1) + length(whichSigned), 1);
            
            fused(1:length(whichSigned), :) = abs(mean(meaned(whichSigned, :), 2));     % Signed
            fused((length(whichSigned) + 1):end, :) = mean(abs(meaned), 2);             % Unsigned
        end

        function fused = purgeSourceMatrix(sim, idx)
            % Removes the source sample which is equivalent to the target
            % sample from the source data comparison matrix
            % It is included in both the source-sample-out and the
            % source-sample-in, need to purge from both

            simnu = sim(1:length(sim) ~= idx);          % Remove the source-sample-out
            numSourceLeft = length(simnu);
            fusedcell = cell(numSourceLeft, 1);         % Remove the source-sample-in
            for m = 1:numSourceLeft
                if m < idx
                    fusedcell(m) = {PRISM.applyFusion(simnu{m}(:, :, 1:numSourceLeft ~= (m - 1)))};     % The sample is already removed once so this counting is necessary
                else
                    fusedcell(m) = {PRISM.applyFusion(simnu{m}(:, :, 1:numSourceLeft ~= m))};
                end
            end
            fused = cell2mat(fusedcell');
        end

        function sim = getXSimilarity(Xs, Xt, PCs)
            % Calculate the spectral similarity merits for a given source
            % and target (or pseudo-target/source-sample-out) set.
            % Parallelized and optimized for maximum efficiency.
            [~, ~, nSim, ~] = PRISM.getNumSimilarity;
            sim = cell(nSim, size(Xt, 1), size(Xs, 1));

            xc = mean(Xs)';

            Xi = [Xs; Xt];
            
            simCent = cell(nSim, size(Xi, 1));
            parfor m = 1:size(Xi, 1)
                simTmp = cell(nSim, 1);
                xi = Xi(m, :)';

                simTmp(1) = {(xi' * xc) / (norm(xi) * norm(xc))};       % Cosine theta
                simTmp(2) = {norm(xi - xc)};                            % Euclidean distance
                simTmp(3) = {det([xi'; xc'] * [xi, xc])};               % Determinant

                simTmp(4) = {PRISM.EISC(xi, xc)};                       % EISC-centroid     
                simTmp(5) = {PRISM.EISC(xc, xi)};                       % EISC-target

                simTmp(6) = {PRISM.UCPAF(xi, xc)};                      % Unconstrained procrustes analysis

                simTmp(7) = {PRISM.CPAmat(xi, xc)};                     % Constrained procrustes analysis

                simTmp(8) = {(xi - xc)' * xc / (norm(xc) * norm(xi - xc))};     % Centered cosine omega
                simTmp(9) = {PRISM.MahalanobisDistance(Xs, xi, PCs)};           % Mahalanobis distance

                [r1, D] = PRISM.IPC(Xs, xi, PCs);                       
                simTmp(10) = {r1};                                      % Inner product correlation
                simTmp(11) = {D};                                       % Divergence criterion

                [Q, s] = PRISM.Qresidual(Xs, xi, PCs);      
                simTmp(12) = {Q};                                       % Q-residual
                simTmp(13) = {s};                                       % Sin phi

                simCent(:, m) = simTmp;
            end

            for i = 1:size(Xs, 1)
                for j = 1:size(Xt, 1)
                    jprime = j + size(Xs, 1);
                    for k = 1:nSim
                        sim(k, j, i) = {simCent{k, i} - simCent{k, jprime}};        % Execute the delta-similarity between target sample-out and source sample-in
                    end
                end
            end

        end

        function sim = getYSimilarity(Xs, Ys, Xt, LVs)
            % Calculate the analyte-informed similarity merits for a given 
            % source and target (or pseudo-target/source-sample-out) set.
            % Parallelized and optimized for maximum efficiency, although
            % not all measures fit the delta-similarity framework.

            [~, ~, ~, nSim] = PRISM.getNumSimilarity;
            sim = cell(nSim, size(Xt, 1), size(Xs, 1));

            B = PLS1ErikB(Xs - mean(Xs), Ys - mean(Ys), max(LVs));

            Xi = [Xs; Xt];
            
            % The first five similarity measures are best achieved without
            % using delta-similarity because merits 1 and 4 require Ys.

            yhats = (Xs - mean(Xs)) * B(:, LVs) + mean(Ys);
            yhatt = (Xt - mean(Xs)) * B(:, LVs) + mean(Ys);

            simLV = NaN(5, size(Xt, 1), size(Xs, 1), size(yhats, 2));

            for k = 1:size(yhats, 2)
                simLV(1, :, :, k) = (Ys').*((1./yhats(:, k))' - 1./yhatt(:, k));
                simLV(2, :, :, k) = (2./yhats(:, k))' - 2./yhatt(:, k);
                simLV(3, :, :, k) = yhats(:, k)' - yhatt(:, k);
                simLV(4, :, :, k) = Ys' - yhatt(:, k);
                simLV(5, :, :, k) = (yhats(:, k)'./vecnorm(Xs')) - (yhatt(:, k)'./vecnorm(Xt'))';
            end

            sim(1:5, :, :) = num2cell(simLV, 4);

            % Loop through the remaining similarity measures as delta-SMs

            B = B(:, LVs);

            simCent = cell(nSim - 5, size(Xi, 1));
            for m = 1:size(Xi, 1)
                simLV = NaN(nSim - 5, length(LVs));
                xi = Xi(m, :)';
                
                % Local sensitivity analysis (SA) merits

                eps = [0.01, 0.05, 0.1, 0.5, 5];
                del = eps * (max(Ys) - min(Ys));
                yhateps = NaN(length(del), length(LVs));
                yhatnegeps = NaN(length(del), length(LVs));
                for k = 1:length(LVs)
                    b = B(:, k);
                    yhati = (xi' - mean(Xs)) * b + mean(Ys);
                    for i = 1:length(del)
                        tmpY = [Ys; yhati + del(i)];
                        tmpX = [Xs; xi'];
                        Btmp = PLS1ErikB(tmpX - mean(tmpX), tmpY - mean(tmpY), LVs(k));
                        yhateps(i, k) = (xi' - mean(tmpX)) * Btmp(:, LVs(k)) + mean(tmpY);

                        tmpY = [Ys; yhati - del(i)];
                        Btmp = PLS1ErikB(tmpX - mean(tmpX), tmpY - mean(tmpY), LVs(k));
                        yhatnegeps(i, k) = (xi' - mean(tmpX)) * Btmp(:, LVs(k)) + mean(tmpY);
                    end

                    simLV(1, k) = mean(yhateps - yhatnegeps);                                   % SAvar
                    simLV(2, k) = mean(abs(yhateps - yhati) + abs(yhatnegeps - yhati)) / 2;     % SAbias-a
                    simLV(3, k) = mean(yhateps - yhati + yhatnegeps - yhati) / 2;               % SAbias-s

                    projmult = (eye(length(b)) - b * b') / (b' * b);
                    simLV(4, k) = norm(projmult * (xi - mean(Xs)'));                            % Orthogonally projected Euclidean distance
                    simLV(5, k) = norm(projmult * (xi - mean(Xs)')) / norm(xi - mean(Xs)');     % Sin psi
                end

                simCent(:, m) = num2cell(simLV, 2);
            end

            for i = 1:size(Xs, 1)
                for j = 1:size(Xt, 1)
                    jprime = j + size(Xs, 1);
                    for k = 1:size(simCent, 1)
                        sim(k + 5, j, i) = {simCent{k, i} - simCent{k, jprime}};        % Delta similarity for the last five merits
                    end
                end
            end


        end

    end

    methods (Static)        % Similarity merit equations
        function [bd] = EISC(x1, x2)
            % Calculation of extended inverted scatter correction
            % Correction term based off of x2

            warning('off', 'all');                                              % Inverse not always well conditioned so warnings off
            len = length(x1) - 2;
            Xc = [ones(len, 1), x2(1:len), x2(1:len).^2, x2(1:len).^3, diff(x2(1:(len + 1))), diff(x2, 2)];     % Calculate the correction matrix
            ivXc = (Xc' * Xc) \ Xc';                                                                            % Invert the correction matrix using Moore-penrose inverse
            bd = norm(ivXc * (x1(1:len, :) - x2(1:len, :)));                    
            warning('on', 'all');
        end

        function F = UCPAF(x1, x2)
            % Unconstrained procrustes analysis
            % Based on SVDs of outer products and therefore they are rank 1
            F = norm((x2 * x2' / (x2' * x2)^2) * (x1 * x1' - x2 * x2'), 'fro');
        end

        function [H] = CPAmat(x1, x2)
            % Constrained procrustes analysis
            % Closed form solution for this particular setup
            H = norm((x2 * x2' / (x2' * x2)) * ((norm(x2) / norm(x1)) * (x1 * x1') / (x2' * x1) - 1), 'fro');
        end

        function [MD] = MahalanobisDistance(Xr, x1, PCs)
            % Mahalanobis distance
            mnxr = mean(Xr);
            [U, S, V] = svds(((Xr - mnxr)' * (Xr - mnxr)) / (size(Xr, 1)), max(PCs));
            MD = NaN(length(PCs), 1);
            for j = 1:length(PCs)
                Cinv = V(:, 1:PCs(j)) / S(1:PCs(j), 1:PCs(j)) * U(:, 1:PCs(j))';
                MD(j) = sqrt((x1' - mnxr) * Cinv * (x1' - mnxr)');
            end
        end

        function [r1, Div] = IPC(Xr, x2, PCs)
            % Inner product correlation and divergence criterion
            % Combined together since they use the same SVD

            X2 = x2 * x2';
            X2inv = (x2 * x2') / (norm(x2)^2);
            [U, S, V] = svds(Xr' * Xr / (size(Xr, 1) - 1), max(PCs));

            r1 = NaN(length(PCs), 1);
            Div = NaN(length(PCs), 1);
            for j = 1:length(PCs)
                r1(j) = 1 - (trace(X2 * (U(:, 1:PCs(j)) * S(1:PCs(j), 1:PCs(j)) * V(:, 1:PCs(j))')) / sqrt(norm(x2)^(4) * trace(S(1:PCs(j), 1:PCs(j)).^(2))));

                Div(j) = abs(0.5 * trace((X2 - U(:, 1:PCs(j)) * S(1:PCs(j), 1:PCs(j)) * V(:, 1:PCs(j))') * (X2inv - V(:, 1:PCs(j)) / S(1:PCs(j), 1:PCs(j)) * U(:, 1:PCs(j))')) + ...
                    0.5 * trace((X2inv + V(:, 1:PCs(j)) / S(1:PCs(j), 1:PCs(j)) * U(:, 1:PCs(j))') * (x2' - mean(Xr))' * (x2' - mean(Xr))));
            end
        end

        function [Q, sine] = Qresidual(Xr, x1, PCs)
            % Q residual and sine phi
            % Combined together since they use the same SVD

            [~, ~, V] = svds(Xr, max(PCs));
            Q = NaN(length(PCs), 1);
            sine = NaN(length(PCs), 1);
            for j = 1:length(PCs)
                xpreproj = (eye(length(x1)) - (V(:, 1:PCs(j)) * V(:, 1:PCs(j))'));
                x1proj = xpreproj * x1;
                Q(j) = (norm(x1proj))^2;
                sine(j) = norm(x1proj) / norm(x1);
            end
        end

    end

end
 %---------------------------------------------------------------------------------------
 function [out] = SRDlightweight(input, target)
%
%   Code originally authored by Robert Spiers, Idaho State University
%   Last Updated: June 21, 2022
%
% SUMMARY:
%
% A lightweight version of the original SRD code in order to decrease the
% computational time associated with LAFR. Removes bells and whistles, 
% such as the cross validation, fusion, and special cases on the CRRN. 
% Original theory presented in:
    % Heberger, K. Sum of ranking differences compares methods or models
    % fairly. Trends in Analytical Chemistry. 2010, 29, 101-109. 
% Comparison of Ranks by random numbers is also used, as presented in: 
    % Heberger, K. Sum of ranking differences for method discrimination and 
    % its validation: comparison of ranks with random numbers. 
    % Journal of Chemometrics. 2011, 25, 151-158.
%   
%==========================================================================
%   INPUTS:
%==========================================================================
%
%       -> input: SRD input matrix, arranged with the merits being each
%                   row, and the spectra are columnwise. Thus, mxn, with
%                   m being number of similarity merits, and n being number
%                   of samples.
%
%       -> target: SRD target string, specified by any of: 'max', 'min',
%                   'median', or 'mean'
%
%==========================================================================
%   OUTPUTS:
%==========================================================================
%
%       -> out: Struct file containing fields of output information
%
%           :: .mnCRRN: Mean rank of the CRRN distribution
%           :: .stdCRRN: Standard deviation of the CRRN distribution
%           :: .srdRank: Raw SRD rankings of the input samples
%           :: .input: Input matrix 'input'
%           :: .target: SRD target 'target'
%
%==========================================================================

    calculateCRRN = false;

    % Determine Target Vector
    if ischar(target)
        switch target
            case 'min'
                T = min(input, [], 2);
            case 'max'
                T = max(input, [], 2);
            case 'median'
                T = median(input, 2);
            case 'mean'
                T = mean(input, 2);
        end
    else
        T = target;
    end
    
    % Generate SRD Ranking
    
    [~, ~, rT] = unique(T);
    rI = NaN(size(input));
    for m = 1:size(input, 2)
        [~, ~, rI(:, m)] = unique(input(:, m));
    end
    
    srdRANK = sum(abs(rT - rI));

    mx = srdMax(size(input, 1));
    
    if calculateCRRN
        % Determine CRRN Distribution
        CRRN = NaN(10000, 1);
        for n = 1:10000
            CRRN(n, 1) = sum(abs(randperm(length(rT)) - (1:length(rT))));
        end
        out.mnCRRN = mean(CRRN);
        out.stdCRRN = std(CRRN);
        out.normCRRN = CRRN / mx;
    end
    
    out.srdRank = srdRANK;
    out.input = input;
    out.target = target;
    out.normSRD = srdRANK / mx;
end

function m = srdMax(nrow)
    % Given  the number of rows, compute the maximum SRD
    %==========================================================================
    if rem(nrow,2) == 1
        k = (nrow - 1) / 2;
        m = 2 * k * (k + 1);
    else
        k = nrow / 2;
        m = 2 * k^2;
    end
    %==========================================================================
end
