function varargout = validateDec(dec, varargin)
% validate decomposition of valid cone
%
% Syntax: dec = validateDecomposition(dec)
%
% This function checks whether dec is a proper decomposition of the set of
%   facet normals in dec.A. The structuring elements are double-checked for
%   (constrained) indecomposability and the unified partitions are checked
%   against the cone of valid polytopes. The simplex partitions are not
%   verified as they are only used internally.
%
% The decomposition structure is extended with the fields:
%   dec.validation.tolerance - chosen zerotol
%   dec.validation.result - scalar logical 1 (passed) or 0 (failed)
%   dec.timing.validateDec - total time that was taken by validation
%
% validateDec(.., 'zerotol', zerotol) sets the zerotol. Default: elkZerotol
%
% validateDec(.., 'verbose', verbose) sets the verbosity. For verbose=0,
%   no messages are displayed, verbose=1 only states which tests went wrong
%   and verbose=2 also displays the succesful tests including some
%   statistical data.
%
% validateDec(.., 'output', out) specifies the output that is returned.
%   With out='dec' the decomposition structure is returned. With
%   out='result' only the test result is returned. With out='both' both
%   results are returned using [dec result] = validateDec(..).
%
% Quick checks that are redundant to the more time consuming ones:
%   (Q1) All structuring elements represent an extreme ray of at least one
%          unified partition. The algorithm stores indecomposable 
%          polytopes, even if the partition is rejected. A failure here
%          may result from the Minkowski decomposition (the polytope is
%          inproper) or from an incomplete decomposition (Q4 should fail
%          as well).
%   (Q2) All structuring elements are part of the cone of valid polytopes.
%   (Q3) No volume of a unified partition is zero. Together with (L2), this
%          ensures that the intersection of two unified partitions has at 
%          most dimension (nC-1).
%   (Q4) The sum of all volumes of the unified partitions equals the volume
%          of the valid cone. Together with (L2) this ensures that the 
%          decomposition of the valid cone is complete.
%
% Tests that take some more time:
%   (L1) All structuring elements are (constrained) indecomposable. This
%          double-checks the validity of the structuring elements. A
%          failure here maps to a failure in Minkowski decomposition.
%   (L2) All intersections of two unified partitions have the volume 0. 
%          Together with (Q3), this ensures that the intersection of two 
%          unified partitions has at most dimension (nC-1).
%   (L3) If the Minkowski decomposition for an inner point (not the
%          boundary) of a unified partition is computed. The resulting 
%          summands equal the extreme rays of the unified partition. The
%          partition was found this way. A failure here maps to a failure
%          in Minkowski decomposition.
%
% See also: obtainDec, validateMeasureDataDec

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

% timing
tic

%% Input
% options
opt = mapOptionStruct(varargin, ...
    'verbose', 2, ...
    'zerotol', elkZerotol, ...
    'output', 'undefined', ...
    'L1', 1);

if strcmpi(opt.output, 'undefined')
    if nargout == 0
        opt.output = 'none';
    elseif nargout == 1
        opt.output = 'dec';
    elseif nargout >= 2
        opt.output = 'both';
    end
end
        

% if no test fails, accept decomposition
acceptDec = 1;

% Get volumes in a vector
partitionVolumeVector = zeros(1, dec.nPU);
unifiedSeIndexMatrix = zeros(dec.nPU, dec.nS);
for iPartOne = 1:dec.nPU
    partitionVolumeVector(iPartOne) = dec.unifiedPartition{iPartOne}.volumeScaled;
    unifiedSeIndexMatrix(iPartOne, dec.unifiedPartition{iPartOne}.seIndexVector) = 1;
end

%% Output volume
if opt.verbose == 2
    disp(['Domain volume is:                   ' num2str(dec.confinementCone.volume)]);
    disp(['Sum of individual scaled partition volumes: ' ...
           num2str(sum(partitionVolumeVector))]);
    disp(['Minimal scaled partition volume is: ' ...
           num2str(min(partitionVolumeVector))]);
    disp(['Maximal scaled partition volume is: ' ...
           num2str(max(partitionVolumeVector))]);
end

%% Q1 - all structuring elements used
seUsedCount = sum(unifiedSeIndexMatrix, 1);
if any(seUsedCount == 0)
    if opt.verbose
        disp(['[-] Q1: The structuring element(s) ' ...
              num2str(find(seUsedCount == 0)) ...
              ' are no extreme ray of any partition']);
    end
    acceptDec = 0;
else
    if opt.verbose == 2
        disp('[+] Q1: All structuring elements are used for at least one partition');
    end
end

