function dec = addMeasureDataDec(dec, measure, par, varargin)
% add mixed volumes for measure calculation
%
% Syntax: dec = addMixedVolumeData(dec, measure)
%   or    dec = addMixedVolumeData(dec, measure, par)
%
% This function adds the coefficients that are required to calculate the
%   given measure based on a lambda-vector in a partition or a hC-vector.
%
% The parameters in par are usually directions and contain one parameter
%   for each row.
% 
% Use obtainMeasureInfo for a listing of measures that are supported.
%
% addMeasureDataDec(.., 'recalculate') enforces the algorithm to
%   recalculate the coefficients, if the measure is already present.
%
% See also: obtainDec, computeMeasureDec, validateMeasureData, validateDec

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 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/>.

%% Input Handling
[recalculateFlag verbose zerotol] = mapOption(varargin, ...
    'recalculate', 0, 'verbose', 0, 'zerotol', elkZerotol);

% make parameter (par) existing
if ~exist('par', 'var')
    par = [];
end

% collect measure information (name, handle, dim, ispar)
if isstruct(measure)
    measureStruct = measure;
else
    measureStruct = obtainMeasureInfo(measure, dec.A);
end

%% Create measure
if ~isfield(dec, 'measureData')
    dec.measureData = struct();
end

if ~isfield(dec.measureData, measureStruct.name) || recalculateFlag
    dec.measureData.(measureStruct.name) = measureStruct;
end

%% prepare seMatrix as h-vectors
hMatrix  = (dec.mappingReducedToFull*dec.se')';
%! DEV: This is not nice but this correction is necessary in the
%       scope of projectDec. And it is still valid in the scope of
%       crystal decompositions because the crystals are always
%       centered around zero.
% hMatrix(abs(hMatrix) < sqrt(zerotol)) = 0;
% OK - exactly this line causes a problem for decompositions with fixed
%   positions that WILL have facets close to zero. In particular, for
%   generic lactose, case 5, with the feret diameter in direction [1 1 1].
% Hence, this line is removed and measure calculation for projected
%   decompositions should rather rely on the proper embedding than on a new
%   compilation of the measure data.
%/

%% Handle tensor calculation
% Collect tensor
if ~measureStruct.isPar
    
    %dec.measureData.(measureStruct.name).par = [];
    measureStruct.par = [];
    
    if ~isfield(dec.measureData.(measureStruct.name), 'tensorSe')
        % compute tensor
        tic
        measureStruct.tensorSe = computeMixedVolumeTensor(...
            hMatrix, measureStruct.handleH, measureStruct.dim, ...
            dec.seIncidenceMatrix, verbose);
        measureStruct.time.mixedVolumeTensor = toc;   
        
        % convert tensor
        tic
        if verbose
            disp('Cycling through simplex partitions: ')
            fprintf('......');
        end
        for iPart = 1:length(dec.simplexPartition)
            if verbose
                fprintf('\b\b\b\b\b\b%5.1f%%', ...
                    (iPart/length(dec.simplexPartition))*100);
            end
            measureStruct.tensorHc{iPart} = convertMeasureTensorSe(...
                measureStruct.tensorSe, ...
                dec.simplexPartition{iPart}.seIndexVector, ...
                dec.simplexPartition{iPart}.mappingHcToLam, ...
                measureStruct.dim); 
        end
        if verbose, fprintf('\b\b\b\b\b\b');end
        measureStruct.time.convertedTensor = toc;
    end
    
    dec.measureData.(measureStruct.name) = measureStruct;
    
else
    
    % initialize tensorSe
    if ~isfield(dec.measureData.(measureStruct.name), 'tensorSe')
        dec.measureData.(measureStruct.name).tensorSe = cell(0);
    end
    
    % initialize par
    if ~isfield(dec.measureData.(measureStruct.name), 'par')
        dec.measureData.(measureStruct.name).par = [];
    end
    
    % shortcut
    measureStruct = dec.measureData.(measureStruct.name);
    
    % we handle parameters subsequently
    for iPar = 1:size(par, 1)
        % only so something, if measure for parameter does not yet exist
        if ~isempty(...
                eliminateRows(par(iPar,:), ...
                measureStruct.par, 0)...
                )
            
            % add parameter
            measureStruct.par(end+1, :) = par(iPar, :);
            
            % add seTensor
            tic
            measureStruct.tensorSe{end+1} = ...
                computeMixedVolumeTensor(hMatrix, ...
                    @(h)(measureStruct.handleH(h, par(iPar, :))), ...
                    measureStruct.dim, ...
                    dec.seIncidenceMatrix, verbose);
           measureStruct.time(iPar).mixedVolumeTensor = toc;   

           % add tensors for hC
           tic
           if verbose
            disp('Cycling through simplex partitions: ')
            fprintf('......');
           end
           for iPart = 1:length(dec.simplexPartition)
               if verbose
                   fprintf('\b\b\b\b\b\b%5.1f%%', ...
                       (iPart/length(dec.simplexPartition))*100);
               end
               measureStruct.tensorHc{iPar}{iPart} = convertMeasureTensorSe(...
                   measureStruct.tensorSe{end}, ...
                   dec.simplexPartition{iPart}.seIndexVector, ...
                   dec.simplexPartition{iPart}.mappingHcToLam, ...
                   measureStruct.dim);
           end
           if verbose, fprintf('\b\b\b\b\b\b');end
           measureStruct.time(iPar).convertedTensor = toc; 
        end
    end
    
    % store shortcut
    dec.measureData.(measureStruct.name) = measureStruct;
end

end

