function partitionData = generateConfinementPartition(dec, hC, zerotol)
% generate a unified partition with confinement condition
%
% THIS IS NO USER FUNCTION

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

%% initialize acuumulating array structure
if exist('hC', 'var')
    partitionData.sampleHc = hC;
else
    partitionData.sampleHc = [];
end
partitionData.flagResult = [];
partitionData.hrep = [];
partitionData.isValid = [];
partitionData.isConfined = [];
partitionData.unifiedIndex = [];
partitionData.validityProjectionMatrix = [];
partitionData.validityProjectionConstraintMatrix = [];
partitionData.violatedConstraintIndexVector = [];
% partitionData.confinementProjectionMatrix = [];
% partitionData.confinementConstraintMatrix = [];
partitionData.debug = [];
    
% quit if no input is given; decomposeCone uses this call to initialize
%   the accumulating array
if ~exist('dec', 'var')
    partitionData(:) = [];
    return
end

%% precalculation
% vector in proper embedding space
hE = dec.properData.mappingNewToProper * hC;

% projection data
projectionData = obtainProjectionDataForVector(dec.properData, hE, 0);
% save the validity projection for later measure mapping
partitionData.validityProjectionMatrix = projectionData.projectionMatrix;
partitionData.validityProjectionConstraintMatrix = ...
    projectionData.constraintMatrix;
% also save the violated constraints
partitionData.violatedConstraintIndexVector = ...
    projectionData.violatedConstraintIndexVector;

% what is the unified partition for hE' ?
unifiedIndex = computePartitionIndexDec(dec.properData, ...
    projectionData.projectionMatrix*hE, zerotol, 1);

%% Cases
% Everything that follows is still in hE-space. The transformation to
%   hC-space is equal for all cases.
if all((dec.properData.validityCone.A * hE) < 0)
    %% CASE:  hC is valid
%     disp('DBG: valid');
    % this will be a valid parittion
    partitionData.flagResult = 1;
    % but let's remember that its an originally valid one
    partitionData.isValid = 1;
    partitionData.isConfined = 1;
    
    % the hrep is the unified partition itself; projected to hC later
    thisA = dec.properData.unifiedPartition{unifiedIndex}.A;
    
    % we save the unified index
    partitionData.unifiedIndex = unifiedIndex;
     
elseif rank(projectionData.projectionMatrix * dec.properData.mappingNewToProper, ...
        zerotol) == dec.nC
    %% CASE: confined hC vector but not valid
%     disp('DBG: confined');
    % this will be a valid partition
    partitionData.flagResult = 1;
    % but not valid
    partitionData.isValid = 0;
    partitionData.isConfined = 1;
    
    % add the constraint such that the valid hE always resides in THIS
    %   unified partition.
    thisA = [projectionData.constraintMatrix; ...
             dec.properData.unifiedPartition{unifiedIndex}.A*...
             projectionData.projectionMatrix];  
         
    % we save the unified index
    partitionData.unifiedIndex = unifiedIndex;
    
else
    %% unconfined hC vector
%     disp('DBG: unconfined');
    % this partition will also be accepted, but is is neither valid nor
    %   confined. But we will store it together with the precalculated
    %   projection data.
    partitionData.flagResult = 1;
    partitionData.isValid = 0;
    partitionData.isConfined = 0;
    
    %! Dev: This was the code previously used. But since we now include
    %    unconfined partitions in the decomposition, we need to consider
    %    all constraints like for the confined partitions above.
    % removed part:
%     % add the constraint that the projection is valid, but do not store the
%     %   target unified partition as it is not relevant.
%     thisA = projectionData.constraintMatrix;
    % new part
    thisA = [projectionData.constraintMatrix; ...
             dec.properData.unifiedPartition{unifiedIndex}.A*...
             projectionData.projectionMatrix];
    % we save the unified index
    partitionData.unifiedIndex = unifiedIndex;
    %/
    
end

%% Final conclusion
% save original constraint matrix (in proper embedding)
partitionData.debug.embeddingConstraintMatrix = thisA;

% project constraints to new space
thisA = thisA*dec.properData.mappingNewToProper;
% while usually a lot of redundant facets will appear, especially
%   considering the transformation of the nonempty-cone, these are no
%   problem. Problematic are close-to zero projections where the facet
%   of the cone is almost parallel to the projection.
%
% remove zero rows and trim
thisA = trimConstraintMatrix(thisA, zerotol, 1);
% remove redundant rows
thisA = unique(thisA, 'rows');

% accumulate now:
partitionData.hrep.A = thisA;
partitionData.hrep.h = zeros(size(thisA, 1), 1);
partitionData.debug.projectionData = projectionData;


end

function [projectionMatrix constraintMatrix] = ...
    concatenateProjectionData(projectionData, filterVector)
% concatenate projection matrices according to a filter vector
%
% Input are the cell arrays of the single projection matrices and
%   constraints for this projection in hE-space.
%
% THIS IS NO USER FUNCTION
                          
% input
if ~exist('filterVector', 'var')
    filterVector = (1:length(projectionData.projectionList))>0;
end
properNc = size(projectionData.projectionList{1}, 1);

% init data
projectionMatrix = eye(properNc);
constraintMatrix = zeros(0, properNc);

% concatenate data
for iProj = find(filterVector)
    constraintMatrix = [constraintMatrix; ...
                        projectionData.constraintList{iProj}];
    projectionMatrix = projectionMatrix * ...
                       projectionData.projectionList{iProj}; 
end

end