clear all

%% startup
dec = obtainCrystalDec('generic_potash_alum', 1);
dec = addMeasureDataDec(dec, 'volume');
dec = addMeasureDataDec(dec, 'surface');
dec = addMeasureDataDec(dec, 'meanWidth');
decExt = computeUnifyingDec(dec);
dec = addMeasureDataDec(dec, 'support', decExt.A);
decExt = computeUnifyingDec(dec);
extendedSupportDirections = decExt.A([1 9 15], :);

%% write pbeDef
pbeDef.dec = dec;
% pbeDef.seeding.positionMatrix = mapToValidityCone([1 1.4 10]', dec, 0) - [0 0 3e-10]';
pbeDef.seeding.positionMatrix = [1 0.75]';
pbeDef.seeding.numberVector = 1;
pbeDef.growthRate = @(H, P, y, t)(bsxfun(@plus, 0*H, ...
    [0.5 1.25]'));
pbeDef.timeSpan = 0:0.01:1;
% solve
result = solvePbe(pbeDef, 'relTol', 1e-5, 'distanceLimit', 1e-3);
% read ot to single matrix
dataOriginal = readResultData(result, pbeDef);
heMatrixOriginal = mapToExtendedEmbedding(dataOriginal.hcMatrix, decExt);

%% approximation
% scaling
heMatrixScaled = bsxfun(@times, heMatrixOriginal, ...
    1./ computeMeasureDec(dec, 'hC', dataOriginal.hcMatrix, 'volume'));
% get approximation
% dataApprox = approximateDec(heMatrixScaled, decExt, ...
%     'targetDim', 2, ...
%     'method', 'PCA');
dataApprox = approximateDec(heMatrixScaled, decExt, ...
    'targetDim', 2, ...
    'method', 'nonlinear', ...
    'measureList', 'volume', ...
    'errorType', 'relativeSquared')
heMatrixProjected = dataApprox.projectionMatrix*heMatrixOriginal;
% project decomposition
decApprox = projectDec(decExt, dataApprox.mappingReducedToFull);
%% approximation (dynamics)
% build approx pbeDef
pbeDefApprox.dec = decApprox;
pbeDefApprox.seeding.positionMatrix = dataApprox.mappingFullToReduced*...
    heMatrixProjected(:,1);
pbeDefApprox.seeding.numberVector = 1;
pbeDefApprox.timeSpan = pbeDef.timeSpan;
% pbeDefProj.pivotProperty = @(H, up)(up);
% growth rate: MIMIC PERFECT
% limitingRay = dataApprox.mappingFullToReduced*decExt.se(3,:)'
% limitingNull = null(limitingRay')
% partitionFromhO = @(hO)((limitingNull'*hO < 0) + 1);
% testOne = partitionFromhO([1 1; 2 0.5])
% pbeDefApprox.growthRate = @(H, P, y, t)(dataApprox.mappingFullToReduced*...
%     mapToExtendedEmbedding(pbeDef.growthRate(H(1:2, :), P,y,t), ...
%       decExt, [], partitionFromhO(H)));
% growth rate:  FOURTH APPROACH (mean projection)
nPU = dec.nPU;
projection = cell(1,nPU);
partitionVector = computePartitionIndexDec(dec, dataOriginal.hcMatrix, [], 1)
mappingOtoE = zeros(decExt.nC, dec.nC);
for iPart = 1:nPU
    simplexIndex = dec.unifiedToSimplexVector(iPart);
    
    lamMapping = zeros(dec.nS, dec.nC);
    lamMapping(dec.simplexPartition{simplexIndex}.seIndexVector, :) = ...
        dec.simplexPartition{simplexIndex}.mappingHcToLam;
    
%     projection{iPart} = decExt.se' * lamMapping;
    mappingOtoE = mappingOtoE + ...
        sum(partitionVector==iPart)/length(partitionVector)*...
        decExt.se' * lamMapping;
end
% mappingMatrix = decProper.mappingProperToNew * ...
%     (0.5*projection{2} + 0.5*projection{3});
mappingMatrix = decApprox.properData.mappingProperToNew * mappingOtoE;
pbeDefApprox.growthRate = @(H, P, y, t)(mappingMatrix * ...
    pbeDef.growthRate(H(1:2, 1), P,y,t));

% solve
resultApprox = solvePbe(pbeDefApprox, 'relTol', 1e-5, 'distanceLimit', 1e-6, ...
    'updateUnifiedLimit', 2);
% read ot to single matrix
dataApproxSim = readResultData(resultApprox, pbeDefApprox);
heMatrixApprox = dataApprox.mappingReducedToFull*dataApproxSim.hcMatrix;

%% plot 3D
% close all;
figure, hold on
title('Extended Embedding')
% view extended embedding cone (unified parition)
viewPolytope(decExt.unifiedPartition{1}, 'infinityBoxSize', 10)
% view points of original simulation
plot3(heMatrixOriginal(1,:), heMatrixOriginal(2,:), heMatrixOriginal(3,:), ...
    'k.');
% projected originals
plot3(heMatrixProjected(1,:), heMatrixProjected(2,:), heMatrixProjected(3,:), ...
    'g-');
%! Debug
% hcMatrixProjected = dataApprox.mappingFullToReduced*heMatrixProjected;
% filter = partitionFromhO(hcMatrixProjected) == 1;
% plot3(heMatrixProjected(1,filter), ...
%       heMatrixProjected(2,filter), ...
%       heMatrixProjected(3,filter), 'm.');
% filter = partitionFromhO(hcMatrixProjected) == 2;
% plot3(heMatrixProjected(1,filter), ...
%       heMatrixProjected(2,filter), ...
%       heMatrixProjected(3,filter), 'c.');
% limitinRayMatrix = dataApprox.projectionMatrix*decExt.se(3,:)'*3;
% plot3(limitinRayMatrix(1, :), ...
%       limitinRayMatrix(2, :), ...
%       limitinRayMatrix(3, :), 'b.');
%/
% view points of dynamic evolution
plot3(heMatrixApprox(1,:), heMatrixApprox(2,:), heMatrixApprox(3,:), ...
    'r-');

% 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; 1 1 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(result)
%     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-3);
% plot3(heMatrix(1,:), heMatrix(2,:), heMatrix(3, :), 'r.');
% plot3(heProjMatrix1(1,:), heProjMatrix1(2,:), heProjMatrix1(3, :), 'k.');
% plot3(heProjMatrix2(1,:), heProjMatrix2(2,:), heProjMatrix2(3, :), 'kx');
% view trajectory for matching hC-vectors
% heMatrixProjectedHc = decProj.properData.mappingNewToProper*...
%     projectVectorDec(decProj, heMatrix);
% heMatrixProjectedHc = mapToProperEmbedding(...
%     projectVectorDec(decProj, heMatrix), decProj, 1e-3);
% 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)

%% plot time evolution (projection)
supportMatrixOriginal = computeMeasureDec(decExt, 'hC', heMatrixOriginal, ...
    {'support', extendedSupportDirections});
supportMatrixProjected = computeMeasureDec(decExt, 'hC', heMatrixProjected, ...
    {'support', extendedSupportDirections});
printTimeEvolution(dataOriginal.timeVector, supportMatrixProjected, ...
    supportMatrixOriginal, {'1', '2', '3'})
title('projection')
% ylim([-0.1 0.1]);
measMatrixOriginal = computeMeasureDec(decExt, 'hC', heMatrixOriginal, ...
    {'volume', 'surface', 'meanWidth'});
measMatrixProjected = computeMeasureDec(decExt, 'hC', heMatrixProjected, ...
    {'volume', 'surface', 'meanWidth'});
printTimeEvolution(dataOriginal.timeVector, measMatrixProjected, ...
    measMatrixOriginal, {'volume', 'surface', 'meanWidth'})
title('projection')
% ylim([-0.1 0.1]);

%% plot time evolution (dynamics)
supportMatrixApprox = computeMeasureDec(decExt, 'hC', heMatrixApprox, ...
    {'support', extendedSupportDirections});
printTimeEvolution(dataOriginal.timeVector, supportMatrixApprox, ...
    supportMatrixOriginal, {'1', '2', '3'})
title('dynamic')
% ylim([-0.1 0.1]);
measMatrixApprox = computeMeasureDec(decExt, 'hC', heMatrixApprox, ...
    {'volume', 'surface', 'meanWidth'});
printTimeEvolution(dataOriginal.timeVector, measMatrixApprox, ...
    measMatrixOriginal, {'volume', 'surface', 'meanWidth'})
title('dynamic')
% ylim([-0.1 0.1]);


%% plot some polytopes
% 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))