% this script loads model fits for each patient (first two weeks) and
% outputs patient_pain_reports
clear;
load('PatientData_Final.mat');
load('all_fitted_params_fixk0.mat')
% matrix to hold residuals, AIC, log likelihood, white noise fit
all_patient_params = [params, zeros(39,9)];

for ii=1:39
    ii
    % load correct patient information
    filename = patient_IDs{ii};
    patient_data = eval(patient_IDs{ii});

    M = csvread(strcat(filename,'_fixk0.csv'));
    patient_data = [patient_data(:,1),patient_data(:,10),patient_data(:,30),patient_data(:,31),patient_data(:,32)];
    % for patients with initial weird report, skip first report
    if (strcmp(filename,'PatientA0004') || strcmp(filename,'PatientA0009')) || (strcmp(filename,'PatientB1009') || strcmp(filename,'PatientB1010'))
        % FOR PATIENTA0004/9, PATIENTB1009/10 ONLY!!! Delete first report time (outlier)
        patient_data(:,1) = patient_data(:,1) - patient_data(2,1);
        patient_data(1,:) = [];
    end
    % only look at first 2 weeks of reports (to speed up test and avoid time-dependence)
    patient_data = patient_data(patient_data(:,1)<336,:);

    num_params = length(M);
    uhat = M(1);
%     k0hat = M(2);
    epshat = M(2);
    khat = M(3:num_params-1)';
    err = M(end);
    patient_info = [nanmean(patient_data(:,2)), uhat, epshat, log(2)/0.5, khat];
    drug_info = [log(0.5)/(-12.0),1; log(0.5)/(-4.0),1; log(0.5)/(-1.0),1];
    % Set drug times for each class of drug (will give empty vector if no drugs taken)
    la = patient_data(:,3); sa = patient_data(:,4); no = patient_data(:,5);
    la_drugtimes = patient_data(la>0,1); la_numdrugs = length(la_drugtimes); dose_times1 = la_drugtimes;
    sa_drugtimes = patient_data(sa>0,1); sa_numdrugs = length(sa_drugtimes); dose_times2 = sa_drugtimes;
    no_drugtimes = patient_data(no>0,1); no_numdrugs = length(no_drugtimes); dose_times3 = no_drugtimes;
    pain_report_real = patient_data(:,1:2)';
    pain_report_real(:,isnan(pain_report_real(2,:))) = [];
    % NEW AS OF 9/18/16 - SEND ONLY RELEVANT DRUG INFO (ONES ACTUALLY TAKEN)
    if no_numdrugs < 1
        drug_info(3,:) = [];
    end
    if sa_numdrugs < 1
        drug_info(2,:) = [];
        % shift untaken sa drug to last spot
        dose_times2 = no_drugtimes;
        dose_times3 = sa_drugtimes;
    end
    if la_numdrugs < 1
      %  patient_info(5) = [];
        drug_info(1,:) = [];
        % shift untaken la drug to last spot
        dose_times1 = dose_times2;
        dose_times2 = dose_times3;
        dose_times3 = [];
    end
    % vector containing number of drug doses for each class
    num_drug_vector = [la_numdrugs, sa_numdrugs, no_numdrugs];
    % total number of drugs taken
    num_drugs_taken = (la_numdrugs>0)+(sa_numdrugs>0)+(no_numdrugs>0);

    tvec = (patient_data(1,1):0.5:(patient_data(end,1)+1))';

    report = run_painsim_drugs3_ode(tvec,patient_info,drug_info,...
             dose_times1,dose_times2,dose_times3,pain_report_real);
     
    %csvwrite(strcat(filename,'_report_fixk0.csv'),report); 
    
    % calculate AIC for model (up to 5 tuning parameters)
    residuals = report(2,:)-report(3,:); % real - sim pain reports
    % Now AIC calculation: Assume Gaussian centered on deterministic model
    % prediction.
    sigma = nanstd(residuals);
    mu = 0; % expect mean of residuals to be zero for ideal model
    loggaussfunc = @(x)(-log(2*pi)/2-log(sigma)-(mu-x).^2./(2*sigma^2));
    loglike = sum(loggaussfunc(residuals));
    % The AIC is defined as 2k - 2 log L, where k is the number of
    % parameters and L is the likelihood.:
    k = 1+length(khat); % number of parameters in model (u+drug parameters)
    n = length(residuals); % number of data points for this patient
    AIC = 2*k - 2*loglike;
    if (n>=k+1)
        AICc = AIC + 2*k*(k+1)/(n-k-1);
    else
        AICc = nan;
    end
    
    
    % calculate AIC for white noise model (1 tuning parameter)
    res_noise = report(2,:)-mean(report(2,:)); % real - average pain reports
    % Now AIC calculation: Assume Gaussian centered on deterministic model
    % prediction.
    sigma_noise = nanstd(res_noise);
    mu = 0; % expect mean of residuals to be zero for ideal model
    loggaussfunc_noise = @(x)(-log(2*pi)/2-log(sigma_noise)-(mu-x).^2./(2*sigma_noise^2));
    loglike_noise = sum(loggaussfunc_noise(res_noise));
    % The AIC is defined as 2k - 2 log L, where k is the number of
    % parameters and L is the likelihood.:
    k_noise = 0; % number of parameters in model (u)
    n_noise = length(res_noise); % number of data points for this patient
    AIC_noise = 2*k_noise - 2*loglike_noise;
    if (n_noise>=k_noise+1)
        AICc_noise = AIC_noise + 2*k_noise*(k_noise+1)/(n_noise-k_noise-1);
    else
        AICc_noise = nan;
    end
    
    % store all calculations 
    % sigma should be same as error in col 7, give or take rounding errors
    all_patient_params(ii,7) = sigma;
    all_patient_params(ii,8) = loglike;
    all_patient_params(ii,9) = k;
    all_patient_params(ii,10) = AIC;
    all_patient_params(ii,11) = AICc;
    all_patient_params(ii,12) = sigma_noise;
    all_patient_params(ii,13) = loglike_noise;
    all_patient_params(ii,14) = AIC_noise;
    all_patient_params(ii,15) = AICc_noise;
    
    % rearrange k_i values to align with known drugs taken
    k_vals = zeros(1,3)*nan;
    k_vals(find(num_drug_vector>0)) = khat;
    all_patient_params(ii,3:5) = k_vals;
    
end

%% plots to compare model versus white noise (null)
figure
plot(all_patient_params(:,10), all_patient_params(:,14),'o')
xlabel('AIC full model')
ylabel('AIC white noise')
grid on

figure
plot(1:39,all_patient_params(:,10)-all_patient_params(:,14),'o')
xlabel('patient number')
ylabel('AIC model - white noise')

figure
plot(1:39,all_patient_params(:,11)-all_patient_params(:,15),'o')
xlabel('patient number')
ylabel('AIC model - white noise')