%% Snyder Cuboid
% ...
% TODO % TEst1
clc
clear all

caseNameBase = 'ureaGenericGrowth';
caseList = {'original', 'reduced', 'reducedScaled', 'scaled', 'scaledReduced'};
caseListNonlinear = {'scaled'};
measureList = {'volume', 'surface', 'meanWidth'};
% TODO % Test 2

%% NON-AP: executePbe
executeLocalPbeCase(caseNameBase);

%% NON-AP: load data
% dec = obtainCrystalDec('snyder_cuboid');
[result pbeDef] = loadLocalPbeResult(caseNameBase);


%% PRE-AP: generate point sets for approximation
generatePointSets(result, pbeDef, ...
    struct('nBin', 20, 'nSample', 50, 'nBatch', 1, ...
           'viewStatistics', 0, 'caseName', caseNameBase)...
    )
%% PLOT: variance of PCA-coordinates
% plotStyle = {'k-', 'b.-', 'b.--', 'r--', 'r.--'};
% figure,hold on
% yMax = 0;
% for iCase = 1:length(caseList)
%     clear data
%     load([cd filesep 'data' filesep caseNameBase '_pointSet_' ...
%         caseList{iCase} '.matx'], '-mat');
%     % svd for pca variances
%     [u,s,v] = svd(dataPointSet.hcMatrix * dataPointSet.hcMatrix' * ...
%         (1/size(dataPointSet.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{:});
% 
%% AP-PCA: generate projections
% make projection(s)
% generateProjections(pbeDef, 2, caseNameBase, caseList, {'volume', ...
%     'volume'});
generateProjections(pbeDef, 2, caseNameBase, caseList, {'volume', ...
    {'volume', 'support', pbeDef.dec.A(1,:)}});
% write pbe definition(s)
writePbeDefinitionForProjection(caseNameBase, 'PCA', caseList, 2);
writePbeDefinitionForProjection(caseNameBase, 'volume', caseList, 2);

% thisCase = [caseNameBase '_2_PCA_scaledReduced'];
% thisCase = [caseNameBase '_2_PCA_original'];
% thisCase = [caseNameBase '_2_PCA_scaled'];
thisCase = [caseNameBase '_2_volume_scaled'];
%% evaluate a case
executeLocalPbeCase(thisCase)
[resultApprox pbeDefApprox] = loadLocalPbeResult(thisCase);
dataReference = readResultData(result, pbeDef);
dataApprox = readResultData(resultApprox, pbeDefApprox);

%% compare to direct projection
heMatrixReference = dataReference.hcMatrix;
heUnifiedVector = computePartitionIndexDec(pbeDef.dec, heMatrixReference, 1e-3, 1);
hcMatrix = mapFromProperEmbedding(projectVectorDec(...
    heMatrixReference, pbeDefApprox.dec, 1e-4), pbeDefApprox.dec);
hcUnifiedVector = computePartitionIndexDec(pbeDefApprox.dec, hcMatrix, 1e-10, 1);
heMatrix = mapToProperEmbedding(hcMatrix, pbeDefApprox.dec);
% hcMatrix = mapFromProperEmbedding(heMatrixReference, pbeDefApprox.dec, 1e-3);
% heMatrix = pbeDefApprox.dec.properData.mappingNewToProper * hcMatrix;
% coordinate space
printTimeEvolution(dataReference.timeVector, [heMatrix; hcUnifiedVector], ...
    [heMatrixReference; heUnifiedVector], {'hc1', 'hc2', 'hc3', 'iPU'})
ylim([-0.1 0.1]);
% printTimeEvolution(dataReference.timeVector, heMatrix, ...
%     heMatrixReference, {'hc1', 'hc2', 'hc3'})
% measures
measMatrixReference = computeMeasureDec(pbeDef.dec, 'hC', ...
    heMatrixReference, {'volume', 'surface', 'meanWidth'}, 1e-4);
measMatrix = computeMeasureDec(pbeDef.dec, 'hC', ...
    heMatrix, {'volume', 'surface', 'meanWidth'}, 1e-4, [], 1);
printTimeEvolution(dataReference.timeVector, ...
    [measMatrix; isValidDec(...
       pbeDefApprox.dec.properData.projectionMatrix*heMatrixReference, ...
       pbeDefApprox.dec.properData)], ...
    [measMatrixReference; measMatrix(1,:)*0+1], ...
    {'volume', 'surface', 'meanWidth', 'valid'})
ylim([-0.1 0.1]);

%% get some number results
errorStruct = computeMatrixError(measMatrix, measMatrixReference, 'all')
errorStruct.relativeSquared'
errorStruct.relativeMean'

%% compare to time evolution
hcMatrix = dataApprox.hcMatrix;
% hcUnifiedVector = computePartitionIndexDec(pbeDefApprox.dec, hcMatrix, 1e-3, 1);
heMatrix = mapToProperEmbedding(hcMatrix, pbeDefApprox.dec, 1e-3);
% coordinate space
% printTimeEvolution(dataReference.timeVector, [heMatrix; hcUnifiedVector], ...
%     [heMatrixReference; heUnifiedVector], {'hc1', 'hc2', 'hc3', 'iPU'})
printTimeEvolution(dataReference.timeVector, ...
    [heMatrix; isValidDec(hcMatrix, pbeDefApprox.dec)], ...
    [heMatrixReference; heMatrixReference(1,:)*0], ...
    {'hc1', 'hc2', 'hc3', 'valid'})
ylim([-0.1 0.1]);
title('time evolution');
% measures
measMatrix = computeMeasureDec(pbeDef.dec, 'hC', ...
    heMatrix, {'volume', 'surface', 'meanWidth'}, 1e-4, [], 1);
printTimeEvolution(dataReference.timeVector, measMatrix, ...
    measMatrixReference, {'volume', 'surface', 'meanWidth'})
ylim([-0.1 0.1]);
title('time evolution');

%% get some number results
errorStruct = computeMatrixError(measMatrix, measMatrixReference, 'all')
errorStruct.relativeSquared'
errorStruct.relativeMean'