function errorVector = computeMatrixError(matrixOne, matrixTwo, ...
    errorType, groupVector)
% compute error vector from matrices
%
% Syntax: errorVector = computeMatrixError(referenceMatrix, ...
%                         sampleMatrix, errorType)
%
% The following error types are supported, while a detailed definition can
%   be obtained <a href="matlab: help computeMatrixError>helpErrorDefinition">here</a>: relativeSquared, absoluteSquared, relativeMean, 
%   absoluteMean, relativeTotal, absoluteTotal.
%
% With errorType='all', all error vectors are returned in a structure with
%   fields according to the error types given above.
%
% See also: approximateDec

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

%% Info
% Note that weights for the elements might be applied before calling this
%   function. A weighted aggregation of all error components can be 
%   obtained by a subsequent scalar product between errorVectors and a 
%   weight vector.
%/

if ~exist('ratio', 'var') || isempty(ratio)
    ratio = 1;
end

n = size(matrixOne, 2);
nReal = n*ratio;
errorMatrix = matrixOne - matrixTwo;
switch lower(errorType)
    case 'relativesquared'
    errorVector = 1/nReal * sqrt(...
        sum((errorMatrix ./ matrixOne).^2, 2) );
    
    case 'relativemean'
    errorVector = 1/nReal * sum(abs(errorMatrix ./ matrixOne), 2);
    
    case 'relativemax'
    errorVector = max(abs(errorMatrix ./ matrixOne), [], 2);
    
    case 'absolutesquared'
    errorVector = (1/ratio) * sum(errorMatrix.^2, 2);

    case 'absolute'
    errorVector = (1/ratio) * sum(abs(errorMatrix), 2);
    
    case 'absolutemax'
    errorVector = max(abs(errorMatrix), [], 2);
    
    case 'relativetotal'
    errorVector = (1/ratio) * (sum(matrixOne, 2) - sum(matrixTwo, 2)) ./ ...
        sum(matrixOne, 2);
    
    case 'absolutetotal'
    errorVector = (1/ratio) * (sum(matrixOne, 2) - sum(matrixTwo, 2));
    
    case 'relativemaxgroupedmean'
    errorVector = 0*errorMatrix(:,1);
    for iGroup = 1:max(groupVector)
        errorVector = max(errorVector, abs((...
            sum(matrixOne(:, groupVector==iGroup), 2) - ...
            sum(matrixTwo(:, groupVector==iGroup), 2)) ./ ...
            sum(matrixOne(:, groupVector==iGroup), 2)));
    end
    
    case 'all'
    errorVector.relativeSquared = ...
        computeMatrixError(matrixOne, matrixTwo, 'relativeSquared');
    
    errorVector.relativeMean = ...
        computeMatrixError(matrixOne, matrixTwo, 'relativeMean');
    
    errorVector.absoluteSquared = ...
        computeMatrixError(matrixOne, matrixTwo, 'absoluteSquared');
    
    errorVector.absoluteMean = ...
        computeMatrixError(matrixOne, matrixTwo, 'absolute');
    
    errorVector.relativeTotal = ...
        computeMatrixError(matrixOne, matrixTwo, 'relativeTotal');
    
    errorVector.absoluteTotal = ...
        computeMatrixError(matrixOne, matrixTwo, 'absoluteTotal');
    
    otherwise
    error('elk:approximation:wrongInput', ...
        ['Given error type ''' errorType '''is not supported']);
    
end

end

function helpErrorDefinition()
% %%% Detailed Error definitions for <a href="matlab: help computeMatrixError">computeMatrixError</a>
% n  - total number of columns in matrixOne and matrixTwo
% xi - reference value (column of matrixOne)
% yi - sample value (column of matrixTwo)
% e  - aggregated error vector
%
% %% absoluteSquared
% e = sum( (xi - yi).^2 )
%
% %% relativeSquared
% e = 1/n * sqrt(sum( ((xi - yi) ./ xi )^2 ))
%
% %% absolute
% e = sum( abs(xi - yi) )
%
% %% relativeMean
% e = 1/n * sum( abs(xi - yi) )
%
% %% absoluteTotal
% e = ( sum( xi ) - sum( yi ) )
%
% %% relativeTotal
% e = ( sum( xi ) - sum( yi ) ) / sum( xi )
%
% %% 
end
