%% Paracetamol Growth
clc
dynamicCaseName = 'paracetamolGrowth';

%% simulation
% execute simulation
% executeLocalPbeCase(dynamicCaseName);
% load simulation data
[result pbeDef] = loadLocalPbeResult(dynamicCaseName);

%% graphs
% close all
% propertyList = {'total crystal volume (m^3)', ...
%     'temperature', 'concentration', 'saturation'};
% % bulk properties
% for iProp = 1:4
%     [timeVector totalVolumeVector] = extractDataFromPbeResult(...
%     'bulkProperty', iProp, result, pbeDef);
%     figure, plot(timeVector, totalVolumeVector);
%     title(propertyList{iProp})
% end
% 
% % particle trajectory
% iParticle = 1;
% nTime = length(result)
% sizeMatrix = zeros(4,nTime)
% distanceVector = zeros(1,nTime);
% for iTime = 1:nTime
%     sizeMatrix(:,iTime) = result(iTime).hcMatrix(:,iParticle);
%     distanceVector(iTime) = max(max(...
%         pbeDef.dec.validityCone.A * result(iTime).hcMatrix));
% end
% figure,hold on,title('facet distances')
% plot(timeVector, sizeMatrix(1,:), 'k-')
% plot(timeVector, sizeMatrix(2,:), 'r-')
% plot(timeVector, sizeMatrix(3,:), 'b-')
% plot(timeVector, sizeMatrix(4,:), 'm-')
% 
% figure,hold on,title('relative facet distances')
% plot(timeVector, sizeMatrix(1,:)./sizeMatrix(1,:), 'k-')
% plot(timeVector, sizeMatrix(2,:)./sizeMatrix(1,:), 'r-')
% plot(timeVector, sizeMatrix(3,:)./sizeMatrix(1,:), 'b-')
% plot(timeVector, sizeMatrix(4,:)./sizeMatrix(1,:), 'm-')
% 
% % distance from validity cone
% figure
% plot(timeVector, distanceVector, 'k-')
% title('distance from validity cone')

%% Generate approximations

% evenly resample approach
rng(1)
generateSingleProjection('paracetamolGrowth_2_PCA_r1000', pbeDef, result, ...
    {'targetDim', 2, ...
     'nResample', 1000, ...
     'method', 'PCA'}, 0)
rng(1)
generateSingleProjection('paracetamolGrowth_2_volume_relmax_r1000', pbeDef, result, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'nResample', 1000, ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'numStageOnePoints', 50 ...
     'numtrialPoints', 500, ...
     'errorType', 'relativeMax', ...
     'measureWeightVector', []}, 0)
 
% time based resample approach
[heMatrix groupingVector] = resampleData(pbeDef, result, 10*10^-6, 100);
rng(1)
generateSingleProjection('paracetamolGrowth_2_PCA_t100_10', pbeDef, heMatrix, ...
    {'targetDim', 2, ...
     'method', 'PCA'}, 0)
rng(1)
generateSingleProjection('paracetamolGrowth_2_volume_relmax_t100_10', pbeDef, heMatrix, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'numStageOnePoints', 50 ...
     'numtrialPoints', 500, ...
     'errorType', 'relativeMax', ...
     'errorGroupVector', groupingVector, ...
     'measureWeightVector', []}, 0)
rng(1)
generateSingleProjection('paracetamolGrowth_2_volume_int_t100_10', pbeDef, heMatrix, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'numStageOnePoints', 50 ...
     'numtrialPoints', 500, ...
     'errorType', 'relativeMaxGroupedMean', ...
     'errorGroupVector', groupingVector, ...
     'measureWeightVector', []}, 0)
 
% 1D approximation
rng(1)
generateSingleProjection('paracetamolGrowth_1_PCA_t100_10', pbeDef, heMatrix, ...
    {'targetDim', 1, ...
     'method', 'PCA'}, 0)
rng(1)
generateSingleProjection('paracetamolGrowth_1_volume_int_t100_10', pbeDef, heMatrix, ...
    {'targetDim', 1, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'numStageOnePoints', 50 ...
     'numtrialPoints', 500, ...
     'errorType', 'relativeMaxGroupedMean', ...
     'errorGroupVector', groupingVector, ...
     'measureWeightVector', []}, 1)


