%% Passive pressure volume curve generation
% This function plots the passive pressure volume curve for all studies in 
% the analysis. 
% Author: ZJW
% Date: 21 Feb 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define colours
close all
clear all

green = [112 173 71]/255;
dark_green = [0, 149, 67]/255;
blue = [91 155 213]/255;
dark_blue = [0, 71, 189]/255;
light_blue = [7, 185, 252]/255;
red = [234, 0, 52]/255;
dark_red = [178, 31, 53]/255;
orange = [237 125 49]/255;
eggplant = [130, 0, 172]/255;
purple = [182, 16, 191]/255;
egg_yolk = [255, 206, 0]/255;
fuchsia = [255, 0, 128]/255;
charcoal = [30, 30, 30]/255;
black = [0,0,0]/255;
study_names =  {'STF_16', 'STF_17', 'STF_18','STF_19', 'STF_20',  'MR_250250'...
    , 'MR_293293', 'MR_119119', 'MR_054054', 'MR_104104', 'MR_236236', 'MR_269269','MR_087087', ...
    'MR_091091', 'MR_124124', 'MR_126126', 'STF_01', 'STF_02', 'STF_08', 'STF_09', 'STF_13', ...
    'STF_10', 'STF_12', 'MR_160160'};

figure_offset = figure;
hold on
figure_no_offset = figure;
hold on
count = 1;
dpdv= zeros(1, length(study_names));
fw = fopen('PassivedPdV.txt', 'w');
fprintf(fw, 'Study name\tdPdV linear\terror\tPV exponential b normalised\terror\n');
for i = study_names
    load([char(i), '_mri.mat']);
    load([char(i), '_lvp.mat']);
    % Get optimised LV wall volume. 
    dir = '/hpc/zwan145/Heart_Failure_Project_STF_Workspace/ParameterEstimation/PVSimulate/';
    opt_v = load([dir, char(i), '_volumes.txt']);

    % Get left ventricular wall volume. 
    fid = fopen(['../../CIM_Models/Studies/',mri.volume_path],'r');
    for j = 1:12
        line = fgetl(fid);
    end
    line = strsplit(line);
    V_w = str2num(line{6});
    figure(figure_offset);
    char(i)
    %h_offset(count) = plot(opt_v', output.LVP_average([mri.DS:end, mri.ED])*7.5, 'k.-', 'LineWidth', 1, 'MarkerSize', 10);
    if count < 6
        h_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_average([mri.DS:end, mri.ED])*7.5, '.-', 'Color', blue, 'LineWidth', 1.5, 'MarkerSize', 20);
    elseif count < 17 
        h_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_average([mri.DS:end, mri.ED])*7.5,'.-', 'Color', orange, 'LineWidth', 1.5, 'MarkerSize', 20);
    else 
        h_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_average([mri.DS:end, mri.ED])*7.5,'.-', 'Color', green, 'LineWidth', 1.5, 'MarkerSize', 20);
    end
    h_text = text(mri.V(mri.ED), output.LVP_average(mri.ED)*7.5, i);
    set(h_text, 'Interpreter', 'none', 'FontSize', 10);
    figure(figure_no_offset);
    h_no_offset(count) = plot(opt_v', output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5, 'k.-', 'LineWidth', 1, 'MarkerSize', 10);
    if count < 6
        h_no_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5,'.-', 'Color', blue, 'LineWidth', 1.5, 'MarkerSize', 20);
    elseif count < 17 
        h_no_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5,'.-', 'Color', orange, 'LineWidth', 1.5, 'MarkerSize', 20);
    else 
        h_no_offset(count) = plot(mri.V([mri.DS:end, mri.ED]), output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5,'.-', 'Color', green, 'LineWidth', 1.5, 'MarkerSize', 20);
    end
    %h_text = text(mri.V(mri.ED), output.LVP_no_offset_average(mri.ED)*7.5, i);
    [f, f_gof] = fit(mri.V([mri.DS:end, mri.ED])', output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5, 'poly1');
    %[fexp, fexp_gof] = fit(mri.V([mri.DS:end, mri.ED])', output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5, 'exp1');
    [fexp, fexp_gof] = fit(mri.V([mri.DS:end, mri.ED])', output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5, 'a + b*exp(c*(x-d))', 'StartPoint', [0,1,0.2,0], 'lower', [-10, 0.1, 0, -200]);
    
    %exp_curve = fexp.a * exp(fexp.b*mri.V([mri.DS:end, mri.ED])');
    exp_curve = fexp.a + fexp.b*exp(mri.V([mri.DS:end, mri.ED])' - fexp.c);
    lin_curve = f.p1*mri.V([mri.DS:end, mri.ED]) + f.p2;
    plot(mri.V([mri.DS:end, mri.ED])', lin_curve, '.-', 'Color', black, 'LineWidth', 1.0, 'MarkerSize', 10);
    opt_v*0.001
    plot(opt_v*0.001, output.LVP_no_offset_average([mri.DS:end, mri.ED])*7.5, '.-', 'Color', blue, 'LineWidth', 1.0, 'MarkerSize', 10);
    fprintf(fw, '%s\t%f\t%f\t%f\t%f\n', char(i), f.p1,f_gof.rmse, fexp.b*V_w, fexp_gof.rmse);
    %set(h_text, 'Interpreter', 'none', 'FontSize', 10);
    count = count + 1;
end
h_no_offset(count) = plot(mri.V([mri.DS:end, mri.ED])', lin_curve, '.-', 'Color', black, 'LineWidth', 1.0, 'MarkerSize', 10);
fclose(fw);
figure(figure_offset);
title('Passive pressure-volume curve');
ylabel('Pressure (mmHg)');
xlabel('Volume (mL)');
set(gca, 'DefaultTextInterpreter', 'none');
legend(h_offset([1, 2, 7, 18]), 'Model', 'Control', 'HFpEF', 'HFrEF', 'Location', 'northwest');

figure(figure_no_offset);
title('Passive pressure-volume curve');
ylabel('Pressure (mmHg)');
xlabel('Volume (mL)');
set(gca, 'DefaultTextInterpreter', 'none');
h_no_offset
legend(h_no_offset([1, 6, 17, 25]), 'Control', 'HFpEF', 'HFrEF','Fitted line', 'Location', 'northwest');

set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', 16);
set(gca,'DefaultAxesFontSize', 14);
set(gca,'DefaultTextFontSize', 8);
set(figure_offset, 'Position', [100, 100, 700, 450]);
set(figure_no_offset, 'Position', [100, 100, 700, 550]);

hold on 
[h] = plotyy(400, 5, 400, 5);
ylabel(h(2), 'Pressure (kPa)');
axis(h(2), [0 400 0 5.33333333]);
h(2).YTick = [0, 1, 2, 3, 4, 5];
h(1).YColor = 'black';
h(2).YColor = 'black';
%print(figure_offset, '-djpeg', 'PassivePressureVolume_offset');
%print(figure_no_offset, '-djpeg', 'PassivePressureVolume_no_offset');
%print(figure_offset, '-djpeg', '../../../Reports/HaemoAnalysisFigures/PassivePressureVolume_offset');
%print(figure_no_offset, '-djpeg', '../../../Reports/HaemoAnalysisFigures/PassivePressureVolume_no_offset');
%print(figure_no_offset, '-djpeg', '../../../../Writings/Manuscripts/JACCPassiveEstimationManuscript/PV_curve');
print(figure_no_offset, '-djpeg', 'PV_curve_with_model');
