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

% write pbeDef
pbeDef.dec = dec;
% pbeDef.seeding.positionMatrix = mapToValidityCone([1 1.4 10]', dec, 0) - [0 0 3e-10]';
pbeDef.seeding.positionMatrix = [0.1 0.72 1]';
pbeDef.seeding.numberVector = 1;
pbeDef.growthRate = @(H, P, y, t)(bsxfun(@plus, 0*H, ...
    [0.1 0.58 0.9]'*(t<=1) + ...
    [0.1 0.65 1.1]'*(t>1)));
pbeDef.timeSpan = 0:0.01:2;
% solve
result = solvePbe(pbeDef, 'relTol', 1e-5, 'distanceLimit', 1e-3);

%% project
% % mapNewToProper = [1 0; 0 1; 0.6 0.6];   
mapNewToProper = [   -0.6236   -0.0000; ...
   -0.7804    0.0573; ...
    0.0448    0.9984];

% mapNewToProper = [dec.se(2,:)', dec.se(4,:)']
decProj = projectDec(dec, mapNewToProper);
% write pbeDef
pbeDefProj.dec = decProj;
pbeDefProj.seeding.positionMatrix = mapFromProperEmbedding(...
    projectVectorDec(pbeDef.seeding.positionMatrix, decProj), decProj);
pbeDefProj.seeding.numberVector = 1;
pbeDefProj.pivotProperty = @(H, up)(up);
pbeDefProj.growthRate = @(H, P, y, t)(...
    pbeDef.growthRate([H; H(1,:)],P,y,t));
pbeDefProj.timeSpan = pbeDef.timeSpan;
% pbeDefProj.timeSpan = 0:0.01:0.58;
% solve
resultProj = solvePbe(pbeDefProj, 'relTol', 1e-5, 'distanceLimit', 1e-6, ...
    'updateUnifiedLimit', 2);

%% 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; 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)

%% more
% figure,viewCrystal(dec.cdef, result(1).hcMatrix)
% figure,viewCrystal(dec.cdef, result(101).hcMatrix)
% figure,viewCrystal(dec.cdef, result(end).hcMatrix)

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