function [all_info_crit] = calculate_aic_given_report(model_index,report_skewed)
% this function loads model fits for each patient and
% outputs patient_pain_reports
load('patient_data_all.mat');
% load modeling frameworks
load('model_names.mat');
modelname = model_names{model_index};
num_patients = length(patient_IDs);
% matrix to hold residual info, AIC, log likelihood, white noise fit
all_info_crit = zeros(num_patients,9);

for ii=1:num_patients
    % load patient data into convenient form (just useful info on specified time)
    filename = patient_IDs{ii};
    
    % read in model parameters for each patient to assess degree of freedom
    M = csvread(strcat('Patient_fixk0_',modelname,'/',filename,'_fixk0.csv'));
    num_params = length(M); khat = M(3:num_params-1)';
    report = csvread(strcat('Patient_fixk0_',modelname,'/',filename,'_report_fixk0.csv'));

    % calculate likelihood and AIC given report,degrees of freedom, and model number
    % for both model and white noise
    [all_info_crit(ii,:)] = calculate_info_crit(report,1+length(khat),report_skewed);
    
end

%% plots to compare model versus white noise (null) -> comment out if no plots desired
% figure
% plot(1:num_patients,all_info_crit(:,4)-all_info_crit(:,8),'o')
% hold on
% plot(1:num_patients,zeros(num_patients,1),'k-')
% xlabel('patient number')
% ylabel('AIC model - white noise')
% title(modelname,'Interpreter', 'none')
% 
% figure
% plot(1:num_patients,all_info_crit(:,5)-all_info_crit(:,9),'o')
% hold on
% plot(1:num_patients,zeros(num_patients,1),'k-')
% xlabel('patient number')
% ylabel('AICc model - white noise')
% title(modelname,'Interpreter', 'none')

end

function [info_crit] = calculate_info_crit(report,num_params,report_skewed)
% extract real and simulated pain from reports
real_pain = report(2,:); real_pain(real_pain==0) = 10^(-10);
sim_pain = report(3,:);
% calculate residuals
residuals = real_pain - sim_pain; % real - sim pain reports
% if the models assumes equal probability of any particular pain value
% being reported by patient...
if report_skewed==0
    % calculate AIC for provided model
    [sigma,loglike,AIC,AICc] = cic(residuals,num_params);
    
    % calculate AIC for white noise model (1 tuning parameter)
    res_noise = real_pain-mean(real_pain); % real - average pain reports
    [sigma_noise,loglike_noise,AIC_noise,AICc_noise] = cic(res_noise,1);

% if the model assumes that patients are more likely to report higher pain
% levels...
else
    % calculate AIC for provided model
    [sigma,loglike,AIC,AICc,fsolve_exitflag] = cic_skew(real_pain,sim_pain,num_params);
    if fsolve_exitflag~=1
        fsolve_exitflag
    end
    
%     % calculate AIC for white noise model (1 tuning parameter)
%     [sigma_noise,loglike_noise,AIC_noise,AICc_noise,fsolve_exitflag_noise] = cic_skew(real_pain,ones(1,length(real_pain))*mean(real_pain),1);
end

info_crit = zeros(1,9);
% store all calculations 
info_crit(1) = num_params;
info_crit(2) = mean(sigma);
info_crit(3) = loglike;
info_crit(4) = AIC;
info_crit(5) = AICc;
% info_crit(6) = mean(sigma_noise);
% info_crit(7) = loglike_noise;
% info_crit(8) = AIC_noise;
% info_crit(9) = AICc_noise;

end

