% Function uses fminsearch to optimize the drug timing in the simulation
%     time period, for a given number of drug doses
% Option average pain plots for 2 drug doses
function [opt_times,ave_pain] = opt_drug1(num_doses)
% read in simulation info
sim_info = csvread('sim_info.csv');

% read in patient info
patient_info = csvread('patient_info.csv');

% read in drug info
drug_info = csvread('drug_info.csv');

% set IC to random dosage times
tmax = sim_info(2);
%dose_times = tmax*rand(1,num_doses);
dose_times = linspace(0,tmax-tmax/num_doses,num_doses);

func_handle = @(doses)(average_pain(doses,patient_info, sim_info, drug_info));
%[opt_times,~] = fminsearchbnd(func_handle, dose_times, [0,24],[1,23], optimset('TolX',1e-3));

% METHOD 1: OPTIMIZE OVER ALL TIME
%[opt_times,ave_pain] = fminsearch(func_handle, dose_times, optimset('TolFun',1e-2));
% opt_times = sort(opt_times);

% METHOD 2: OPTIMIZE OVER DISCRETE TIME SET
[opt_times,ave_pain] = fminsearch(func_handle, dose_times, optimset('TolFun',0.5,'TolX',0.5));

%% Check average pain function (assuming 2 doses, first dose time is 0)
% test_times = linspace(0,24,150);
% ave_pain = zeros(1,length(test_times));
% for ii=1:length(test_times) 
%     ave_pain(ii) = average_pain([0,test_times(ii)],patient_info, sim_info, drug_info);
% end
% 
% figure
% plot(test_times,ave_pain,'-o')

%% Investigate average pain contour plot for 2 drug doses
% t1_range = 0:0.2:24;
% t2_range = 0:0.2:24;
% ave_pain = zeros(length(t1_range),length(t2_range));
% for t1=1:length(t1_range)
%     for t2=1:length(t2_range)
%         ave_pain(t1,t2) = average_pain([t1_range(t1),t2_range(t2)],patient_info, sim_info, drug_info);
%     end
%     t1_range(t1)
% end
% 
% save ave_pain.mat
% 
% figure
% contourf(t1_range,t2_range,ave_pain,75)
% colorbar('location','southoutside')
% xlabel('dose time 1')
% ylabel('dose time 2')

end

%% run pain simulation with one drug intervention
function ave_pain = average_pain(dose_times,patient_info, sim_info, drug_info)
    % rate of decrease of pain without drugs or disease
    relax_rate = patient_info(1);
    % baseline pain, equilibrium for patient in absence of drugs
    base_pain = patient_info(2);
    % initial pain
    init_pain = patient_info(3);

    % initial time
    t0 = sim_info(1);
    % end time
    tmax = sim_info(2);
    % time resolution
    dt = sim_info(3);
    % time vector
    tvec = t0:dt:tmax;
    
    % Round dose times to nearest grid point to avoid numerical noise
    dose_times_new = t0+round((dose_times-t0)/dt)*dt;
    [tvec, P, ~] = pain_sim_intervention1(tvec,init_pain,relax_rate,base_pain,drug_info,dose_times_new);
    % calculate average pain
    ave_pain = trapz(tvec,P)/(tmax-t0);
    
%     figure
%     plot(tvec, P)
end