clear all;
caseName = 'growthDensitySampleMetroMix';

%% load result
[result pbeDefinition] = loadCase(caseName);
close all;

%% parameters
m = 1.125e-6;
off = 0.7*m;
d = m/100;

% %% initial distribution
figure(1),hold on
set(1, 'renderer', 'zbuffer')
gridSpacing = (m-off):d:(m+off);

% add validity cone
hrep.A = [pbeDefinition.dec.validityCone.A; ...
          1*eye(2)];
hrep.h = [pbeDefinition.dec.validityCone.h;...
          m+off; m+off];
hrep.A = [hrep.A zeros(length(hrep.h), 1); ...
          0 0 1; 0 0 -1];
hrep.h = [hrep.h;0; 0];
pHandle = viewPolytope(hrep, 'color', [0.5 1 0.5], 'shade', 0.2);
for iPatch = 1:length(pHandle)
    set(pHandle(iPatch), 'zData', ...
        get(pHandle(iPatch), 'zData') + 1e21) 
end

% sampled distribution
[X, Y, Z] = computePatchDataFromPointCloud(...
        [result(1).hcMatrix; ...
         result(1).numberVector .* result(1).probabilityVector]);
patch(X, Y, Z, Z, 'lineWidth', 0.5)

% add sampling points
plot3(result(1).hcMatrix(1,:), ...
      result(1).hcMatrix(2,:), ...
      result(1).numberVector .* result(1).probabilityVector, ...
      'rd', 'markerSize', 5, 'markerEdgeColor', [1 0.9 0.9], ...
      'markerFaceColor', [1 0 0], 'lineWidth', 1);

set(gca, 'PlotBoxAspectRatioMode', 'auto', ...
      'DataAspectRatioMode', 'auto')
view([10 25])
xlim([m-off, m+off]);ylim([m-off, m+off])
zlim([0 2e23])
 
% print
print -dpng -r300 figure/metroMix/initialDistributionSampled.png
close all

%% sampled distributions
timeVector = [1 150 250 396];
disp('starting')
off = 2e-6;
for iTime = timeVector
    thisTime = result(iTime).time
    thisSimFraction = iTime / length(result)
    
    [X, Y, Z] = computePatchDataFromPointCloud(...
        [result(iTime).hcMatrix; ...
        result(iTime).numberVector .* result(iTime).probabilityVector]);
    mHere = mean(result(iTime).hcMatrix, 2);
    
    hFigure = figure, set(hFigure, 'renderer', 'zbuffer');
    patch(X, Y, Z, Z, 'lineWidth', 0.5)
    xlim([mHere(1)-off, mHere(1)+off]);
    ylim([mHere(2)-off, mHere(2)+off]);
    zlim([0 3e23]);set(gca, 'clim', [0 3e23]);
    view([0 90])
    axis off
    box off
    
    print('-dpng', '-r300', ['figure/metroMix/density_' num2str(iTime) '.png']);
    close all
end

%% plot volume over time
[X, Y] = extractDataFromPbeResult('integralProperty', 'volume', result, pbeDefinition);
%%
figure,hold on
set(1, 'paperPosition', [0 0 30 15])

plot(X, Y(1,:), 'k-', 'lineWidth', 3)
plot(X, Y(1,:)-Y(2,:), 'r.'), plot(X, Y(1,:)+Y(2,:), 'r.')
set(gca,'FontSize',20)
box off
legend({'total volume (left axis)', 'lower error (std)', 'upper error (std)'},...
    'Location', 'NorthWest');
% xlabel('time in s')
ylabel('total volume in m^3')
% error extra y-axis
axisHandle = axes;
plot(X, Y(2,:)./Y(1,:), 'k--', 'parent', axisHandle, 'lineWidth', 3)
ylim([0 0.025])
set(gca,'FontSize',20)
box off
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])
legend('coef. of variation (right axis)', 'Location', 'NorthEast')
ylabel('coefficient of variation')

% print
print -dpng -r300 figure/metroMix/volumeOverTime.png
close all
