function plotIL4ModelResults()
% Use fitted parameter values to simulate all experimental
% conditions and plot together with experimental data

% Load fitted parameters
load('fittedParameters.mat')

% Load simbiology model
proj = sbioloadproject('il4ModelDetailed.sbproj');
m = proj.m1;

% Read experimental data from the respective files
[pStatTimeExp, pStatMean, pStatSE] = readData('pStat6Data.dat');
[tStatTimeExp, tStatMean, tStatSE] = readData('totalStat6Data.dat');
[socsTimeExp, socsMean, socsSE] = readData('socs3Data.dat');

% List of pretreatment conditions and stimuli
pretreat = {'PBS', 'MG132', 'CHX', 'DPI'}; 
stim = {'IL4'}; 

% This dictionary maps experimental condition to the appropriate index in
% the experimental data cell arrays, e.g., the mean values of pStat under
% experimental condition correspoding to DPI+IL4 can be accessed using
% pStatTimeMean{pStatIdx.('DPI_IL4')}
pStatIdx = struct('PBS_PBS', 1, 'PBS_IL4', 2, 'DPI_PBS', 3, 'DPI_IL4', 4, ...
    'MG132_PBS', 5, 'MG132_IL4', 6, 'CHX_PBS', 7, 'CHX_IL4', 8, ...
    'PBS_H2O2', 9, 'PBS_H2O2IL4', 10);
tStatIdx = struct('PBS_PBS', 1, 'PBS_IL4', 2, 'MG132_PBS', 3, 'MG132_IL4', 4);
socsIdx = struct('PBS_PBS', 1, 'PBS_IL4', 2, 'CHX_PBS', 3, 'CHX_IL4', 4);

pStatFig = figure; hold on;
tStatFig = figure; hold on;
socsFig = figure; hold on;

% Simulate all pretreatment and stimulation conditions
for i = 1:length(pretreat)
    for j = 1:length(stim)
        cond = [pretreat{i} '_' stim{j}];
        disp(['Simulating ' pretreat{i} ' + ' stim{j}]);
        
        % Simulate model and get output
        try
            [t, pStatModel, tStatModel, socsModel] = simulateModel(...
                m, fittedParam, pretreat{i}, stim{j});
        catch err
            disp(err);
            return;
        end
        
        if t(end) < 120 %simulation stopped before completing full time span
            error('Failed to integrate over 2 hr period');
        end

        % Set scaling factors for all output variables. The values are
        % estimated to optimize the fits visually for training data.
        if strcmp(pretreat{i}, 'PBS') && strcmp(stim{j}, 'IL4')
            [~, tStata, tStatb] = scale2data(tStatModel, tStatMean{2});
            [~, socsa, socsb] = scale2data(socsModel, socsMean{2});
            socsb = socsb - 0.04; % Adjusting for better visual fit
            pStata = 1/pStatModel(1);
            pStatb = 0;
        end
        
        % Scale model output using the estimated scaling factos
        pStatModel = pStata*pStatModel + pStatb;
        tStatModel = tStata*tStatModel + tStatb;
        socsModel = socsa*socsModel + socsb;
        
        % Add scaled model output to appropriate figure along with
        % experimental data depending on pretreament condition
        plotData(pStatFig, 'pStat6', cond, t, pStatModel, ...
            pStatTimeExp{pStatIdx.(cond)}, pStatMean{pStatIdx.(cond)}, ...
            pStatSE{pStatIdx.(cond)});
        switch pretreat{i}
            case 'PBS' %additional error for total Stat and socs
                plotData(tStatFig, 'tStat6', cond, t, tStatModel, ...
                    tStatTimeExp{tStatIdx.(cond)}, ...
                    tStatMean{tStatIdx.(cond)}, tStatSE{tStatIdx.(cond)});                
                plotData(socsFig, 'socs', cond, t, socsModel, ...
                    socsTimeExp{socsIdx.(cond)}, socsMean{socsIdx.(cond)}, ...
                    socsSE{socsIdx.(cond)});
            case 'MG132' %add error for total Stat
                plotData(tStatFig, 'tStat6', cond, t, tStatModel, ...
                    tStatTimeExp{tStatIdx.(cond)}, ...
                    tStatMean{tStatIdx.(cond)}, tStatSE{tStatIdx.(cond)});
           case 'CHX'
                plotData(socsFig, 'socs', cond, t, socsModel, ...
                    socsTimeExp{socsIdx.(cond)}, socsMean{socsIdx.(cond)}, ...
                    socsSE{socsIdx.(cond)});
       end            
    end
end

% Special case of PBS + H2O2/IL4 treatment
cond = 'PBS_H2O2IL4';
disp('Simulating PBS + H2O2/IL4');
try
    [t, pStatModel, ~, ~] = simulateModel(...
        m, fittedParam, 'PBS', 'H2O2IL4');
catch err
    error(err);
end
if t(end) < 120 %simulation stopped before completing full time span
    error('Failed to integrate over 2 hr period');
end
% Scale output and add to figure
pStatModel = pStata*pStatModel + pStatb;
plotData(pStatFig, 'pStat6', cond, t, pStatModel, ...
    pStatTimeExp{pStatIdx.(cond)}, pStatMean{pStatIdx.(cond)}, ...
    pStatSE{pStatIdx.(cond)});
figure(socsFig); legend('show');
figure(tStatFig); legend('show');
figure(pStatFig); legend('show');
%---


function plotData(figH, name, cond, tModel, yModel, tExp, yExp, seExp)
% Add output data to figure fighH

persistent c;  % Color index
if isempty(c)
    c = 1;
end
% Cycle through these colors
col = [ 1 0 0;
        0 1 0;
        0 0 1;
        0 0 0;
        .5 .5 0;
        0 .5 .5;
        .5 0 .5;
        1 .5 0;
        0 1 .5;
        1 0 .5;
        .5 1 0;
        0 .5 1;
        .5 0 1];

% Open figH, add model output and experimental data
figure(figH);
l = plot(tModel, yModel, 'LineWidth', 3);
e = errorbar(tExp, yExp, seExp);
set(e, 'Color', col(c, :), 'LineStyle', 'none', 'MarkerSize', 8, 'LineWidth', 2);
set(l, 'DisplayName', [name ' (' strrep(cond, '_', '+') ')'], 'Color', col(c, :));

% Cycle to next color
c = c + 1; 
if c > length(col), c = 1; end


function [Y, a, b] = scale2data(y, dataPts)
% Scales the curve y linearly so that the max(y) = max(dataPts) and 
% min(y) = min(dataPts). This means that the dataPts and the curve y can
% be fitted with the same bounds on the y axis. Note that max/min(y) may
% not occur at the same time point at max/min(dataPts). Therefore, this is
% not forcing any points on the two data sets to coincide all the time. 
% Y = scaled curve

% Scaling: Y = a*y + b
% Condition 1: Ymin = dataPtmin = a*ymin + b
% Condition 2: Ymax = dataPtmax = a*ymax + b

Ymax = max(dataPts); Ymin = min(dataPts);
ymax = max(y); ymin = min(y);

a = (Ymax - Ymin)/(ymax - ymin);
b = (Ymax*ymin - Ymin*ymax)/(ymin - ymax);

Y = a*y + b;