%% Run approximation simulation
% executeLocalPbeCase([dynamicCaseName '_2_PCA_r'])
% executeLocalPbeCase([dynamicCaseName '_2_PCA_t'])
% executeLocalPbeCase([dynamicCaseName '_1_PCA_t'])
% executeLocalPbeCase([dynamicCaseName '_2_volRelMax_r'])
% executeLocalPbeCase([dynamicCaseName '_2_volRelMax_t'])
% executeLocalPbeCase([dynamicCaseName '_2_volInt_t'])
executeLocalPbeCase([dynamicCaseName '_1_volInt_t'])


%% Print some figures
% close all
% extraName = '_2_PCA_r';
% extraName = '_2_PCA_t';
% extraName = '_1_PCA_t';
% extraName = '_2_volRelMax_r';
extraName = '_2_volRelMax_t';
% extraName = '_2_volInt_t';
% extraName = '_1_volInt_t';

[resultTest pbeDefTest] = loadLocalPbeResult([dynamicCaseName extraName]);
% dataApprox = load(['./database/' dynamicCaseName extraName '.matx'], '-mat');
% heMatrix = dataApprox.dataApprox.heMatrix;
% heMatrixApprox = dataApprox.dataApprox.projectionMatrix * heMatrix;

%% meanWidth over time
% [X, Y] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
%     result, pbeDef);
% [X, Yapprox] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
%     resultTest, pbeDefTest);
% % Y = computeMeasureDec(pbeDef.dec, 'hC', heMatrix, 'meanWidth', ...
% %     [], [], 1);
% % Yproj = computeMeasureDec(pbeDef.dec, 'hC', heMatrixApprox, 'meanWidth', ...
% %     [], [], 1);
% figure,hold on
% plot(Y(1,:), 'b-x');
% % plot(Yapprox(1,:), 'r.-');
% plot(Yproj(1,:), 'g-');
% ylabel('meanWidth')
% figure
% % plot(X, (Y(1,:)-Yapprox(1,:))./Y(1,:), 'b-')
% % hold on
% plot((Y(1,:)-Yproj(1,:))./Y(1,:), 'm-')

%% volume over time
[X, Y] = extractDataFromPbeResult('integralProperty', 'volume', ...
    result, pbeDef);
[X, Yapprox] = extractDataFromPbeResult('integralProperty', 'volume', ...
    resultTest, pbeDefTest);
resultProj = result;
for iTime = 1:length(result)
    resultProj(iTime).hcMatrix = pbeDefTest.dec.properData.projectionMatrix * ...
        result(iTime).hcMatrix;
end
[X, Yproj] = extractDataFromPbeResult('integralProperty', 'volume', ...
    resultProj, pbeDef);

% Y = computeMeasureDec(pbeDef.dec, 'hC', heMatrix, 'volume', ...
%     [], [], 1);
% Yproj = computeMeasureDec(pbeDef.dec, 'hC', heMatrixApprox, 'volume', ...
%     [], [], 1);
figure,hold on
plot(X, Y(1,:), 'b-x');
plot(X, Yapprox(1,:), 'r.-');
plot(X, Yproj(1,:), 'm.-');
% plot(Yproj(1,:), 'r.');
ylabel('volume')
figure, hold on
plot(X, (Yapprox(1,:)-Y(1,:))./Y(1,:))
plot(X, (Yproj(1,:)-Y(1,:))./Y(1,:), 'm-')

%% volume distribution
% minSize = 100*10^-6;
% maxSize = 300*10^-6;
% deltaSize = 2*10^-6;
% minSize = 0;
% maxSize = 90*10^-13;
% deltaSize = 0.01*(maxSize-minSize);
% sizeVector = minSize:deltaSize:maxSize;
% sizeMatrix = [sizeVector-0.5*deltaSize; sizeVector+0.5*deltaSize];
% [X, Y] = extractDataFromPbeResult('propertyDistribution', ...
%     {'volume', 'volume', sizeMatrix}, result, pbeDef);
% [X, Yapprox] = extractDataFromPbeResult('propertyDistribution', ...
%     {'volume', 'volume', sizeMatrix}, resultTest, pbeDefTest);
% 
% %% show off
% figure,hold on
% nSample = length(sizeVector);
% for iTime = [1,25,35,45,60,80,350]
%     plot(sizeVector, Y(1:nSample, iTime), 'k.-')
%     plot(sizeVector, Yapprox(1:nSample, iTime), 'r-')
% end