function grassStruct = obtainGrassStruct(hcMatrix, par)
% collect information required for computeProjectionFromGrass
%
% 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/>.

%% Coordinate transformation
% We use PCA to align the main variances to the coordinate axes. The 
%   diagonal elements of S represent the variance along the
%   PCA-coordinates. While U contains the PCA coordinates as column
%   vectors.
if numel(par.manualScaling) == 1
    % singular value decomposition
    [U,S] = svd(hcMatrix, 0);
    % We want to represent the original points in positionMatrix such that 
    %   they fill a well defined large area. Therefore, we map them to the 
    %   scaled PCA-space, with:
    %     y = S^(-1) * U' * x
    %   The variance in each PCA component is now 1.
    % But beforehand S^(-1) is problematic for close-to-zero eigenvalues.
    %   So that we fix that manually:
    
    %! DEV: try to relax the scaling:
%     S = diag(diag(S));
%     S = 0.9*S + 0.1*eye(size(hcMatrix, 1));
    S = S.^par.manualScaling;
    %/
    
    S(abs(S(:)) < par.zerotol) = 0;
    Sinf = diag(1./diag(S));
    Sinf(isinf(Sinf(:))) = 0;
    % mapping to PCA space:
    grassStruct.mapToPcaSpace = Sinf * U';
    % To project coordinates in the scaled PCA-space back to original
    %   coordinates, we use:
    grassStruct.mapToRealSpace = U * diag(diag(S));
    %! DEV: is mapToPcaSpace required at all??
    %/
    % The initial range is now simply represented by the first targetDim
    %   coordinate axes in scaled space.
    grassStruct.baseRange = eye(size(hcMatrix, 1));
    grassStruct.baseRange(:, (par.targetDim+1):end) = [];
    % and the null space are the remaining coordinates
    grassStruct.baseNull = eye(size(hcMatrix, 1));
    grassStruct.baseNull(:, 1:par.targetDim) = [];
    
else
    % if the scaling matrix is already provided, we have:
    grassStruct.mapToPcaSpace = par.manualScaling;
    % the scaling might not have full rank (see above PCA scaling approach
    %   when the date itself is already singular); hence, we schould use
    %   pinv to be safe, here:
    grassStruct.mapToRealSpace = pinv(par.manualScaling, par.zerotol);
        
    % the initial subspace must now be selected from PCA projection and
    %   manually mapped to the scaled space
    pcaData = approxPca(hcMatrix, par);
    
    % Let y be a point in scaled-PCA space and x the corresponding point in
    %   the original real space. Let further denote xP and yP the corresponding
    %   projected points. We also have the original projection matrix P:
    %     xP = P * x
    %   and an unknown projection matrix Ps in scaled-PCA space:
    %     yP = Ps y
    %   This unknown projection matrix can be determined, by:
    %     xP = P *x
    %     mapToRealSpace yP = P * mapToRealSpace * y
    %     yP = mapToPcaSpace * P * mapToRealSpace * y
    %     Ps = mapToPcaSpace * P * mapToRealSpace
    scaledProjectionMatrix = grassStruct.mapToPcaSpace * ...
        pcaData.projectionMatrix * grassStruct.mapToRealSpace;
    
    % from these some base vectors for the range and null space
    %! DEV: this code was originally used, here but turned out to be
    %       numerically instable.
    grassStruct.baseRange = orth(scaledProjectionMatrix);
    grassStruct.baseNull = null(scaledProjectionMatrix);
    % The following is adopted from orth/null, but enforced to use the right
    %   dimensions
    % [U, S, V] = svd(initProjectionScaledPca, 0)
    % grassStruct.baseRange = U(:, 1:par.targetDim);
    % grassStruct.baseNull = V(:, (par.targetDim+1):end);
end

  
%! DEV: new approach using manual initial condition (deprecated): did only
%       shift the Grass space in directions but did not yield any
%       distortion.
%/ DEV: the below is very unfortunate for the PCA approach since the base
%       and null matrices are actually party of the unity matrix.
% grassStruct.baseRange = orth(grassStruct.mapToPcaSpace*par.manualInitialCondition);
% grassStruct.baseNull = null(grassStruct.baseRange');
%/

% grassmann dimensions
grassStruct.nDim = size(grassStruct.baseRange, 1);
grassStruct.targetDim = size(grassStruct.baseRange, 2);
grassStruct.grassDim = (grassStruct.nDim - grassStruct.targetDim) * ...
    grassStruct.targetDim;
grassStruct.zerotol = par.zerotol;


%% Data bounding cone
%! DEV: is this still required?? In high dimensions this part will take
%       extremely long. And so far, I think the nonlinear constraint wa
%       swtiched off, anyway.
% >> need a switch for turning on/off this feature - at least the
%    documentation requires the nonlincon to be built. But it is absolutely
%    awquard
% if isempty(par.manualBoundingCone)
%     
%     % these lines are seperated:
%     [grassStruct.extremePointMatrix, ...
%      grassStruct.coneNormalVectorMatrix] = ...
%         computeDataBoundingCone(grassStruct, dataOriginal.positionMatrix);
% 
% else
%     if isstruct(par.manualBoundingCone)
%         % hrep provided
%         grassStruct.coneNormalVectorMatrix = par.manualBoundingCone.A;
%         vrep = convertHrepToVrep(par.manualBoundingCone);
%         grassStruct.extremePointMatrix = vrep.V';
%     else
%         % generating rays provided
%         grassStruct.extremePointMatrix = par.manualBoundingCone;
%         hrep = convertVrepToHrep(struct('V', ...
%             [grassStruct.extremePointMatrix'; ...
%              zeros(grassStruct.nDim, 1)] ...
%             ));
%         grassStruct.coneNormalVectorMatrix = hrep.A;
%     end
end
%/
    
%! Debug
% figure(5)
% viewPolytope(struct('A', grassStruct.coneNormalVectorMatrix', ...
%                     'h', zeros(size(grassStruct.coneNormalVectorMatrix, 2), 1)), ...
%              'infinityBoxSize', 1)
% plot3(grassStruct.extremePointMatrix(1,:), ...
%       grassStruct.extremePointMatrix(2,:), ...
%       grassStruct.extremePointMatrix(3,:), 'rx', 'markerSize', 15);
% figure(12),hold on
% viewPolytope(struct('A', [coneNormalVectorMatrix'; 1 0 0], ...
%                     'h', ...
%      [zeros(size(grassStruct.coneNormalVectorMatrix, 2), 1); 1]), ...
%              'infinityBoxSize', 1)
% % 
% set(gca, 'PlotBoxAspectRatioMode', 'auto', ...
%          'DataAspectRatioMode', 'auto');
% xlim([-1.5 0.5]);ylim([-1 1]);zlim([-1 1]);
%/                              