%% Development example
% Run solver and display results for an initially selected example
clear all
% close all

data = struct([]);
% data(end+1).caseName = 'roundingRatioCommon';
data(end+1).caseName = 'roundingRatio_kLow_aNeg';
data(end+1).caseName = 'roundingRatio_kHigh_aNeg';
data(end+1).caseName = 'roundingRatio_kLow_aZero';
data(end+1).caseName = 'roundingRatio_kHigh_aZero';
data(end+1).caseName = 'roundingRatio_kLow_aPos';
data(end+1).caseName = 'roundingRatio_kHigh_aPos';
obtainData('init')

% execute
for iCase = 1:length(data)
    pbeDefinition = obtainData(data(iCase).caseName);
    result = solvePbe(pbeDefinition);
    save(['data' filesep data(iCase).caseName '.mat'], 'result');
end

%% load
for iCase = 1:length(data)
    iCase
    data(iCase).pbeDefinition = obtainData(data(iCase).caseName);
    result = load(['data' filesep data(iCase).caseName '.mat']);
    data(iCase).result = result.result;
end

%% summarize/evaluate
for iCase = 1:length(data)
    % generate data
    [data(iCase).time, data(iCase).mw] = extractDataFromPbeResult(...
        'bulkProperty', 1, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).surf] = extractDataFromPbeResult(...
        'bulkProperty', 2, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).vol] = extractDataFromPbeResult(...
        'bulkProperty', 3, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).muB1] = extractDataFromPbeResult(...
        'bulkProperty', 4, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).muB2] = extractDataFromPbeResult(...
        'bulkProperty', 5, ...
        data(iCase).result, data(iCase).pbeDefinition);
    [data(iCase).time, data(iCase).muB3] = extractDataFromPbeResult(...
        'bulkProperty', 6, ...
        data(iCase).result, data(iCase).pbeDefinition);
    data(iCase).filter = data(iCase).mw > 0;

    % plot
    figure,hold on,
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).mw(data(iCase).filter), 'k-')
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).surf(data(iCase).filter)/data(iCase).surf(1), 'k-')
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).vol(data(iCase).filter)/data(iCase).vol(1), 'k-')
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).muB1(data(iCase).filter), 'm-')
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).muB2(data(iCase).filter), 'c-')
    plot(data(iCase).time(data(iCase).filter), ...
         data(iCase).muB3(data(iCase).filter), 'k-', 'lineWidth', 3)
    title(data(iCase).caseName)
    ylim([0 1])
    xlim([0 1])
end

save('data/roundingRatioSummary.mat', 'data')

