%% Snyder Cuboid
% In this case, we construct the approximation by hand, leading to model
%   reduction.
% We also compare it to a "simple" approach that will cease working when
%   growth / dissolution takes place.

clc
% clear all

dynamicCaseName = 'snyderCuboidGrowthDissolution';

%% preconditioning
% execute simulation
executeLocalPbeCase(dynamicCaseName);
% load simulation data
[result pbeDef] = loadLocalPbeResult(dynamicCaseName);
% generate point sets for approximation
generatePointSets(result, pbeDef, ...
    struct('nBin', 50, 'nSample', 1000, ...
           'viewStatistics', 0, 'caseName', dynamicCaseName)...
    )

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

shapeMatrix = [eye(3) 100*ones(3, 4)]';
% shapeMatrix = [1 2 3 100 100 100 100]';
shapeMatrix = mapToValidityCone(shapeMatrix, pbeDef.dec, 0);
% GROWTH
spanMatrix = [1./pbeDef.dec.cdef.relativeGrowthCaseTwo ...
    shapeMatrix];
% DISSOLUTION
% spanMatrix = [pbeDef.dec.cdef.relativeGrowthCaseTwo ...
%     shapeMatrix];
% GROWTH + DISSOLUTION
% spanMatrix = [shapeMatrix ...
%     1./pbeDef.dec.cdef.relativeGrowthCaseTwo ...
%     pbeDef.dec.cdef.relativeGrowthCaseTwo
%     ];

% orthogonalize
spanMatrix = orth(spanMatrix);

% check shapes
% figure,viewPolytope(struct('A', pbeDef.dec.A, 'h', ...
%     pbeDef.dec.mappingReducedToFull*spanMatrix(:,1)));
% figure,viewPolytope(struct('A', pbeDef.dec.A, 'h', ...
%     pbeDef.dec.mappingReducedToFull*spanMatrix(:,2)));
% figure,viewPolytope(struct('A', pbeDef.dec.A, 'h', ...
%     pbeDef.dec.mappingReducedToFull*spanMatrix(:,3)));
% figure,viewPolytope(struct('A', pbeDef.dec.A, 'h', ...
%     pbeDef.dec.mappingReducedToFull*spanMatrix(:,4)));

thisFile = ['./database/' dynamicCaseName '_4_handmade_.matx'];
if ~exist(thisFile, 'file')
    dec = projectDec(pbeDef.dec, spanMatrix, 'zerotol', 1e-8);
    dataApprox.dec = dec;
    dataApprox.cdef = pbeDef.dec.cdef;
    save(thisFile, '-mat', 'dataApprox');
% else
%     data = load(thisFile, '-mat');
end
    
caseList = {''};

%% Evaluate approximation (projection)
% dataApprox = evaluateConfinedProjection(dynamicCaseName, ...
%     [dynamicCaseName '_4_handmade_'])

%% Evaluate "simple" approach
resultApprox = evaluateSimpleApproach(pbeDef);
evaluateConfinedSimulation({result, pbeDef}, {resultApprox, pbeDef}, 0)

%% Compute dynamical evolution for approximation
% writePbeDefinitionForProjection(dynamicCaseName, 'handmade', caseList, 4);
thisCase = [dynamicCaseName '_4_handmade_'];
[resultApprox pbeDefApprox] = executeLocalPbeCase(thisCase);

evaluateConfinedSimulation({result, pbeDef}, {resultApprox, pbeDefApprox}, 1)
