%% Snyder Cuboid
% ...
clc

%% Generate dynamic analysis data
% dynamicCaseName = 'snyderCuboidDissolution';
dynamicCaseName = 'snyderCuboidGrowth';

% run simulation
% executeLocalPbeCase(dynamicCaseName);

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

% % check covariance at t=0
% disp('### Checking covariance at t=0')
% tmpHcMatrix = result(1).hcMatrix;
% initialCovarianceMatrix = cov(tmpHcMatrix')
% sqrt(abs(diag(initialCovarianceMatrix)))'
% [U,S,V] = svd(tmpHcMatrix', 0);num2str(diag(S), 10)
% disp('#/')

% generate point sets for approximation
generatePointSets(result, pbeDef, ...
    struct('nBin', 50, 'nSample', 1000, ...
           'viewStatistics', 0, 'caseName', dynamicCaseName)...
    )

%% PLOT: variance of PCA-coordinates
% caseList = {'original', 'reduced', 'reducedScaled', 'scaled', 'scaledReduced'};
% plotStyle = {'k-', 'b.-', 'b.--', 'r--', 'r.--'};
% figure,hold on
% yMax = 0;
% for iCase = 1:length(caseList)
%     clear data
%     load([cd filesep 'data' filesep 'snyderCuboid_pointSet_' ...
%         caseList{iCase} '.matx'], '-mat');
%     % svd for pca variances
%     [u,s,v] = svd(data.hcMatrix * data.hcMatrix' * (1/size(data.hcMatrix, 2)), 0);
% %     s = s * (1/s(1,1));
%     yMax = max(yMax, s(2,2));
%     % plot variances
%     plot(diag(s), plotStyle{iCase}, 'markerSize', 15);
% end
% ylim([0 1.5*yMax])
% legend(caseList{:});

%% Approximation
caseList = {'original', 'reduced', 'reducedScaled', 'scaled', 'scaledReduced'};
caseListCore = {'scaledReduced'};
% make projection(s)
generateProjections(pbeDef, 1, dynamicCaseName, caseList);
% generateProjections(pbeDef, 3, dynamicCaseName, caseListCore, ...
%     {'volume', 'volume'});
generateProjections(pbeDef, 4, dynamicCaseName, caseListCore, ...
    {'volume', 'volume'});

%% Evaluate approximation (projection)
% clc
% dataApprox = evaluateConfinedProjection(dynamicCaseName, ...
%     [dynamicCaseName '_1_PCA_scaledReduced']);
dataApprox = evaluateConfinedProjection(dynamicCaseName, ...
    [dynamicCaseName '_4_volume_scaledReduced'])

%% Compute dynamical evolution for approximation
% write pbe definition(s)
% writePbeDefinitionForProjection(dynamicCaseName, 'PCA', caseList, 3);
% writePbeDefinitionForProjection(dynamicCaseName, 'PCA', caseList, 2);
writePbeDefinitionForProjection(dynamicCaseName, 'volume', caseListCore, 4);
%% execute a case
thisCase = [dynamicCaseName '_4_volume_scaledReduced'];
executeLocalPbeCase(thisCase)
[resultTest pbeDefTest] = loadLocalPbeResult(thisCase);

%% number
[X, Y] = extractDataFromPbeResult('integralProperty', 'number', ...
    result, pbeDef);
[X, Yapprox] = extractDataFromPbeResult('integralProperty', 'number', ...
    resultTest, pbeDefTest);
figure,hold on
plot(X, Y(1,:));
plot(X, Yapprox(1,:), 'r.-');
% plot(X, Y(1,:) + Y(2,:), 'r.');
ylabel('number')

%% meanWidth
[X, Y] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
    result, pbeDef);
[X, Yapprox] = extractDataFromPbeResult('integralProperty', 'meanWidth', ...
    resultTest, pbeDefTest);
figure,hold on
plot(X, Y(1,:));
plot(X, Yapprox(1,:), 'r.-');
% plot(X, Y(1,:) + Y(2,:), 'r.');
ylabel('meanWidth')