function M = computeMeasureDec(dec, vecType, vec, measureList, ...
    zerotol, useSimplex, enforceValid)
% calculate measure from decomposition data
%
% Syntax: measure = computeMeasureDec(dec, 'hc', hC, measure)
%   or    measure = computeMeasureDec(dec, 'se', {p, lam}, measure)
%
% computeMeasureDec(dec, 'hC', hC) calculates the measure based on the 
%   hC-vector of the polytope (Hc-representation). This vector must have
%   the length dec.nC.
%
% computeMeasureDec(dec, 'hC', {p, hC}) also provides the simplex
%   partitions for the vectors in hC.
%
% computeMeasureDec(SR, 'se', {p, lam}) calculates the measure bases on the 
%   vector lam that represents the linear combination of the structuring 
%   elements in the simplex partition that is given by p. The vector lam 
%   must have the length dec.nC in this case.
%
% The simple forms above are sufficient and simply calculate all measures
%   that are defined in dec.measureData.
% 
% getMeasure(.., measure, ..) specifies the measure name (the fields
%   in dec.measureData) and eventually the parameter par that is required.
%   If this list starts with 'grad' or 'gradient' the gradient of the
%   subsequent measures are calculated. Calculating measures and measure
%   gradients cannot be combined.
%   Examples: {'volume'}, {'volume', 'feret', [1 0 0;1 1 0]}, {'grad',
%             'surface'}, {'grad', 'volume', 'support', [1 0 0], 'feret',
%             [1 1 1; 0 0 1], 'facet', [1 0 0]}
%
% getMeasure(..., zerotol, simplex) lets you control the tolerance
%   (zerotol) that is used for algorithms (default elkZerotol) and enforces
%   the algorithm to search and use the proper simplex partition 
%   (simplex=1) instead of the unified partition (simplex=0). This also 
%   effects how p in {p, hc} is interpreted (default is simplex=0).
%
% See also: obtainMeasureInfo, validateMeasureDataDec, obtainDec,
%   validateDec

% NOT YET SUPPORTED AND, HENCE, NOT DOCUMENTED:
%
% Symbolic vectors hc or se can be used to obtain the polynomials that
%   correspond to the specified partitions (columns) and the specified 
%   measures (rows). The output is a symbolic matrix. For the input, it
%   should be a cell string {p, vec} where p is a vector of partition
%   indices and vec is a single symbolic vector.
% Status: draft

% 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/>.

%% make sure mex-file exists
try
    evaluateMeasurePolynome([], [], 1);
% if ~exist([mexPath 'evaluateMeasurePolynome.' mexext], 'file')
catch noError
    mexPath = mfilename('fullpath');
    mexPath = [mexPath(1:(end-17)) 'private' filesep];
    try
        compileMessage = evalc(['mex ' mexPath 'evaluateMeasurePolynome.c -outdir ' mexPath]);
    catch thisError
%         disp(compileMessage)
        rethrow(thisError)
    end
    try
        rehashMessage = evalc('rehash');
    catch thisError
        disp(rehashMessage)
        rethrow(thisError)
    end
end

%% Input handling
% set zerotol
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = elkZerotol;
end
% set simplex option
if ~exist('useSimplex', 'var') || isempty(useSimplex)
    useSimplex = 0;
end
% set enforceValid option
if ~exist('enforceValid', 'var') || isempty(enforceValid)
    enforceValid = 0;
end
% enforce measureName to cell
if ~iscell(measureList)
    measureList = {measureList};
end
% if input decomposition is not proper, we fall back to the proper
%   embedding
if isfield(dec, 'properData') && ~isfield(dec, 'measureData')
    if strcmpi(vecType, 'hc')
        if iscell(vec)
            warning('elk:decomposition:input', ['You have provided the ' ...
                'unified partition for an improper decomposition with no ' ...
                'measureData field. I compute the measure based on the ' ...
                'proper embedding, but ignore the provided partition data.']);
            vec = vec{2};
        end
        M = computeMeasureDec(dec.properData, vecType, ...
            dec.properData.mappingNewToProper*vec, ...
            measureList, zerotol, useSimplex, enforceValid);
        return
    else
        error('elk:decomposition:wrongInput', ['The input decomposition ' ...
            'does not contain measure data. Utilizing the proper embedding ' ...
            'is only supported for hC-vectors.']);
    end
