function dataApprox = approximateDec(heMatrix, dec, varargin)
% shape approximation based on sample points and decomposition data
%
% Syntax: dataApprox = approximateDec(hcMatrix, dec, 'option', value, ...)
%
% This function approximates the points in heMatrix (one sample point each 
%   column) by an orthogonal projection to a linear subspace of the given
%   Hc-representation (dec). The approximated hE vectors are computed by:
%       heApprox = P * hE
%   where the projection matrix is computed from
%       P = baseMatrix * baseMatrix'.
%   The base matrix contains nTargetDim columns that span the desired
%   linear subspace in the embedding hE-space. Both results are returned in
%   dataApprox as the fields
%     .projectionMatrix - this is P
%     .mappingReducedToEmbedding - this is baseMatrix
%     .mappingEmbeddingToreduced - this is pinv(baseMatrix) or simply
%        the baseMatrix' since the base vectors are orthonormal
%
% The method of the approximation is determined by: 
%       approximateDec(.., 'method', method)
%   where method='PCA' uses principal component analysis (a quick first 
%   result), 'nonlinear' uses a global nonlinear optimization (usually the
%   thing you wan't to do) and 'nonlinearLocal' finds only one local
%   optimum, starting from the PCA optimum.
%
% For the nonlinear approach, you can specify options a string cell array
%   of measures that shalle be consideres ('measureList', {'volume',
%   'meanWidth'}). A corresponding vector assigns the weights
%   ('measureWeightVector', [1 0]).
% 

%! ToDo: rewrite documentation


%
% The output structure contains
%   .time - computation time for approximation and conversion
%   .optimalGrassmanVector - defines the optimal projection subspace based
%      on the PCA subspace - different optimal projections might be compared
%      according to this data
%/
%
%   
% see also: obtainDec, obtainCrystalDec, addMeasureDec

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

%% input handling
par = mapOptionStruct(varargin, ...
    'method', 'PCA', ...
    'targetDim', 1, ...
    'nResample', 0, ...
    'measureList', 'volume', ...
    'measureWeightVector', [], ...
    'errorType', 'absoluteSquared', ...
    'errorGroupVector', [], ...
    'projectionType', 'true', ...
    'manualScaling', 0.5, ...
    'manualBoundingCone', [], ...
    'numTrialPoints', 500, ...
    'numStageOnePoints', 10, ...
    'zerotol', elkZerotol);

%% Scaling / Resampling
if par.nResample > 0
    disp('Reducing set of sample points..')
    resampleData = resampleHeMatrix(heMatrix, par.nResample);
    heMatrix = resampleData.heMatrixResampled;
    disp('  done.')
end

par = validateInput(heMatrix, dec, par);

% add measure values for original hE-vectors
% measure value matrix
par.measureValueMatrix = computeMeasureDec(dec, 'hC', heMatrix, ...
        par.measureList, par.zerotol, 0, 1);

%% handle approximation
tic
switch lower(par.method)
    case 'pca'
        dataApprox = approxPca(heMatrix, par);
    case 'nonlinear'
        dataApprox = approxNonlinear(heMatrix, dec, par);
    case 'nonlinearlocal'
        % call fmincon instead of global search based on fmincon to provide
        %   a local minimum (right thig to use this for testing)
        dataApprox = approxNonlinear(heMatrix, dec, par);
    otherwise
        error('elk:approximation', 'This should not happen.');
end

% dataApprox.mappingEmbeddingToReduced = orth(dataApprox.projectionMatrix)';
% dataApprox.mappingReducedToFull = pinv(dataApprox.mappingFullToReduced);
dataApprox.heMatrix = heMatrix;
if strcmpi(par.projectionType, 'true')
    dataApprox.heMatrixApprox = dataApprox.projectionMatrix*heMatrix;
else 
    dataApprox.heMatrixApprox = ...
        estimateHeMatrix(heMatrix, dataApprox.projectionMatrix, ...
                         dataApprox.nonlinearData.sampleProjectionData);
end
dataApprox.parameterStruct = par;
if exist('resampleData', 'var')
    dataApprox.resampleData = resampleData;
end
dataApprox.time = toc;
