%% This script finds parameters to fit model to patient data
clear;
% seed random number generator
rng(1);
% load patient data from mobile health app database
load('patient_data_all.mat');
% load modeling frameworks
load('model_names.mat')
% number of proposed models
num_models = length(model_names);
% select model to fit this round
model_index = 8;
% get the model name (e.g. 'threshold model')
modelname = model_names{model_index};
% number of patients
num_patients = length(patient_IDs);
% number of hours we want to inspect (336 is two weeks)
total_time = 336;
% time resolution for integration and plotting (0.5 is half hour)
time_res = 0.5;
% number of initial grid search locations for each patient
num_runs = 3;
% fixed k0 for patients (how fast patients rebound from reduced or increased pain, absent drugs)
% we take rebound half life to be half an hour
fixed_k0 = log(2)/0.5;
% indicator to merge all drug classes together (1 for merge, 0 for not)
merge_all_drugs = 0; common_drug_info = [log(0.5)/(-4.0),1];
% indicator of drug classes to include in fitting (LA,SA,NO) 
drugs_to_include = [1 1 1];
% threshold for drug dose times to be included in model (default is 1)
drug_threshold = 1;
% indicator that pain reporting probability is skewed towards higher pain
report_skew = 0;

% table to hold number of reports and drug doses within two week period
patient_table = zeros(num_patients, 4);

%% run through each patient's data
tic
for jj=1:num_patients
    % print patient number so we know where we're at
    jj
    % 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{jj};
    patient_data = load_patient_data_from_filename(filename,eval(patient_IDs{jj}),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,drug_threshold);
    
    % 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,:))) = [];
    
    % total number of drug classes taken (between 0 and 3; 0 or 1 if merged)
    num_drugs_taken = (~isempty(dose_times1))+(~isempty(dose_times2))+(~isempty(dose_times3));
    
    
    patient_table(jj,:) = [length(pain_report_real), length(la_drugtimes), length(sa_drugtimes), length(no_drugtimes)];
    
    
    
    
%     % use pseudo grid search to find best fit parameters for patient
%     [params,err] = find_patient_params(tvec,dose_times1,dose_times2,dose_times3,pain_report_real,patient_data,drug_info,fixed_k0,num_drugs_taken,num_runs,report_skew,nan);
%     % best fit patient info
%     patient_info = [nanmean(patient_data(:,2)),params(1),params(2),fixed_k0,params(3:length(params)-1)'];
%     % generate patient report (both real and simulated pain
%     report = run_painsim_drugs3_ode(tvec,patient_info,drug_info,dose_times1,dose_times2,dose_times3,pain_report_real);
%     
%     % plot results of model fit -> comment out if no plot desired
%     plot_model_fit(tvec,params,err,patient_data,drug_info,dose_times1,dose_times2,dose_times3,pain_report_real,fixed_k0);
%     
%     % plot raw data (pain reports and drugs taken) -> comment out if you don't want to see the raw time series data
%     %plot_raw_data(patient_data);
%     
%     % save results (if desired) -> comment out if you don't want to overwrite other csv or eps data
%     csvwrite(strcat(filename,'_fixk0.csv'),params);
%     csvwrite(strcat(filename,'_report_fixk0.csv'),report);
%     saveas(gcf,strcat(filename,'_fixk0.eps'),'eps2c'); close
end
toc

%% If you didn't save report csv files, use this
%calculate_aic(total_time,time_res,fixed_k0,merge_all_drugs,common_drug_info,drugs_to_include)

% reset number of models during testing - TEMP
num_models = 7;
% information criterion for all models and patients
info_crit_all_models = zeros(num_models,num_patients,9);
% cycle through models to compare AIC
for ii=1:num_models
    % collect information criterion for both model and white noise
    [info_crit_all_models(ii,:,:)] = calculate_aic_given_report(ii,0);
end

% plot model AIC - white noise AIC for all proposed models
figure
plot(1:num_patients,info_crit_all_models(:,:,4),'o-')
hold on
plot(1:num_patients,zeros(num_patients,1),'k-')
xlabel('patient number')
ylabel('AIC model - white noise')
legend('full','null','merge','LA','SA','NO','threshold')
%-info_crit_all_models(:,:,8)

% plot model AICc - white noise AICc for all proposed models
figure
plot(1:num_patients,info_crit_all_models(:,:,5),'o-')
hold on
plot(1:num_patients,zeros(num_patients,1),'k-')
xlabel('patient number')
ylabel('AICc model - white noise')
legend('full','null','merge','LA','SA','NO','threshold')