function dec = obtainDec(A, varargin)
% generates a decomposition for a fixed set of facet normals
%
% Syntax: dec = computeDecomposition(A, 'mapping', groupMappingMatrix)
%   or    dec = computeDecomposition(A, 'constraint', constraintMatrix)
%
% The returned decomposition structurure contains data that is required to
%   map invalid/unconfined vectors to the validity/confinement cone. It
%   also contains a partitioning of the validity/confinement cone for
%   efficient measure calculations. Please, use this <a href="matlab:help helpDec>decDetail">Link</a> for detailed 
%   information on used entities, or this <a href="matlab: help helpDec>decOutput">Link</a> for a description of the
%   output structure.
%
% The procedure has three steps:
%   1) Checking, if the mapping/contraint fixes the position of the
%      polytope. If this is not the case, the mapping/constraint is
%      corrected to fix the position, dependent on the parameter
%      allowFixPosition.
%   2) Obtaining facet validity information in H-representation and
%      checking whether the decomposition is proper. If it is improper, a
%      proper decomposition is constructed (dec.properData) and projected
%      accordingly by projectDec, dependent on the parameter allowImproper.
%   3) Obtaining the validity cone and the required data to map invalid
%      hC-vectors to valid hC-vectors without a change in shape.
%   4) Decompose the validity cone into unified and simplex partitions to
%      prepare the output for efficient measure calculations.
%   5) If the original representation was improper, the constructed proper
%      representation is projected to yield the improper representation.
%
% Input options:
%   verbose - changes verbosity level from 0 up to 2. Default: 0
%   zerotol - tolerance to determine if numbers are close to zero. Zerotol
%     is a relative tolerance, as entities are scaled. Default: elkZerotol
%   constraint - a constraint matrix C*h = 0. Default: []
%   mapping    - a mapping matrix h = M*hC. Default: []
%   allowFixPosition - allow modification of the mapping/constraint matrix
%     to fix the position of represented polytopes. A value (0) throws an
%     error, if position of input is not fixed. A value (1) throws a
%     warning, corrects the mapping/constraints and proceeds. A value (-1)
%     performs as (1), but without a warning. Default: 0
%   allowImproper - allow the processing of improper representations. A
%     value (0) throws an error for improper representations. A value (1)
%     throws a warning, but constructs a proper representation. A value
%     (-1) proceeds as (1), but without a warning. Default: 0
%
% See also: validateDec, computePartitionIndexDec, addMeasureDataDec,
%           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
% input parameters
par = mapOptionStruct(varargin, ...     
       'verbose', 2,...
       'zerotol', elkZerotol,...
       'allowFixPosition', 0,...
       'allowImproper', 0,...
       'constraint', [],...
       'mapping', []);
   
% normalize facet normals
for iFacet = 1:size(A, 2)
    A(iFacet, :) = A(iFacet, :) / norm(A(iFacet, :), 2);
end
dec.A = A;
   
% general information
dec.nDim = size(A, 2);
dec.nH = size(A, 1);
dec.inputData = par;


%% (1) Handle constraint and mapping input; fixPosition
if par.verbose, disp('Generate mapping and constraint matrices (includes fixing position) .. ');end
tic
dec.constraintMatrix = par.constraint;
dec.mappingReducedToFull = par.mapping;
dec = computeMappingAndConstraint(dec, par.zerotol, par.allowFixPosition);
dec.nC = size(dec.mappingReducedToFull, 2);
dec.timingData.fixPosition = toc;
if par.verbose, disp('done.');end

% Result: we now have a constraint matrix, both mapping matrices and all 
%   three are consistent. They also imply that the position of every 
%   polytope is fixed that fulfills the constraints or is generated by the
%   mapping.

%% (2) Cone of valid h/hC-vectors
% get the cone in h-space - this function only generates inequalities for
%   (nDim+1)-tuplets of the facet normals in A
if par.verbose, disp('Gather inequalities for validity cone .. ');end
tic
  validityDataOriginal = gatherValidityConeIneq(dec.A, par.zerotol, ...
      max(par.verbose-1,0)); 
  dec.timingData.generateCone = toc;
if par.verbose, disp('done.');end

% decide on proper/improper decomposition - this call contains the
%   error/warning handling for allowImproper - it does create a field
%   dec.improperData if input was improper.
tic
  dec = constructProperRepresentation(dec, validityDataOriginal, ...
      par.zerotol, par.verbose, par.allowImproper);
  dec.timingData.propernessHandling = toc;

% transfer validity cone to hC and simplify
if par.verbose, disp('Simplify cone representation and extract initial structuring elements ..');end
tic
  % transfer cone to hC-space
  validityDataOriginal.A = validityDataOriginal.A * dec.mappingReducedToFull;
  % remove zero-rows and trim rows
  validityDataOriginal.A = trimConstraintMatrix(validityDataOriginal.A, ...
      par.zerotol, 1);
  
  dec = addValidityCone(dec, validityDataOriginal.A, par.zerotol);
  % the validity cone is idedntical to the confinement cone (validate dec
  %   will look at the decomposition of the confinement cone.
  dec.confinementCone = dec.validityCone;
  
  dec.timingData.simplifyCone = toc;
if par.verbose, disp('done.');end

if par.verbose, disp('Adding facet validity and mapping data ..');end
tic
  dec = addFacetValidityData(dec, validityDataOriginal, par.zerotol, ...
      max(par.verbose-1,0));
  % validity mapping and confinement mapping is equivalent for proper
  %   reprentations so that even if the function is named
  %   addValidityMappingData, it actually adds a field
  %   confinementMappingData to be compatible with improper representations
  %   for which validity mapping requires a confinement mapping (no shape
  %   change) and a validity mapping (shape changes).
  dec = addValidityMappingData(dec, par.zerotol);
  dec.timingData.addBoundaryData = toc;
  if par.verbose, disp('done.');end

% Result: cone of valid polytopes including extreme rays and volume.
%   The starting set of structuring elements that are scaled to the
%   definition of the cone closure.


%% Decompose this cone
% This algorithm consists of three parts. decomposeCone handles the
%   decomposition of the validity cone in a generic way. This means, it
%   handles the selection of sample hC-vectors and the handling of the
%   remaining domain.
% To generate and evaluate a sample, it uses generateMinkowskiPartition.
%   This function stores all partition informations in a struct that is
%   accumulated by decomposeCone.
% Finally, the accumulated information is written appropriately to the
%   decomposition structure by addMinkowskiPartitionData.
if par.verbose, disp('Decompose cone of valid polytopes .. ');end
tic
  partitionData = decomposeCone(dec, dec.validityCone.A, ...
      @generateMinkowskiPartition, par.zerotol, max(par.verbose-1, 0));
  dec = addMinkowskiPartitionData(dec, partitionData, par.zerotol);
  dec.timingData.decomposeCone = toc;
  if par.verbose, disp('done.');end
  
%% Finalize results
% the following will always be the case (so far dec is a proper
%   representaion)
dec.isComplete = 1;
dec.isOrdinary = 1;

% now we have to distinguish between proper and improper representation
if isfield(dec, 'improperData')
    % extract and remove improper data
    improperData = dec.improperData;
    dec = rmfield(dec, 'improperData');
    % intermediate beautify this proper decomposition
    dec = cleanupDecStructure(dec);
    % project to improper decomposition
    dec = projectDec(dec, improperData.mappingImproperToProper, ...
        'verbose', par.verbose);    
    
else
    % only sort fields and ensure all fields are present
    dec = cleanupDecStructure(dec);
end

end