function [cosineTheta, C, modelIndices] = MDPS(B, Xs, RMSECs, omega, returnFormat)
%
%   Code originally authored by Robert Spiers, Idaho State University, July
%   2020
%
% SUMMARY:
%
% MATLAB Code to select models from the matrix B using the Model Diversity
% and Prediction Similarity (MDPS) model selection framework developed at
% Idaho State University by Robert Spiers and Dr. John Kalivas, doi: 
%   
%==========================================================================
%   INPUTS:
%==========================================================================
%
%       -> B (m x n): Matrix of models, where m is the number of models, 
%                       and n is the number of variables
%                       (wavelengths for most chemometric purposes)
%
%       -> Xs (k x n): Matrix containing the k spectra drawn from the
%                       secondary distribution. Should contain the target
%                       desirable predicted samples (transductive model
%                       selection), and other spectra from the secondary
%                       distribution if available
%               :: NOTE: Secondary spectra MUST be preprocessed according
%                           to how B was formed.
%                           Non-processed spectra interact incorrectly with
%                           B and thus the yhats will be incorrect also.
%
%       -> RMSECs (m x 1): Vector of RMSECs corresponding to the prediction
%                       error of each model in 'B' with respect to the
%                       calibration samples that the model was formed with
%                       (calibraiton samples with reference values)
%
%       -> omega (z x 1): Vector of weights for the model over-fitting and
%                       under-fitting penalty in the MDPS equation.
%                       Different models will be selected for each omega
%                       weight
%
%       -> returnFormat: String identifying whether the user prefers the
%                       return values in matrix format or vector format.
%                       Vector format is simpler for basic model selection,
%                       but both provide the same fundamental information,
%                       just with different post-processing.
%
%==========================================================================
%   OUTPUTS IF returnFormat == 'MATRIX'
%==========================================================================
%
%       -> cosineTheta (m x m): Symetric matrix of the cosine of the angle between
%                       all possible model combinations, where the column
%                       and row indices correspond to the two models
%                       compared (model diversity measure). Matrix has 1's
%                       on the diagonal
%
%       -> C (m x m x z): Matrix of the composite C values representing the
%                       fused prediction similarity measures. Incorporates
%                       the secondary prediction difference (SPD), average
%                       2-norm (for overfitting), and average RMSEC (for
%                       underfitting). Each layer on the third dimension
%                       corresponds to an omega value
%
%       -> modelIndices (m x m x 2): Indices of each of the models in each
%                       model combination. This matrix is functionally
%                       redundant since returnFormat == 'MATRIX'
%
%==========================================================================
%   OUTPUTS IF returnFormat == 'VECTOR'
%==========================================================================
%
%       DEFINE s = 0.5 * (s^2 - s), the number of possible combinations
%                                   of two models in B
%
%       -> cosineTheta (s x 1): Cosine of the angle between all possible 
%                       model combinations
%
%       -> C (s x z): Vector of the composite C values representing the
%                       fused prediction similarity measures. Incorporates
%                       the secondary prediction difference (SPD), average
%                       2-norm (for overfitting), and average RMSEC (for
%                       underfitting). Each layer on the second dimension
%                       corresponds to an omega value
%
%       -> modelIndices (s x 2): Indices of each of the models in every
%                       model combination. 
%
%==========================================================================

    % Verification of Correct Inputs
    
    switch nargin
        case 4
            returnFormat = 'matrix';
        case 3
            omega = 0.4;
    end
    
    if size(RMSECs, 1) == 1
        RMSECs = RMSECs';
    end
    
    if size(B, 1) ~= size(RMSECs, 1)
        if size(B, 2) ~= size(RMSECs, 1)
            error('Number of RMSECs and number of models does not match')
        else
            B = B';
        end
    end
    
    if size(B, 2) ~= size(Xs, 2)
        if size(B, 2) ~= size(Xs, 1)
            error('Number of variables does not match between B and Xs')
        else
            Xs = Xs';
        end
    end
    
    % Computing model diversity and prediction similarity measures

    m = size(B, 1);

    bnorms = vecnorm(B')';

    cosineTheta = (B * B') ./ (bnorms * bnorms');

    yhat = NaN(m, size(Xs, 1));

    for kk = 1:size(Xs, 1)
        yhat(:, kk) = Xs(kk, :) * B';
    end

    SPD = RangeScale(sum(abs(reshape(yhat, size(yhat, 1), 1, size(yhat, 2)) - reshape(yhat, 1, size(yhat, 1), size(yhat, 2))), 3));
    RMSECmat = RangeScale((RMSECs + RMSECs') / 2);
    bnormmat = RangeScale((bnorms + bnorms') / 2);

    C = NaN(m, m, length(omega));

    for ii = 1:length(omega)
        C(:, :, ii) = SPD + omega(ii) * (RMSECmat + bnormmat);
    end
    
    modelIndices = cat(3, repmat((1:m)', 1, m), repmat(1:m, m, 1));
    
    % Post-processing if user desires vector output
    
    if strcmpi(returnFormat, 'vector')
        nonRepeatIndices = logical(reshape(triu(ones(m, m), 1), [], 1));
        
        Cfinal = NaN(sum(nonRepeatIndices), length(omega));
        for ii = 1:length(omega)
            Ctmp = C(:, :, ii);
            Cfinal(:, ii) = Ctmp(nonRepeatIndices);
        end
        C = Cfinal;
        cosineTheta = cosineTheta(nonRepeatIndices);
        
        modelIDXfinal = NaN(sum(nonRepeatIndices), 2);
        for jj = 1:2
            modelIDXtmp = modelIndices(:, :, jj);
            modelIDXfinal(:, jj) = modelIDXtmp(nonRepeatIndices);
        end
        modelIndices = modelIDXfinal;
    end
end

function vrs = RangeScale(v)
    vrs = (v - min(v, [], 'all')) / (max(v, [], 'all') - min(v, [], 'all'));
end