function [sigma,loglike,AIC,AICc] = cic(residuals,k)
% 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.
% number of parameters in model (k = 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
end

function [sigma,loglike,AIC,AICc,fsolve_exitflag] = cic_skew(real_pain,sim_pain,num_params)
% calculuate negative log likelihood
% weighting to determine p_r(PAIN=0) and p_r(PAIN=10)
a = 0.08; b = 0.1;
% turn pain=0 into pain=really small to avoid errors
real_pain(real_pain==0) = 10^(-10);
mu = sim_pain;
mu = mu'; real_pain = real_pain';
[~, ~, sigma, fsolve_exitflag] = solve_mu_sigma(mu, real_pain, a, b);
% vector of all probabilities of reported pain given skewed normal dist   
p_vec = (a*real_pain + b).*exp(-((real_pain-mu).^2)./(2*sigma.*sigma))./...
        (a*exp(-(mu.*mu./(2*sigma.*sigma))).*sigma.*sigma +... 
        sqrt(pi/2)*(b+a*mu).*sigma.*(1 + erf(mu./(sqrt(2)*sigma))));
% log likelihood    
loglike = sum(log(p_vec));

% The AIC is defined as 2k - 2 log L, where k is the number of
% parameters and L is the likelihood.:
k = num_params; % number of parameters in model (u+drug parameters)
n = length(real_pain); % 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
end


%%% function solves system of equations for skewed prob. 
%%%       mu_r = f_1(mu, sigma)
%%%  sigma^2_r = f_2(mu, sigma)
%%%  sigma^2_r = g(mu_r, real_pain)
%%%
%%%  real_pain is a vector of length N
%%%  mu is a vector of length N
%%%  mu_r is a vector of length N
%%%  sigma_r is a scalar (assumed constant over time)
%%%  sigma is a vector of length N
function [mu_r, sigma_r, sigma, exitflag] = solve_mu_sigma(mu, real_pain, a, b)
% number of pain reports
N = length(real_pain);
% function for mu_r (biased reported pain mean)
f1_handle = @(mu,sigma) (exp(-(mu.*mu./(2*sigma.*sigma))).*(b+a*mu).*(sigma.*sigma)+... 
                        sqrt(pi/2)*sigma.*(b*mu+a*(mu.*mu+sigma.*sigma)).*...
                        (1 + erf(mu./(sqrt(2)*sigma))))./(a*exp(-(mu.*mu./(2*sigma.*sigma))).*...
                        sigma.*sigma + sqrt(pi/2)*(b+a*mu).*sigma.*(1+erf(mu./(sqrt(2)*sigma))));
% function for sigma^2_r (biased reported pain variance - assumed constant in time)
f2_handle = @(mu,sigma) ((sigma.^3).*(exp(-(mu.*mu./(sigma.*sigma))).*(-2*b*(b+a*mu).*...
                        sigma+4*a*a*sigma.^3) + pi*sigma.*((b+a*mu).^2 -... 
                        a*a*sigma.*sigma).*(1+erf(mu./(sqrt(2)*sigma))).^2 +... 
                        exp(-(mu.*mu./(2*sigma.*sigma))).*sqrt(2*pi).*(b*mu.*(b+a*mu) -... 
                        a*(b + 3*a*mu).*sigma.*sigma).*(-2 + erfc(mu./(sqrt(2)*sigma)))))./...
                        (2*(a*exp(-(mu.*mu./(2*sigma.*sigma))).*sigma.*sigma +... 
                        sqrt(pi/2)*(b+a*mu).*sigma.*(1+erf(mu./(sqrt(2)*sigma)))).^2);
% function for sigma^2_r (definition of variance)
g_handle = @(mu_r, real_pain) 1.0/N*sum((real_pain-mu_r).^2);

% solve for sigma at each report time using mu at each time (known)
% x contains sigma at each report time
system1_handle = @(x) f2_handle(mu,x) - ones(N,1)*g_handle(f1_handle(mu,x), real_pain);
% guess IC near input for sigma_r
x0 = ones(N,1)*g_handle(mu,real_pain);
% solve for zero
options = optimoptions('fsolve','Display','none');
[sigma,~,exitflag] = fsolve(system1_handle,x0,options);
% solve directly for mu_r
mu_r = f1_handle(mu,sigma);
% solve directly for sigma_r (assumed constant in time)
sigma_r = sqrt(g_handle(mu_r, real_pain));
end