clear all
%% Note for confidence of std
% Below, we show the coefficient of variation from the estimates - these
%   are based on 100 points and the c.o.v based on the 10 simulation runs.
%   The confidence interval of the second estimation is extremely large and
%   goes from 69% to 183% of the estimated value while the first estimate
%   ranges from 88% to 116%
% CODE: 
% LEN=20;alpha=0.05;
% lowBound = sqrt((LEN-1)/chi2inv(1-alpha/2,LEN-1))
% highBound =sqrt((LEN-1)/chi2inv(alpha/2,LEN-1))

%% generate or load results
tmp = executeCase('init');
caseList = {'constGrowth_normal', 'constGrowth_metroNum', ...
    'constGrowth_metroVol', 'constGrowth_metroMix'};
nameList = {'Normal', 'MetroNum', 'MetroVol', 'MetroMix'};
% caseList = {'constGrowth_normal', 'constGrowth_metroVol'};
% nameList = {'Normal','MetroVol'};
% caseList = {'constGrowth_normal'};
% nameList = {'Normal'};
for iCase = 1:length(caseList)
for iSeed = 0:9
    thisCaseName = [caseList{iCase} '_' num2str(iSeed)]
    try
        [result pbeDef] = loadCase(thisCaseName);
    catch
        pbeDef = obtainData(thisCaseName);
        result = executeCase(thisCaseName);
    end
    eval(['result' nameList{iCase} '{' num2str(iSeed+1) '}=result;'])
    eval(['pbeDef' nameList{iCase} '=pbeDef;'])
end
end
%% high-res
thisCaseName = ['constGrowth_metroMix_highRes'];
    try
        [resultMetroMixHR pbeDefMetroMixHR] = loadCase(thisCaseName);
    catch
        pbeDefMetroMixHR = obtainData(thisCaseName);
        resultMetroMixHR = executeCase(thisCaseName);
    end
%% analytic results (must be generated beforehand)
anaStruct = prepareConstGrowth(1, 'volume');

%% selection
% result = resultNormal;pbeDef = pbeDefNormal;
% result = resultMetroNum;pbeDef = pbeDefMetroNum;
result = resultMetroVol;pbeDef = pbeDefMetroVol;
% result = resultMetroMix;pbeDef = pbeDefMetroMix;

%% total volume
X = zeros(10, 31);
Y = zeros(10, 31);
YE = zeros(10, 31);
for iSeed = 1:10
    [Xt, Yt] = extractDataFromPbeResult('integralProperty', 'volume', ...
        result{iSeed}, pbeDef);
    X(iSeed, :) = Xt;
    Y(iSeed, :) = Yt(1,:);
    YE(iSeed, :) = Yt(2,:);
end

figure,hold on
plot(X', Y', 'k-')%, plot(X, Y-YE, 'r.'), plot(X, Y+YE, 'r.')
plot(X, anaStruct.value, 'r-', 'lineWidth', 3)
title('total volume');
% legend({'total volume (left axis)', 'lower error (std)', 'upper error (std)'},...
%     'Location', 'NorthWest');
% set(gca, 'yscale', 'log', 'xscale', 'log')
xlabel('time in s')
ylabel('total volume in m^3')
% error extra y-axis
axisHandle = axes;hold on
plot(X', YE'./Y', 'k--', 'parent', axisHandle)
plot(X', abs(bsxfun(@rdivide, Y, anaStruct.value) - 1)', 'r--', 'parent', axisHandle)
plot(X(1,:), std(Y,[],1)./mean(Y, 1), 'r-', ...
    'parent', axisHandle, 'lineWidth', 3)
plot(X(1,:), std(Y,[],1)./anaStruct.value, 'g-.', ...
    'parent', axisHandle, 'lineWidth', 3)
% plot([0 X(end)], [0 0], 'm-', 'parent', axisHandle, 'lineWidth', 3)
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])
% legend('coef. of variation (right axis)', 'Location', 'NorthEast')
ylabel('coefficient of variation')

%% total volume (convergence)
[X, Y] = extractDataFromPbeResult('integralProperty', 'volume', ...
    resultMetroMixHR, pbeDefMetroMixHR);
    
figure,hold on
plot(X, Y(1,:), 'k-')%, plot(X, Y-YE, 'r.'), plot(X, Y+YE, 'r.')
% plot(X, anaStruct.value, 'r-')
title('total volume');
% 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;hold on
plot(X, Y(2,:)./Y(1,:), 'k--', 'parent', axisHandle)
plot(X, abs(Y(1,:)./anaStruct.value - 1), 'r--', 'parent', axisHandle)
% plot([0 X(end)], [0 0], 'm-', 'parent', axisHandle, 'lineWidth', 3)
set(axisHandle, 'color', 'none', 'yaxisLocation', 'right', ...
    'xTick', [])
% legend('coef. of variation (right axis)', 'Location', 'NorthEast')
ylabel('coefficient of variation')