function calculate_aic(total_time,time_res,fixed_k0,merge_all_drugs,common_drug_info,drugs_to_include)
% this function loads model fits for each patient and
% outputs patient_pain_reports
load('patient_data_all.mat');
num_patients = length(patient_IDs);
% matrix to hold residuals, AIC, log likelihood, white noise fit
% params is just a placeholder for now
all_fitting_error = zeros(num_patients,9);

for ii=1:num_patients
    % print patient number so we know where we're at
    ii
    % Long acting (t_halflife = 10 hr); Short acting (t_halflife = 4 hr); Non-opiod (t_halflife = 4 hr)
    drug_info = [log(0.5)/(-10.0),1; log(0.5)/(-4.0),1; log(0.5)/(-4.0),1];
    % load patient data into convenient form (just useful info on specified time)
    filename = patient_IDs{ii};
    patient_data = load_patient_data_from_filename(filename,eval(patient_IDs{ii}),total_time);
    
    % indicator time vectors of drugs taken (1 for drug taken, 0 for not)
    la = patient_data(:,3); sa = patient_data(:,4); no = patient_data(:,5);
    % drug time vectors for each class of drug (empty vector if no drugs taken)
    la_drugtimes = patient_data(la>0,1); sa_drugtimes = patient_data(sa>0,1); no_drugtimes = patient_data(no>0,1); 

    % collect and reorder drugs taken to eliminate those not taken; drug_info must also change to reflect reordering of drugs
    [dose_times1,dose_times2,dose_times3,drug_info] = get_dose_times(drug_info,la_drugtimes,sa_drugtimes,no_drugtimes,merge_all_drugs,drugs_to_include,common_drug_info);
    
    % time vector (add an extra hour to end for interp1 - it freaks out about rounding errors)
    tvec = (patient_data(1,1):time_res:(patient_data(end,1)+1))';
    % pain report contains time and pain info from patient (remove times when patient does not report actual pain)
    pain_report_real = patient_data(:,1:2)'; pain_report_real(:,isnan(pain_report_real(2,:))) = [];
    
    % read in model parameters for each patient
    M = csvread(strcat('Patient_fixk0 (merge all drugs)/',filename,'_fixk0.csv'));
    num_params = length(M);
    uhat = M(1); epshat = M(2); khat = M(3:num_params-1)';
    % set patient info provided by file
    patient_info = [nanmean(patient_data(:,2)),uhat,epshat,fixed_k0,khat];
    % generate pain report for this patient
    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); 

% THIS IS FOR GAUSSIAN PROB MODEL
    % 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));
% END GUASSIAN PROB MODEL
    
% % THIS IS FOR SKEWED PROB MODEL
%     real_pain = report(2,:); real_pain(real_pain==0) = 10^(-10);
%     sim_pain = report(3,:);
%     % calculuate negative log likelihood
%     % standard deviation of residuals
%     sig = std(real_pain-sim_pain);
%     % vector of all probabilities of reported pain given skewed normal dist
%     p_vec = real_pain.*exp(-(real_pain-sim_pain).^2/(2*sig*sig)) ./ ...
%             (sig*sig*exp(-(sim_pain).^2/(2*sig*sig)) + sqrt(pi/2)*sim_pain*sig.*...
%             (1.0+erf(sim_pain/(sqrt(2)*sig))));
%     loglike = sum(log(p_vec));
% % END SKEWED PROB MODEL 
    
    
    % 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 = 1; % 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_fitting_error(ii,1) = sigma;
    all_fitting_error(ii,2) = loglike;
    all_fitting_error(ii,3) = k;
    all_fitting_error(ii,4) = AIC;
    all_fitting_error(ii,5) = AICc;
    all_fitting_error(ii,6) = sigma_noise;
    all_fitting_error(ii,7) = loglike_noise;
    all_fitting_error(ii,8) = AIC_noise;
    all_fitting_error(ii,9) = AICc_noise;
    
end

%% plots to compare model versus white noise (null)
figure
plot(all_fitting_error(:,4), all_fitting_error(:,8),'o')
xlabel('AIC full model')
ylabel('AIC white noise')
grid on

figure
plot(1:num_patients,all_fitting_error(:,4)-all_fitting_error(:,8),'o')
hold on
plot(1:num_patients,zeros(num_patients,1),'k-')
xlabel('patient number')
ylabel('AIC model - white noise')

figure
plot(1:num_patients,all_fitting_error(:,5)-all_fitting_error(:,9),'o')
hold on
plot(1:num_patients,zeros(num_patients,1),'k-')
xlabel('patient number')
ylabel('AICc model - white noise')

end