% read in patient info
M = csvread('patient_info.csv');
% rate of decrease of pain without drugs or disease
relax_rate = M(1);
% baseline pain, equilibrium for patient in absence of drugs
base_pain = M(2);
% initial pain
init_pain = M(3);

% read in simulation info
M = csvread('sim_info.csv');
% initial time
t0 = M(1);
% end time
tmax = M(2);
% time resolution
dt = M(3);
% time vector
tvec = t0:dt:tmax;

% read in drug info
M = csvread('drug_info.csv');
% marginal effect of drug 1 on pain relaxation rate
k1 = M(1);
% rate of elimination of drug 1
kD1 = M(2);
% amplitude of response to drug 1
a1 = M(3);
% drug information
drug_info1 = [k1, kD1, a1];
%     % number of doses
%     num_doses = 5;
%     wait = (tmax-t0)/(num_doses);
%     drug 1 dose times
%     dosetimes1 = [0.1:wait:tmax]
% read in drug dosage
dosetimes1 = csvread('dose_times.csv');

% read in "real" patient data
M = csvread('pain_report_real.csv');
report_times_real = M(1,:);
report_pain_real = M(2,:);

%% run pain simulation with no drug intervention
% [tvec, P] = pain_sim_nointervention(tvec, init_pain, relax_rate, base_pain);
% figure;
% plot(tvec, P, '-r');
% xlabel('Time');
% ylabel('Pain level');
% title('Pain versus time (no intervention)');

%% run pain simulation with one drug intervention
[tvec, P, D] = pain_sim_intervention1(tvec,init_pain,relax_rate,base_pain,drug_info1,dosetimes1);
figure;
plot(tvec, P, '-r');
hold on;
plot(tvec, D, '--b');
legend('Pain', 'Drug 1');
xlabel('Time');
ylabel('Pain and drug level');
title('Pain and drug concentration versus time (1 drug intervention)');

%% calculate average pain
% average_pain = sum(P)*dt/(tmax-t0)
% average_drug = sum(D)*dt/(tmax-t0)

%% report pain levels
report_pain_sim = interp1(tvec,P,report_times_real);
plot(report_times_real, report_pain_real, 'ok')
plot(report_times_real, report_pain_sim, 'xk')
legend('Pain', 'Drug 1','real pain report','sim pain report');
csvwrite('pain_report.csv',[report_times_real; report_pain_real; report_pain_sim])

%% calculate pain prediction error
stdev = sqrt(1/length(report_times_real)*sum((report_pain_real-report_pain_sim).^2));
dlmwrite('pain_report.csv',stdev,'-append')