%% Snyder Cuboid Dissolution
% This one runs all simulations and stores the results (and shows some
% figure).
clc
dynamicCaseName = 'ureaGenericGrowth';

%% Run original simulation(s)
% run simulation
executeLocalPbeCase(dynamicCaseName);
% run simulation without limits
executeLocalPbeCase([dynamicCaseName '_noGrowthLimit']);

% load data (no limits)
[result pbeDef] = loadLocalPbeResult([dynamicCaseName '_noGrowthLimit']);
[result2 pbeDef2] = loadLocalPbeResult([dynamicCaseName]);

%% view trajectory
% hTrajectory = figure; hold on
% % viewPolytope(pbeDef.dec.validityCone, 'infinityBoxSize', 2, ...
% %     'shade', 0.2)
% heMatrix = nan(3, length(result));
% for iTime = 1:length(result)
%     heMatrix(:,iTime) = result(iTime).hcMatrix;
% end
% plot3(heMatrix(1,:), heMatrix(3,:), heMatrix(3,:))
% xlabel('1'),ylabel('2'),zlabel('3')
% xlim([0 max(heMatrix(:))])
% ylim([0 max(heMatrix(:))])
% zlim([0 max(heMatrix(:))])

%% Generate approximations
rng(1)
generateSingleProjection('ureaGenericGrowth_2_PCA', pbeDef, result, ...
    {'targetDim', 2, ...
     'method', 'PCA'})
rng(1)
generateSingleProjection('ureaGenericGrowth_2_volume_abssq', pbeDef, result, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'errorType', 'absoluteSquared', ...
     'measureWeightVector', []}, 0)
rng(1)
generateSingleProjection('ureaGenericGrowth_2_volume_relmax', pbeDef, result, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'numStageOnePoints', 1000 ...
     'numtrialPoints', 2000, ...
     'errorType', 'relativemax', ...
     'measureWeightVector', []}, 0)

rng(1)
generateSingleProjection('ureaGenericGrowth_2_volume_relsq', pbeDef, result, ...
    {'targetDim', 2, ...
     'method', 'nonlinear', ...
     'projectionType', 'true', ...
     'measureList', 'volume', ...
     'errorType', 'relativeSquared', ...
     'measureWeightVector', []}, 0)
 
%% ratio-wise approx
rng(1)
ratioVector = 0:0.025:1;
nRatio = length(ratioVector);

meanWidthError = nan(1, nRatio);
volumeError = nan(1, nRatio);
totalError = nan(1, nRatio);
grassMatrix = nan(2, nRatio);
for iRatio = 1:nRatio
    thisRatio = ratioVector(iRatio);
    thisName = [dynamicCaseName '_2_ratio_' ...
        num2str(round(thisRatio*100))];
    
    disp(['Generating ratio: ' num2str(iRatio)])
    generateSingleProjection(thisName, ...
        pbeDef, result, {...
        'targetDim', 2, ...
        'method', 'nonlinear', ...
        'projectionType', 'true', ...
        'measureList', {'volume', 'meanWidth'}, ...
        'errorType', 'relativemax', ...
        'numStageOnePoints', 50, ...
        'numtrialPoints', 5000, ...
        'measureWeightVector', [thisRatio (1-thisRatio)]}, 0)

    dataApprox = load(['./database/' thisName '.matx'], '-mat');
    grassMatrix(:, iRatio) = dataApprox.dataApprox.nonlinearData.globalX;
    
    % recompute measures
    heMatrix = dataApprox.dataApprox.heMatrix;
    heMatrixApprox = dataApprox.dataApprox.heMatrixApprox
    
    thisVolume = computeMeasureDec(pbeDef.dec, 'hC', heMatrix, ...
        'volume', [], [], 1);
    thisVolumeApprox = computeMeasureDec(pbeDef.dec, 'hC', heMatrixApprox, ...
        'volume', [], [], 1);
    
    thisMeanWidth = computeMeasureDec(pbeDef.dec, 'hC', heMatrix, ...
        'meanWidth', [], [], 1);
    thisMeanWidthApprox = computeMeasureDec(pbeDef.dec, 'hC', heMatrixApprox, ...
        'meanWidth', [], [], 1);
    
    % recompute errors
    meanWidthError(iRatio) = max(...
        (thisMeanWidth - thisMeanWidthApprox)./thisMeanWidth);
    volumeError(iRatio) = max(...
        (thisVolume - thisVolumeApprox)./thisVolume);
    totalError(iRatio) = thisRatio*volumeError(iRatio) + ...
        (1-thisRatio)*meanWidthError(iRatio);
    grassMatrix(:, iRatio) = dataApprox.dataApprox.nonlinearData.globalX;
end
% save data
save('./data/ureaApproxRatioSummary.mat', 'ratioVector', ...
    'grassMatrix', 'volumeError', 'totalError', 'meanWidthError')

figure,hold on
plot(ratioVector, grassMatrix(1,:)/pi, 'b.-')
plot(ratioVector, grassMatrix(2,:)/pi, 'r.-')
title('grass')

figure,hold on
plot(ratioVector, totalError, 'k.-')
plot(ratioVector, volumeError, 'r.-')
plot(ratioVector, meanWidthError, 'b.-')
title('measError')
 
%% Run approximation simulation
executeLocalPbeCase([dynamicCaseName '_2_PCA'])
executeLocalPbeCase([dynamicCaseName '_2_volume_abssq'])
executeLocalPbeCase([dynamicCaseName '_2_volume_relmax'])
executeLocalPbeCase([dynamicCaseName '_2_volume_relsq'])

%% Print some figures
% extraName = '_2_PCA';
extraName = '_2_volume_abssq';
% extraName = '_2_volume_relmax';
% extraName = '_2_volume_relsq';
[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);
Yproj = computeMeasureDec(pbeDef.dec, 'hC', heMatrixApprox, 'meanWidth', ...
    [], [], 1);
figure,hold on
plot(X, Y(1,:), 'b-x');
plot(X, Yapprox(1,:), 'r.-');
plot(X, Yproj(1,:), 'g-');
% plot(X, Y(1,:) + Y(2,:), 'r.');
ylabel('meanWidth')
figure
figure,plot(X, (Y(1,:)-Yapprox(1,:))./Y(1,:), 'b-')
hold on,plot(X, (Y(1,:)-Yproj(1,:))./Y(1,:), 'm-')

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