%% extra script for unifying approximation
clc
clear all
% close all
dynamicCaseName = 'ureaGenericGrowth';

%% prepare approx data set
% run simulation (should already be available by runParacetamolGrowth)
executeLocalPbeCase(dynamicCaseName);

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

% prepare point set
% [hcMatrix groupingVector] = resampleData(pbeDef, result, 10*10^-6, 100);
hcMatrix = [result(:).hcMatrix];

% obtain required embedding
thisFileName = ['database' filesep dynamicCaseName '_decExt.matx'];
if exist(thisFileName, 'file')
    data = load(thisFileName, '-mat');
    decExt = data.decExt;
else
    decExt = computeUnifyingDec(pbeDef.dec);
%     decExt.measureData = rmfield(decExt.measureData, 'support');
    decExt = addMeasureDataDec(decExt, 'support', decExt.A);
    save(thisFileName, 'decExt');
end

%% perform approximation
% PCA is sufficient going from 4D to 9D back to 4D since the mean
%   projection only spans 4D in the 9D embedding, anyway
rng(1)
generateSingleProjection('ureaGenericGrowth_3u_PCA', pbeDef, hcMatrix, ...
    {'targetDim', 3, ...
     'method', 'PCA'}, 0, decExt)
extraName = '_3u_PCA';
dataApprox = load(['./database/' dynamicCaseName extraName '.matx'], '-mat');

%% simulate monomials
momentData = simulateUreaMoments(pbeDef, dataApprox);

%% Print some figures
% close al

% [resultTest pbeDefTest] = loadLocalPbeResult([dynamicCaseName extraName]);

% prepare projection
% heMatrix = mapToExtendedEmbedding(result.hcMatrix, decExt, []);
% heApproxMatrix = dataApprox.meanMappingMatrix*result.hcMatrix;

%% 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 = dataApprox.meanMappingMatrix * ...
        result(iTime).hcMatrix;
    resultProj(iTime).hcMatrix = ...
        dataApprox.dataApprox.dec.properData.mappingProperToNew*...
        resultProj(iTime).hcMatrix;
    resultProj(iTime).upVector = resultProj(iTime).upVector*0 + 1;
end
pbeDefProj = pbeDef;
% pbeDefProj.dec = decExt;
pbeDefProj.dec = dataApprox.dataApprox.dec;
[X, Yproj] = extractDataFromPbeResult('integralProperty', 'volume', ...
    resultProj, pbeDefProj);

% 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, momentData.volumeVector, 'r.-');
plot(X, Yproj(1,:), 'm.-');
% plot(Yproj(1,:), 'r.');
ylabel('volume')
figure, hold on
plot(X, (momentData.volumeVector-Y(1,:))./Y(1,:), 'b-')
plot(X, (Yproj(1,:)-Y(1,:))./Y(1,:), 'm-')