function varargout = variance_analyser(varargin)
% function varargout = variance_analyser(varargin)
%
% This function calculates the cumulative variance for original and
% projected shapes. In this way we can somehow define how much of the
% variance of the data is captured after the projection. 
%
% Inputs:
%  DATA - a matrix NxM where the columns represent the dimensions and the
%         rows are samples (i.e. 100x2 for a 100 2D shapes).
%  EIGENVECTORS - the rows of this matrix correspond to the eigenvectors to
%                 use to project. The first column is associated with the highest
%                 eigenvalue. See PCA for details.
%  METHOD - either 'Euclidean' or 'Mahalanobis'. This determines how to
%           calculate the distance of each sample from the mean, i.e. its deviation.
%
% Outputs:
%  Vo - variance of the original data with respect to the mean
%  Vp - variance of the projected data with respect to the mean
%  Xp - the projected data
%
% Example:
%
% %Some zero centered 3D data
% X = randn(100,3);
% X(:,3) = X(:,1) + randn(size(X(:,1)))/50;
% X = X - ones(size(X, 1),1)*mean(X,1);
% 
% % Get the eigenvalues of the covariance matrix
% [V, L] = eig(cov(X));
% L = diag(L);
% [val, ind] = sort(L, 'descend');
% L = L(ind);
% V = V(:, ind);
% 
% % Take the 2 with most variance.
% P = V(:,1:2);
% 
% 
% [v_m_o, v_m_p, xmp] = variance_analyser('Data', X, 'eigenvectors', P, 'method', 'mahalanobis');
% [v_e_o, v_e_p, xep] = variance_analyser('Data', X, 'eigenvectors', P, 'method', 'euclidean');
% 
% v_m_r = v_m_p./v_m_o
% v_e_r = v_e_p./v_e_o
% figure(1); clf;
% hold on;
% plot_shape('pts', xmp, 'color', [0 0 1]);
% plot_shape('pts', X, 'color', [1 0 0]);
%
%
% Dr. A. I. Hanna (2005) CMP, UEA, 2006.
error(nargchk(0,inf,nargin));
P = [];
V = [];
method = 'euclidean';
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error('Arguments to variance_analyser must come param/value in pairs.')
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'data'
            X = varargin{i+1};
        case 'eigenvectors'
            P = varargin{i+1};
        case 'method'
            method = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to LISTDLG.  Name was ' varargin{i}])
    end
end
if isempty(X)
    error('Data parameter is required.')
end
if isempty(P)
    error('Eigenvector matrix is required.')
end
mu = mean(X,1);
Xo = X;
matmu = ones(size(X,1),1)*mu;
b = (P'*(X - matmu)')';
Xp = matmu + (P*b')';
switch lower(method)
    case 'euclidean'
        Vo = calc_cum_euclid_model_variance(Xo, mu);
        Vp = calc_cum_euclid_model_variance(Xp, mu);
    case 'mahalanobis'
        Vo = calc_cum_mahal_model_variance(Xo, mu);
        Vp = calc_cum_mahal_model_variance(Xp, mu);
end
varargout{1} = Vo;
varargout{2} = Vp;
varargout{3} = Xp;
return;
%%%%%
%
%
%%%%%
function V = calc_cum_euclid_model_variance(X, mu)
diff = X - ones(size(X,1),1)*mu;
sqdiff = diff.^2;
err = sum(sqdiff, 2);
V = sum(err)/(length(err)+1);
return;
%%%%%
%
%
%%%%%
function V = calc_cum_mahal_model_variance(X, mu)
X = X - ones(size(X,1),1)*mu;
C = cov(X);
% This is the (x-xm)*S*(x-xm)' bit
d = real(sum(X/C.*X,2));
V = sum(d)/(length(d)-1);
return;