%% Q2 - all structuring elements are in the cone of valid polytopes
thisTestAccepted = 1;
for iSe = 1:dec.nS
    if all(dec.confinementCone.A * dec.seScaledToClose(iSe, :)' < opt.zerotol)
        % nothing to do, everything OK
    else
        if opt.verbose
            disp(['[-] Q2: The following structuring element is not in the valid cone:' ...
                  num2str(iSe)]);
        end
        thisTestAccepted = 0;
        acceptDec = 0;
    end  
end
% print something, if everything is alright
if opt.verbose == 2 && thisTestAccepted
    disp('[+] Q2: Every structuring element is contained in the valid cone');
end

%% Q3 - no zero volume of partitions
if any(isZero(partitionVolumeVector, opt.zerotol))
    if opt.verbose
        disp(['[-] Q3: The volume of the following unified partition(s) is zero: '...
              num2str(isZero(partitionVolumeVector, opt.zerotol))]);
    end
    acceptDec = 0;
else
    if opt.verbose == 2
        disp('[+] Q3: No partition volume equals zero');
    end
end

%% Q4 - sum of volumes
if isZero(sum(partitionVolumeVector) - 1, opt.zerotol)
    if opt.verbose == 2
        disp('[+] Q4: The sum of partition volumes equals the domain volume');
    end
else
    if opt.verbose
        disp('[-] Q4: The sum of partition volumes does not equal the domain volume');
        disp(['        ..the scaled difference (sum(volumes) - 1) is: ' ...
              num2str(sum(partitionVolumeVector)-1)]);
    end
    acceptDec = 0;
end

%% L1 - all structuring elements are indecomposable
if opt.L1
thisTestAccepted = 1;
thisHrep.A = dec.A;
for iSe = 1:dec.nS
    thisHrep.h = dec.mappingReducedToFull * dec.seScaledToClose(iSe, :)';
    summandMatrix = decomposeMinkowskiHrep(thisHrep, ...
        'constraint', dec.constraintMatrix, ...
        'result', 'indecomposable', 'matrix', 1, ...
        'zerotol', opt.zerotol);
    summandMatrix = (dec.mappingFullToReduced * summandMatrix')';
    summandMatrix = scaleToClose(summandMatrix, dec);
    if size(summandMatrix, 1) == 1 && ...
        all(isZero(summandMatrix - dec.seScaledToClose(iSe, :), opt.zerotol))
        % nothing to do, everything OK
    else
        if opt.verbose
            disp(['[-] L1: The following structuring element is not (constrained) indecomposable: ' ...
              num2str(iSe)]);
        end
        thisTestAccepted = 0;
        acceptDec = 0;
    end  
end
% print something, if everything is alright
if opt.verbose == 2 && thisTestAccepted
    disp('[+] L1: Every structuring element is (constrained) indecomposable');
end
end

%% L2 - All intersections of two unified partitions have the volume 0.
% totalIntersectionVolume = 0;
thisTestAccepted = 1;
for iPartOne = 1:dec.nPU
    hrepOne = closeCone(dec.unifiedPartition{iPartOne}, dec.centerRay);
    
    for iPartTwo = (iPartOne+1):dec.nPU
        
        hrepTwo = closeCone(dec.unifiedPartition{iPartTwo}, dec.centerRay);
        
        hrepIntersection = computeIntersectionHrep(hrepOne, hrepTwo, 0);
        
        hrepIntersectionRed = reduceHrepDimension(hrepIntersection, opt.zerotol);
        
%         volumeIntersection = computeVolumeHrep(hrepIntersection, zerotol);
        
%         if volumeIntersection > zerotol * dec.confinementCone.volume
        if size(hrepIntersectionRed.A, 2) == size(hrepIntersection.A, 2)
            disp(['[-] L2: Parition ' num2str(iPartOne) ' and ' num2str(iPartTwo) ' have a full-dimensional intersection'])
%             disp(['    ..the volume of this intersection is: ' ...
%                 num2str(volumeIntersection/dec.confinementCone.volume)])
%             disp(['    ..there individual volumes are ' ...
%                   num2str(partitionVolumeVector(iPartOne)) ...
%                   ' and ' ...
%                   num2str(partitionVolumeVector(iPartTwo)) ...
%                   ', repsectively']);
%             totalIntersectionVolume = totalIntersectionVolume + volumeIntersection;
            thisTestAccepted = 0;
            acceptDec = 0;
        else
            % nothing to do, everything is fine
        end
        
    end
end
if opt.verbose == 2 && thisTestAccepted
    disp('[+] L2: No two partitions have a full-dimensional intersection');
end

%% L3 - Minkowski decomposition for an inner point
thisTestAccepted = 1;
thisHrep.A = dec.A;
for iPart = 1:dec.nPU
    thisSummandMatrix = dec.seScaledToClose(...
        dec.unifiedPartition{iPart}.seIndexVector, :);
    
    thisHrep.h = dec.mappingReducedToFull * sum(thisSummandMatrix, 1)';
    
    summandMatrix = decomposeMinkowskiHrep(thisHrep, ...
        'constraint', dec.constraintMatrix, ...
        'result', 'indecomposable', 'matrix', 1, ...
        'zerotol', opt.zerotol);
    
    summandMatrix = (dec.mappingFullToReduced * summandMatrix')';
    
    summandMatrix = scaleToClose(summandMatrix, dec);
    
    if size(summandMatrix, 1) == size(thisSummandMatrix, 1) && ...
        isempty(eliminateRows(thisSummandMatrix, summandMatrix, opt.zerotol))
        % everything OK
    else
        if opt.verbose
            disp(['[-] L3: The Minkowski decomposition of an inner point of ' ...
                  'parition ' num2str(iPart) ' does not result in the ' ...
                  'extreme rays of that partition']);
        end
        
        thisTestAccepted = 0;
        acceptDec = 0;
    end
end
if opt.verbose == 2 && thisTestAccepted
    disp('[+] L3: Every decomposition of an inner point of a unified partition results in that partition');
end

%% Output
dec.validationData.tolerance = opt.zerotol;
dec.validationData.result = acceptDec;
dec.timingData.validateDec = toc;
switch lower(opt.output)
    case 'none'
        % nothing to do
    case 'dec'
        varargout{1} = dec;
    case 'result'
        varargout{1} = acceptDec;
    case 'both'
        varargout{1} = dec;
        varargout{2} = acceptDec;
    otherwise
        error('elk:decomposition:inputError', ...
            'The output option must be one of: ''dec'', ''result'', ''both''.');
end