end


[vecType, vec, partitionIndexVector] = clarifyVectorInformation(...
                                         lower(vecType), vec, useSimplex, dec, ...
                                         zerotol, enforceValid);
% the resulting partition index is always a simplex partition(!)
nVec =length(partitionIndexVector);

% set symbolic handler
isSymbolic =  strcmp(class(vec), 'sym');

%% strip gradient option
% if the measureList starts with 'grad' or 'gradient', the partial
%   derivative of all following measures is calculated.
if strcmpi(measureList{1}, 'grad') || strcmpi(measureList{1}, 'gradient')
    measureList(1) = [];
    gradientOption = 1;
else
    gradientOption = 0;
end

%% determine total number of measures
% the created vector has the number of measures assigned to each
%   position where a string would be in measName.
nMeasureVec = nan(1, length(measureList));
isParVec = nMeasureVec;
for iMeasure = 1:length(measureList)
    if ischar(measureList{iMeasure})
        nMeasureVec(iMeasure) = 1;
        isParVec(iMeasure) = 0;
    else
        nMeasureVec(iMeasure-1) = size(measureList{iMeasure}, 1);
        isParVec(iMeasure-1) = 1;
        nMeasureVec(iMeasure) = 0;
        isParVec(iMeasure) = nan;
    end
end
nMeasure = sum(nMeasureVec);

%% Initialize output
caseList = {'measure', 'gradient', 'symbolic measure', 'symbolic gradient'};
outputCase = caseList{1.0 + logical(gradientOption)*1 + logical(isSymbolic)*2};
switch outputCase
    case 'measure'
        M = zeros(nMeasure, nVec);
    case 'gradient'
        M = zeros(nMeasure, nVec, dec.nC);
    case 'symbolic measure'
        M = sym('tmp', [nMeasure, nVec])*0;
    case 'symbolic gradient'
        error('elk:decomposition:wrongInput', ...
              'Symbolic gradients for measure calculation are not supported.');
    otherwise
        error('elk:decomposition:internalError', ...
              'This should not happen at all.');
end

%% Prepare "easy-way-through (TM)"
% In the following there are numurous decisions to make that include
%   optional extra cycles:
%
% (1) cycle: measure names
% (2) cycle: parameters (if so)
% (3) cycle: gradients (if so)
% (4) cycle: vectors
%
% To simplify the loop and cross-decisions, we summarize the measure and
%   parameter loop. This needs some preparation:
measureIndexVector = [];
parameterIndexVector = [];
for iMeasure = find(nMeasureVec)
    measureIndexVector = [measureIndexVector ...
        iMeasure*ones(1, nMeasureVec(iMeasure))];
    parameterIndexVector = [parameterIndexVector 1:nMeasureVec(iMeasure)];
end

% cycle through measures/parameters AND simplex partition indices
if useSimplex
    partitionIndexCycleVector = 1:dec.nPS;
else
    partitionIndexCycleVector = dec.unifiedToSimplexVector(1:dec.nPU);
    % note that iPart is always a simplex partition (see below) while
    %   partitionIndexVector denotes either simplex or unified partitions.
    %   This is corrected, here.
end

for iMeasPar = 1:length(measureIndexVector)
for iPart = partitionIndexCycleVector    
    vectorFilter = (partitionIndexVector == iPart);
    
    % handle no vectors to deal with
    if ~any(vectorFilter)
        continue;
    end
    
    % handle invalid vectors (they have partitionIndex 0)
