% clear all
%% Dynamic potash alum (growing out facets)
dec = obtainCrystalDec('generic_potash_alum', 3);
dec = addMeasureDataDec(dec, 'volume');
dec = addMeasureDataDec(dec, 'surface');
dec = addMeasureDataDec(dec, 'meanWidth');


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

%% project
% % scaling
% heMatrixScaled = bsxfun(@times, heMatrixOriginal, ...
%     1./ computeMeasureDec(dec, 'hC', heMatrixOriginal, 'volume', [], [], 1));
% % get approximation
% % dataApprox = approximateDec(heMatrixScaled, dec, ...
% %     'targetDim', 2, ...
% %     'method', 'PCA');
% heMatrixScaled = mapToValidityCone(heMatrixScaled, dec, 0);
% dataApprox = approximateDec(heMatrixScaled, dec, ...
%     'targetDim', 2, ...
%     'method', 'nonlinear', ...
%     'measureList', 'volume', ...
%     'errorType', 'relativeSquared');
% dataApprox.fadeOutVector
% construct mapping
dataApprox.mappingReducedToFull = orth(...
    [result(1).hcMatrix [1 3 2]']);
% project decomposition
decApprox = projectDec(dec, dataApprox.mappingReducedToFull);
% project points
heMatrixProjected = projectVectorDec(heMatrixOriginal, decApprox);
hcMatrixProjected = decApprox.properData.mappingProperToNew * heMatrixProjected;
% heMatrixProjected = mapToValidityCone(heMatrixOriginal, dec, 0);
% heMatrixProjected = dataApprox.heMatrixProjected;

%% write pbeDef (approx)
% build approx pbeDef
pbeDefApprox.dec = decApprox;
% pbeDefApprox.seeding.positionMatrix = dataApprox.mappingFullToReduced*...
%     heMatrixProjected(:,1);
pbeDefApprox.seeding.positionMatrix = mapFromProperEmbedding(...
    pbeDef.seeding.positionMatrix, decApprox);
pbeDefApprox.seeding.numberVector = 1;
pbeDefApprox.growthRate = @(H, P, y, t)(pbeDef.growthRate(H([1 1 1], :), P,y,t));
pbeDefApprox.timeSpan = pbeDef.timeSpan;
% solve
resultApprox = solvePbe(pbeDefApprox, 'relTol', 1e-5, 'distanceLimit', 1e-8, ...
    'updateUnifiedLimit', 2);
% read out to single matrix
dataApproxSim = readResultData(resultApprox, pbeDefApprox);
% heMatrixApprox = dataApprox.mappingReducedToFull*dataApproxSim.hcMatrix;
% heMatrixApprox = mapToProperEmbedding(mapToConfinementCone(...
%     dataApproxSim.hcMatrix, decApprox, 0), decApprox);
heMatrixApprox = mapToProperEmbedding(dataApproxSim.hcMatrix, decApprox);
% heMatrixApprox = mapToValidityCone(heMatrixApprox, dec);

%% plot 3D
% close all;
figure, hold on
title('Extended Embedding')
% view extended embedding cone (unified parition)
viewPolytope(dec.validityCone, 'infinityBoxSize', 2.5)
% view points of original simulation
plot3(heMatrixOriginal(1,:), heMatrixOriginal(2,:), heMatrixOriginal(3,:), ...
    'k.');
% projected originals
plot3(heMatrixProjected(1,:), heMatrixProjected(2,:), heMatrixProjected(3,:), ...
    'g-');
% view points of dynamic evolution
plot3(heMatrixApprox(1,:), heMatrixApprox(2,:), heMatrixApprox(3,:), ...
    'rx');

% 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)
printTimeEvolution(dataOriginal.timeVector, heMatrixProjected, ...
    heMatrixOriginal, {'1', '2', '3'})
title('projection')
% ylim([-0.1 0.1]);
measMatrixOriginal = computeMeasureDec(dec, 'hC', heMatrixOriginal, ...
    {'volume', 'surface', 'meanWidth'}, [], [], 1);
measMatrixProjected = computeMeasureDec(dec, 'hC', heMatrixProjected, ...
    {'volume', 'surface', 'meanWidth'}, [], [], 1);
printTimeEvolution(dataOriginal.timeVector, measMatrixProjected, ...
    measMatrixOriginal, {'volume', 'surface', 'meanWidth'})
title('projection')
% ylim([-0.1 0.1]);

%% plot time evolution (dynamics)
printTimeEvolution(dataOriginal.timeVector, heMatrixApprox, ...
    heMatrixOriginal, {'1', '2', '3'})
title('dynamic')
% ylim([-0.1 0.1]);
measMatrixApprox = computeMeasureDec(dec, 'hC', heMatrixApprox, ...
    {'volume', 'surface', 'meanWidth'}, [], [], 1);
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))