function p = computePartitionIndexDec(dec, hcMatrix, zerotol, unified)
% determine the partition of a given hC-vector
%
% Syntax: p = computePartitionIndexDec(dec, hcMatrix)
%
% The simplex partition (the index to dec.simplexPartition) is returned so
%   that dec.simplexPartition{p}.A*hC <= 0 is fulfilled.
%
% getMeasure(..., zerotol, unified) lets you control the tolerance
%   (zerotol) that is used for algorithms (default elkZerotol) and lets the
%   function alternatively return the unified partition (the index to
%   dec.unifiedPartition) for unified=1 (defaulf unified=0).
%
% See also: obtainDec, validateDec, addMeasureDataDec, validateMeasureDec,
%   computeMeasureDec

% 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
% ensure hc being a column vector or matrix of hc-vectors
if size(hcMatrix, 1) ~= dec.nC
    if size(hcMatrix, 2) == dec.nC
        hcMatrix = hcMatrix';
    else
        error('elk:decomposition:wrongInput', ...
            ['the provided vector must be a constrained h-vector that does ' ...
             'have the dimension size(dec.se, 2)']);
    end
end

% set zerotol
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = elkZerotol*ones(1, size(hcMatrix,2));
elseif length(zerotol) < size(hcMatrix, 2)
    zerotol = zerotol*ones(1, size(hcMatrix,2));
end
% set unified option
if ~exist('unified', 'var') || isempty(unified)
    unified = 0;
end

% catch incomplete parititions
if ~dec.isComplete
    p = ones(1, size(hcMatrix, 2));
    if ~unified
        p = dec.unifiedToSimplexVector*p;
    end 
    return
end

%% gather validity/confinement information
% ensure that vector is in cone.. ..the confinement cone was
%   renamed according to the type of hC-representation to get
%   pin-pointed where the type makes a difference.
if dec.isComplete && dec.isOrdinary
    isConfinedVector = max(dec.confinementCone.A * hcMatrix, [], 1) < zerotol;
elseif dec.isComplete && ~dec.isOrdinary
    isConfinedVector = max(dec.confinementConeHull.A * hcMatrix, [], 1) < zerotol;
else
    error('elk:decomposition:internalError', 'This should not happen');
end

%% Determine partition
% the default is: "not in cone"
pU = zeros(1, size(hcMatrix, 2));
pS = zeros(1, size(hcMatrix, 2));
for iVector = 1:size(hcMatrix, 2)
    %! DEV: the partition search is definitely a candidate for MEX
    %       implementation. The numerical integrations in cs-solver with
    %       Vegas took about 0.33 seconds for the tensor computations while
    %       158 seconds (almost all time) went into the identification of
    %       the unified partition. It was a 7D paracetamol example that -
    %       indeed had a lot of partitions to search.
%     if mod(iVector, 100) == 0, disp(iVector/size(hc, 2));end
    %/
    
    if unified && isConfinedVector(iVector)
        % partition invalid if not found:
        pU(iVector) = -1;        
        for iPart = 1:length(dec.unifiedPartition)
            if all(dec.unifiedPartition{iPart}.A * hcMatrix(:, iVector) < zerotol(iVector))
                pU(iVector) = iPart;
                break;
            end
        end
    elseif isConfinedVector(iVector)
        % partition invalid if not found:
        pS(iVector) = -1;
        % find simplex partition
        for iPart = 1:length(dec.simplexPartition)
            if prod(double((dec.simplexPartition{iPart}.A * hcMatrix(:, iVector)) < zerotol(iVector)))
                pS(iVector) = iPart;
                break;
            end
        end
    end
    
end

% the following seems to be odd, but adds some safety to the current
%   implementation
if unified
    p = pU;
else
    p = pS;
end

% error handling
if any(p == 0)
    error('elk:decomposition:wrongInput', ...
        ['the given vector is outside of the confinement cone: ' ...
        ' the vector must fulfill (dec.validityCone.A * vec <= 0)']);
elseif any(p == -1) && dec.isOrdinary
    % this error could mean that the decomposition is not complete
    error('elk:decomposition:numericalProblem', ...
        ['the simplex partition of the given vector could not be found - ' ...
        'this should not happen.']);
elseif any(p == -1) && ~dec.isOrdinary
    warning('elk:decomposition:generalProblem', ...
        ['the unified/simplex partition of the given vector(s) could ' ...
        'not be found - the input decomposition is singular (not ' ...
        'ordinary) so that this can happen.']);
end

end