%     if iPart == 0
%         % there is nothing to do since all fields are prefilled with zeros.
%         continue;
%     end
    
    % determine tensor with/without parameter
    if isParVec(iMeasure)
        % obtain tensor (with parameter)
        [thisTensor tensorDim] = obtainTensorDec(dec, ...
            measureList{measureIndexVector(iMeasPar)}, ...
            measureList{measureIndexVector(iMeasPar)+1}(...
              parameterIndexVector(iMeasPar), : ...
            ), vecType, iPart);  
    else
        % obtain tensor (without parameter)
        [thisTensor tensorDim] = obtainTensorDec(dec, ...
            measureList{measureIndexVector(iMeasPar)}, [], ...
            vecType, iPart);
    end
    
    %! Debug: Here, I would need to distinguish symbolic and normal output
    %    But "M" must be initialized properly in that case. .. enforce this
    %    cast once, whenever a symbolic vector is met.
    %      > if [vector is symbolic] and [M not symbolic]
    %      >   then [cast to symbolic]
    %      > end
    % Hu?! Since vec is a matrix with all vectors, the symbolic thing is a
    %   single decision based on the input-vec only. This can already be
    %   considered with the initialization of M!!
    %
    % 1) Decide on symbolic/nonsymbolic and thisOutputCase in
    %    initialization step
    % 2) Handle calculations case-based like below
    % 3) That's it - Isn't it??
    %/
    switch outputCase
        case 'measure'
            M(iMeasPar, vectorFilter) = evaluateMeasurePolynome(...
                thisTensor, vec(:, vectorFilter), tensorDim);
        case 'gradient'
            for iGradDim = 1:dec.nC
            switch tensorDim
                case 1
                    M(iMeasPar, iVec, iGradDim) = thisTensor(iGradDim);
                case 2
                    tmpTensor = permute(thisTensor(iGradDim, :), [2 1]) ...
                              + thisTensor(:, iGradDim);
                    M(iMeasPar, vectorFilter, iGradDim) = ...
                        evaluateMeasurePolynome(...
                          tmpTensor, vec(:, vectorFilter), 1);
                case 3
                    tmpTensor = permute(thisTensor(iGradDim, :, :), [2 3 1]) ...
                              + permute(thisTensor(:, iGradDim, :), [1 3 2]) ...
                              + thisTensor(:, :, iGradDim);
                    M(iMeasPar, iVec, iGradDim) = evaluateMeasurePolynome(...
                          tmpTensor, vec(:, vectorFilter), 2);
            end
            end
        case 'symbolic measure'
            %! Dev: Obtaining symbolic expressions was implemented for
            %    shape approximation algorithms that used polynomial
            %    optimization. Since this approach did not work
            %    sufficiently well, the following code was not adopted.
            %    Minor changes would be required. The outer loop went over
            %    all vectors (iVec) before and is now running over all
            %    partitions (iPart).
            error('elk:decomposition:wrongInput', ['Extracting symbolic ' ...
                'expressions is currently not supported']);
%             zeroTensor = thisTensor*0;
%             thisTensor = sym('tmp', size(thisTensor(:)'))*0 + thisTensor(:)';
%             switch tensorDim
%                 case 1
%                     for iDim = (1:size(dec.se, 2))
%                         thisTensor(iDim) = thisTensor(iDim) * vec(iDim);
%                     end
%                 case 2
%                     for iDim = (1:size(dec.se, 2))
%                         thisFilter = zeroTensor;thisFilter(iDim,:) = 1;
%                         thisFilter = thisFilter(:) == 1;
%                         thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim);
%                         
%                         thisFilter = zeroTensor;thisFilter(:,iDim) = 1;
%                         thisFilter = thisFilter(:) == 1;
%                         thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim);
%                     end
%                 case 3
%                     for iDim = (1:size(dec.se, 2))
%                         thisFilter = zeroTensor;thisFilter(iDim,:,:) = 1;
%                         thisFilter = thisFilter(:) == 1;
%                         thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim);
%                         
%                         thisFilter = zeroTensor;thisFilter(:,iDim,:) = 1;
%                         thisFilter = thisFilter(:) == 1;
%                         thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim);
%                         
%                         thisFilter = zeroTensor;thisFilter(:,:,iDim) = 1;
%                         thisFilter = thisFilter(:) == 1;
%                         thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim);
%                     end
%             end
%             % the meaning in this symbolic case is more iPartition than iVec
%             M(iMeasPar, iVec) = sum(thisTensor);
        otherwise
            error('elk:decomposition:internalError', 'This should not happen.');
    end
end
end

end