clear all
%% Dynamic potash alum (cube -> octahedron)
dec = obtainCrystalDec('generic_potash_alum', 3);

%% project
mapNewToProper = [1 0; 0 1; 0.6 0.6];
% mapNewToProper = [dec.se(2,:)', dec.se(4,:)']
decProj = projectDec(dec, mapNewToProper);

%% random points in hC
hcMatrix = 2*rand(2, 1000);
heMatrix = decProj.properData.mappingNewToProper*hcMatrix;

%% COPY of validity mapping
heMatrixValid = heMatrix;
zerotol=elkZerotol;
% create filterConstraintMatrix
valueMatrix = decProj.properData.confinementMappingData.A * heMatrixValid;
[boundaryIndexVector hcIndexVector] = find(valueMatrix > zerotol);
invalidConstraintMatrix = [boundaryIndexVector, hcIndexVector];
invalidConstraintMatrix(:, 1) = dec.confinementMappingData.responsibleVector(...
    invalidConstraintMatrix(:, 1));

appliedProjectionMatrix = zeros(size(heMatrixValid, 2), decProj.properData.nC);
for iGroup = 1:decProj.properData.nC
    
    % these vectors are affected:
    thisVectorIndexVector = invalidConstraintMatrix(...
        invalidConstraintMatrix(:, 1) == iGroup, 2);
    % skip, if there is nothing to do
    if isempty(thisVectorIndexVector)
        continue
    end
    
    % the vectors might be chosen several times at once (multiple entries
    %   in responsibleMatrix after trafo to directions. This is corrected,
    %   here:
    [thisVectorIndexVector originalIndexVector] = unique(thisVectorIndexVector);
    
    % project from data (for increase of zerotol, see comments above)
    [heMatrixValid(:, thisVectorIndexVector) ...
     appliedProjectionMatrix(thisVectorIndexVector, iGroup)] = ...
          projectFromData(heMatrixValid(:, thisVectorIndexVector), ...
          decProj.properData.confinementMappingData.projectionData(iGroup), ...
          zerotol + 10*eps);
end
processMatrix = unique(appliedProjectionMatrix, 'rows')
hcMatrixValid = 100*ones(size(hcMatrix));
for iCase = 1:size(processMatrix, 1)
    thisVectorFilter = all(bsxfun(@eq, ...
        appliedProjectionMatrix, processMatrix(iCase, :)), 2)';
    thisProjection = eye(3);
    for iAppend = 1:size(processMatrix, 2)
        if processMatrix(iCase, iAppend)
            thisProjection = thisProjection*...
                decProj.properData.confinementMappingData.projectionData(...
                  iAppend).projectionList{processMatrix(iCase, iAppend)};
        end
    end
    
    hcMatrixValid(:, thisVectorFilter) = pinv(thisProjection * ...
        decProj.properData.mappingNewToProper)*heMatrixValid(:, thisVectorFilter);
end
heMatrixValidFromHc  =decProj.properData.mappingNewToProper*hcMatrixValid;

%% plot
close all;figure, hold on
% view embedding cone
viewPolytope(dec.validityCone, 'infinityBoxSize', 1.5)
% view projected cone
% viewPolytope(struct('V', 2*(mapNewToProper*([0 0;decProj.se]'))' ))
% view unified partitions
colorMatrix = [1 1 0;0 0 1; 1 0 0];
for iPart = 1:length(decProj.unifiedPartition)
    % plot meaning in hE-space
    viewPolytope(decProj.unifiedPartition{iPart}.embedding.hrepStrict, ...
        'color', colorMatrix(iPart, :), 'shade', 0.2, ...
        'infinityBoxSize', 3);
    % also plot the unified partitions as they are
    viewPolytope(struct('V', 4*(mapNewToProper*(...
        [0 0;decProj.se(decProj.unifiedPartition{iPart}.seIndexVector, :)...
        ]'))'), 'color', colorMatrix(iPart, :), 'shade', 0.4)
end

% % view trajectory
% timeVector = zeros(1, length(result));
% heMatrix = zeros(3, length(result));
% for iTime = 1:length(resultProj)
%     timeVector(iTime) = result(iTime).time;
%     heMatrix(:, iTime) = result(iTime).hcMatrix;
%     hcProjMatrix(:, iTime) = resultProj(iTime).hcMatrix;
% end
% heProjMatrix1 = decProj.properData.mappingNewToProper * hcProjMatrix;
% heProjMatrix2 = mapToProperEmbedding(hcProjMatrix, decProj, 1e-5);
% plot3(heMatrix(1,:), heMatrix(2,:), heMatrix(3, :), 'k.');
plot3(heMatrixValid(1,:), heMatrixValid(2,:), heMatrixValid(3, :), 'm.');
plot3(heMatrixValidFromHc(1,:), heMatrixValidFromHc(2,:), heMatrixValidFromHc(3, :), 'r.');
% plot3(heProjMatrix1(1,:), heProjMatrix1(2,:), heProjMatrix1(3, :), 'k.');
% plot3(heProjMatrix2(1,:), heProjMatrix2(2,:), heProjMatrix2(3, :), 'cx');
% % view trajectory for matching hC-vectors
% heMatrixProjectedHc = decProj.properData.mappingNewToProper*...
%     projectVectorDec(decProj, heMatrix);
% plot3(heMatrixProjectedHc(1,:), heMatrixProjectedHc(2,:), ...
%     heMatrixProjectedHc(3, :), 'm.-');
% sampleGrowthRate = 0.5*pbeDef.growthRate([1 1 1]',0,0,0);
% plot3([0 sampleGrowthRate(1)], [0 sampleGrowthRate(2)], ...
%     [0 sampleGrowthRate(3)], 'r-', 'lineWidth', 3)

%% more
% figure,viewCrystal(dec.cdef, result(1).hcMatrix)
% figure,viewCrystal(dec.cdef, heMatrixProjectedHc(:,1))
% figure,viewCrystal(dec.cdef, result(end).hcMatrix)
% figure,viewCrystal(dec.cdef, heMatrixProjectedHc(:,end))
% figure,viewCrystal(dec.cdef, result(12).hcMatrix)
% figure,viewCrystal(dec.cdef, heMatrixProjectedHc(:,12))