% DP Solver Utility - Sensitivity Analysis (Created 09/20/18; Last 
% Modified 04/23/19)

%% Initialization

% NOTE 090918: Stop execution of code by typing CTRL + C.
% NOTE 091518: "Alpha" denotes gender (i.e., male and female); "beta"
% denotes ethnicity (i.e., Afro-Caribbean, Asian-Indian, and Caucasian);
% and "gamma" denotes education level (i.e., TO COME).

clearvars; clc;

LM = "_042019";
BASE = regexprep(matlab.desktop.editor.getActiveFilename, ...
    strcat('Diabetes_Model', LM, '.*$'), ...
    strcat('Diabetes_Model', LM));

Ref_Path = strcat(BASE, "\R_Objects\");
MAT_Path = strcat(BASE, "\MATLAB\");

Sens_Analysis_Types = ["1_Functional_Form", "2_Alpha", ...
    "3_Discount_Factor", "4_Price_Consump", "5_Wage", "6_ARPD", ...
    "7_Savings_Percentage", "8_Dietary_Costs"];

% NOTE 092018: SENS_ANALYSIS_ID refers to the parameter that is currently
% being varied. The vector Person_Type indicates the gender (i.e., male or 
% female), ethnicity (i.e., Afro-Caribbean, Asian-Indian, and Caucasian), 
% and level of educational attainment (i.e., No certificate, diploma or
% degree; High school diploma or equivalent; Postsecondary certificate or 
% diploma below bachelor level; University certificate, diploma or degree
% at bachelor level or above; and, Postsecondary certificate, diploma or 
% degree) of the decision-maker.
Sens_Analysis_ID = 5; Person_Type = [2, 3, 4];

Table_1_Names = cellstr({"Sens_Analysis_ID", "Sens_Analysis_Types", ...
    "Start", "End"});
Table_1 = array2table([(1:8)' transpose(Sens_Analysis_Types) ...
    transpose(repelem(1, 8)) [1 1 1 3 3 1 1 1]']);
Table_1.Properties.VariableNames = Table_1_Names;
Table_1.Sens_Analysis_ID = str2double(Table_1.Sens_Analysis_ID);
Table_1.Start = str2double(Table_1.Start);
Table_1.End = str2double(Table_1.End);

% This table refers to values of the substitution parameter incorporated
% into the sensitivity analysis.
Table_2_Names = cellstr({"Rho", "Value"});
Table_2 = array2table([(1:3)' [-10 0 1]']);
Table_2.Properties.VariableNames = Table_2_Names;

% The discount factor is the parameter that will actually need to be varied
% in M-File itself. Variation in all other parameters is handled in the
% corresponding R scripts.
Disc_Factor = [0.00, 0.50, 0.75, 0.95, 1.00];

% https://www.mathworks.com/help/matlab/ref/readtable.html
% https://www.mathworks.com/help/matlab/ref/textscan.html
% https://www.mathworks.com/matlabcentral/answers/291406-how-to-specify-the-format-of-just-one-column-with-readtable
RefActions = readtable(strcat(Ref_Path, "Ref_Actions.dat"), ...
    'Format', '%u%s%s%s%s');
% https://www.mathworks.com/help/matlab/matlab_prog/clean-messy-and-missing-data-in-tables.html
RefStates = readtable(strcat(Ref_Path, "Ref_States.dat"), ... 
    'TreatAsEmpty',{'.','NA'}, 'Format', '%u%s%s%s%s%s%s%s%s%s%s');

% MDPSolve won't work unless the appropriate folders are included in the
% user's MATLAB path. The following line of code adds MATLAB and all its
% subfolders to the path.
% https://www.mathworks.com/matlabcentral/answers/247180-how-may-i-add-all-subfolders-in-my-matlab-path
addpath(genpath("C:\Users\starg\Documents\MATLAB\"))

for z = 1:5
for a = Table_1{Sens_Analysis_ID, 'Start'}:Table_1{Sens_Analysis_ID, 'End'}
%% Import data from R

% Alpha, beta, and gamma are all integers; before we can concatenate 
% them, we must therefore first convert each from an integer to a string.
Case = strcat(num2str(Person_Type(1)), num2str(Person_Type(2)), ...
    num2str(Person_Type(3)));

% Generate a cell array to hold the results from the XLSX files.
% https://www.mathworks.com/help/matlab/ref/cell.html
TP_List = cell(24, 1); RWD_List = cell(24,1);

% https://www.mathworks.com/help/matlab/matlab_prog/access-data-in-a-cell-array.html
% https://stackoverflow.com/questions/12164752/matlab-convert-array-of-number-to-array-of-strings
TP_01 = (1:24); TP_02 = transpose(strtrim(cellstr(num2str(TP_01'))'));
RWD_01 = (1:24); RWD_02 = transpose(strtrim(cellstr(num2str(RWD_01'))'));

% NOTE 042319: Modified code surrounding import of TP matrices. This was
% motivated by the observation that we do not actually need to produce TP
% matrices for each run of every sensitivity analysis since no element of
% the matrices actually change from one run to the next. Consequently, we
% can simply recycle the same matrices used in the baseline analysis, which
% are of course stored in the baseline folder.
% 
% ARRAYFUN applies a function to each element of an array. In this case,
% all we're doing is concatenating strings to produce a series of file
% pathways/names.
TP_03 = arrayfun(@(x) ...
    strcat(MAT_Path, "Baseline_Results", ...
    "\TP_Matrix_Components\Case_", Case, ...
    "\TP_", x, "_", Case, ".csv"), TP_02);

if (Sens_Analysis_ID == 4 || Sens_Analysis_ID == 5)
    Sens_Anal_String = strcat("\", "Rho_", num2str(Table_2{a, 'Value'}));
else
    Sens_Anal_String = [];
end
    
RWD_03 = arrayfun(@(x) ...
    strcat(MAT_Path, "Sensitivity_Analysis\", ...
    Sens_Analysis_Types(Sens_Analysis_ID), ...
    "\RWD_Matrix_Components\Case_", num2str(z), Sens_Anal_String, ...
    "\RWD_", x, "_", Case, ".csv"), RWD_02);

%% Read TP & RWD Matrices from XLS

% https://www.mathworks.com/help/matlab/ref/arrayfun.html
% https://www.mathworks.com/matlabcentral/answers/274210-xlsread-returns-empty-when-reading-a-cell-with-a-formula
% https://www.mathworks.com/matlabcentral/answers/103299-convert-multiple-csv-to-xls-files

% Create an empty cell array to hold the components of the TP and RWD
% matrices.
TP_04 = cell(24, 1); RWD_04 = cell(24, 1);

% Read the components of the matrices in sequence, assigning them to
% "slots" in the empty cell array.
TP_04 = arrayfun(@(x) xlsread(char(x)), TP_03, 'UniformOutput', false);
RWD_04 = arrayfun(@(x) xlsread(char(x)), RWD_03, 'UniformOutput', false);

% NOTE 091018: While now a relic, the next three lines of code illustrate
% something important about the XLSREAD function, namely that it preserves
% only numeric values unless specifically instructed to do otherwise. In
% particular, bifurcating output as is done below serves to segment it into
% numeric (the 1st argument on the LHS) and text elements (the 2nd
% argument). The clue to the issues I had been experiencing with the
% XLSREAD function (documented in greater detail in the accompanying R
% script) came about when the content of the sole TP matrix component I had
% manually formatted as numeric data type was the only non-null element
% in the first of the two sets of outputs generated through the code below.
% for i = 1:24
%    [TP_04{i,1}, TP_04b{i,1}] = xlsread(TP_03{i,1});
% end

% NOTE 091218: When running the model I receive the following message:
% "P has rows that do not sum to 1: maximum deviation = 1e-1." Adding the 
% contents of each TP matrix element (e.g., TEST = sum(TP_04{1,1}(:,3)))
% suggests the presence of extremely small deviations from the desired
% figure (i.e., 228,357). The error does not appear fatal, in the sense
% that it persists in generating results. I had hoped I could eliminate the
% error by rounding the appropriate columns, but when I do this the model
% actually fails to run.

% for i = 1:24
% 
%     TP_04{i,1}(:,3) = round(TP_04{i, 1}(:,3),6);
%     RWD_04{i,1}(:,2) = round(RWD_04{i,1}(:,2),6);
%     
% end

%% Process TP & RWD Matrices for Incorporation Into MDPSolve

% Now we create sparse matrices around each of the TP and RWD matrices
% within the cell array (i.e., TP_List). The first and second arguments in
% the SPARSE function indicate what row and column each data point should
% occupy within each sparse matrix, while the last is the information we're
% actually trying to store. Surprisingly, this is not a computationally
% taxing operation, taking only seconds to perform.
for i = 1:24

    TP_List{i, 1} = ...
        sparse(TP_04{i,1}(:,1), TP_04{i,1}(:,2), TP_04{i,1}(:,3));

end

for i = 1:24

    % https://www.mathworks.com/help/matlab/ref/sortrows.html
    RWD_04{i,1} = sortrows(RWD_04{i,1}(:,1:2));
    RWD_List{i, 1} = RWD_04{i,1}(:,2);

end

% https://www.mathworks.com/help/matlab/ref/vertcat.html
TP = vertcat(TP_List{:}); RWD = horzcat(RWD_List{:});

%% Solving the Model

% Now we solve the model, using the MDPSolve package: 
% Fackler, P. L. (2014). MDPSOLVE: MATLAB tools for solving Markov 
% Decision Problems (Version 2014a). Raleigh, NC. Retrieved from 
% https://sites.google.com/site/mdpsolve/

%   NOTE 092018: For this particular sensitivity analysis, the parameter
%   being varied is in the MATLAB script itself (i.e., the discount factor)
%   - this is NOT true of the other analyses we are undertaking. In
%   general, delta will be set at its default value of 0.95.
years = 70;
clear model 
model.R = RWD; 
% NOTE 042519: Because the order of the sensitivity analysis components was
% revised (i.e., to enable similar components to be grouped together), the
% value of Sens_Analysis_ID was updated.
if Sens_Analysis_ID == 3
    model.d = Disc_Factor(z);
else
    model.d = 0.95;
end
model.P = TP; model.T = years;

% This line of code carries out backwards iteration to identify the optimal
% policy and corresponding value function. It should take approximately 5
% seconds to run.
[v] = mdpsolve(model);

%% Optimal Policy

% The objective now is to interpret the optimal policy. MATLAB stores the 
% results generated by the backward induction algorithm in the structure
% "v", which include both the optimal policy (i.e., Ixopt) and the value
% function (i.e., v).
Opt_Pol_00 = v.Ixopt;
% Because MDPSOLVE expresses the optimal policy in terms of a combination
% of state and action, to determine what the individual did we must filter
% out the state by dividing by the size of the state space (i.e., 228,357).
% There is a problem, however: the size of the state space will not
% necessarily divide evenly into the ID, resulting in non-integer values
% which are not really meaningful. To address this we apply the FLOOR
% function, which rounds down to the nearest integer. Now there is another
% problem: IDs with values less than 228,357 will, under this approach,
% lead to an action identifier of 0. Given that action identifier ranges
% between 1 and 24, this is also not meaningful. We therefore complete the
% calculation by incrementing by one.
Opt_Pol_00(:,2) = floor(Opt_Pol_00(:,1)/228357) + 1;
% The remainder after division of the ID by the size of the state space
% (i.e., the modulus) is simply the state the individual occupies ...
Opt_Pol_00(:,3) = mod(Opt_Pol_00(:,1),228357);
% ... but the remainder after division of 228,357 by 228,357 is of course
% zero, so we set this value manually (an ID of 0 is meaningless and may
% cause problems later on).
Opt_Pol_00(228357,3) = 228357;

% So far we have been dealing primarily with matrices; however, for the
% purposes of further analysis in R, we'd like to be able to label the
% columns of the construct we have just created; converting the matrix to a
% table enables us to do just this (we note that RefActions and RefStates
% were imported as tables, so this is necessary to enable us to combine
% them with Opt_Pol, as we shall do presently).
% https://www.mathworks.com/matlabcentral/answers/179129-is-it-possible-to-rename-columns-in-a-table
Opt_Pol_01 = array2table(Opt_Pol_00);
% Label the columns of Opt_Pol (the columns of RefActions and RefStates are
% already labelled).
Opt_Pol_01.Properties.VariableNames = {'Policy','ID_Actions','ID_States'};
% Perform two left joins on Opt_Pol_01.
Opt_Pol_01 = outerjoin(Opt_Pol_01, RefActions, 'Type', 'Left');
Opt_Pol_01 = outerjoin(Opt_Pol_01, RefStates, 'Type', 'Left');

% NOTE 042519: This has been modified to account for the incorporation of
% two-sensitivity analyses in the context of examining the impact of price
% changes.
writetable(Opt_Pol_01, strcat(MAT_Path, ...
    "\Sensitivity_Analysis\", ...
    Sens_Analysis_Types(Sens_Analysis_ID), ...
    "\MATLAB_Output\Optimal_Policies\Case_", ...
    num2str(z), Sens_Anal_String, "\OptPolicy_", Case, ".dat")); 

%% Value Function

% Now to deal with the value function. This is simpler, because we do not
% need to "interpret" the value function in the same way as the optimal
% policy. We should however be aware that while the latter was the action
% an individual should take in each state, the former is the reward they
% would then receive.
Val_Fun_00 = v.v;
% This simply creates a sequence of integers to serve as an identifier,
% analogous to what we have for the optimal policy table.
Val_Fun_00(:,2) = 1:228357;
% Convert the matrix to a table.
Val_Fun_01 = array2table(Val_Fun_00);
% Rename the columns appropriately.
Val_Fun_01.Properties.VariableNames = {'Value_Fun', 'ID_States'};

% Val_Fun_01 = outerjoin(Val_Fun_01, RefActions, 'Type', 'Left');
Val_Fun_01 = outerjoin(Val_Fun_01, RefStates, 'Type', 'Left');

% NOTE 042519: This has been modified to account for the incorporation of
% two-sensitivity analyses in the context of examining the impact of price
% changes.
writetable(Val_Fun_01, strcat(MAT_Path, ...
    "\Sensitivity_Analysis\", ...
    Sens_Analysis_Types(Sens_Analysis_ID), ...
    "\MATLAB_Output\Value_Functions\Case_", ...
    num2str(z), Sens_Anal_String, "\ValFun_", Case, ".dat"));

clearvars -except a z BASE Ref_Path MAT_Path Sens_Analysis_Types ...
    Sens_Analysis_ID Person_Type Disc_Factor RefActions RefStates ...
    Table_1 Table_2
    
end
end