% DP Solver Utility (Created 05/05/18; Last Modified 04/21/19)

%% Initialization

% NOTE 042319: Updated script to support extension of sensitivity analysis.
% 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).
% NOTE 011219: Added new alternative strategy in which people non-adherent
% to LTPA and healthy eating habits align their health-related behaviours
% with the optimal policy following diagnosis with T2DM.
% NOTE 031719: Compiled reward vectors corresponding to the optimal policy
% and the three alternate policies.

clearvars; clc;

% https://stackoverflow.com/questions/37830233/matlab-path-to-the-current-script
% https://stackoverflow.com/questions/2720140/find-location-of-current-m-file-in-matlab
% https://stackoverflow.com/questions/6091064/matlab-is-there-a-way-to-get-the-path-of-the-current-script
% https://www.mathworks.com/help/matlab/ref/strrep.html
% https://stackoverflow.com/questions/50151667/matlab-remove-everything-after-the-final-blank-in-a-character-array
% https://stackoverflow.com/questions/830855/what-regex-would-capture-everything-from-mark-to-the-end-of-a-line
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\");
R_alpha = 1:2; R_beta = 1:3; R_gamma = 1:4;

% 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\"))

% tic
for alpha = R_alpha
    for beta = R_beta
        % https://stackoverflow.com/questions/5080191/loop-through-certain-integers-and-in-the-order-that-i-specify
        for gamma = R_gamma
%% Import data from R

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

% 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'))'));

% 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);
RWD_03 = arrayfun(@(x) ...
    strcat(MAT_Path, "\Baseline_Results\RWD_Matrix_Components\Case_", ...
    Case, "\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);

%% 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/
delta = 0.95; years = 70;
clear model 
model.R = RWD; model.d = delta; 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 just under 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');

writetable(Opt_Pol_01, strcat(MAT_Path, ...
    "\Baseline_Results\MATLAB_Output\Optimal_Policies\Case_", Case, ...
    "\", "OptPolicy_", Case, ".dat"));

% 101318: Calculating the pstar matrix requires selecting the rows of the
% transition probability matrix that correspond to the state/action
% combination the individual chooses in each state they could occupy. There
% are a number of ways to do this. The approach used in previous versions
% of this script, which is presented in the scrapyard, appears to have been
% incredibly inefficient, requiring nearly half an hour to complete, and
% nearly a full hour if the same approach was applied to both the optimal
% transition probability matrix and the alternate transition probability
% matrix. By contrast, the script below runs in less than a fifth of a
% second. The lesson, once again, is not to use a for-loop unless
% absolutely necessary.
pstar_opt = TP(Opt_Pol_01.Policy, :);
% NOTE 031719: We now endeavour to isolate the rewards associated with the
% pursuit of the optimal policy and conceivable alternate policies. In the
% first step, we stack the reward vector associated with the 24 possible
% courses of action available to the decision-maker, one atop the other.
% This is directly analogous to how we prepare the transition probability
% matrix prior to applying backward induction to solve the dynamic
% optimization problem. This implies the result will have 228357 * 24 = 
% 5480568 rows.
rwdstar = vertcat(RWD_List{:}); 
% Now we pick out the 228357 rows of the construct we just created, that
% correspond to the choice of actions the individual should take in each
% starting state (i.e., the optimal policy). This, again, is directly
% analogous to the derivation of the pstar matrix.
rwdstar_opt = rwdstar(Opt_Pol_01.Policy, :);
% We save the optimal and alternate probability transition matrices for
% later - we'll need this when we implement simulations using our results.
% https://www.mathworks.com/help/matlab/ref/save.html
save(strcat(MAT_Path, ...
    "Simulation_Results\", "pstar_opt_", Case, ".mat"), 'pstar_opt')
save(strcat(MAT_Path, ...
    "Simulation_Results\", "rwdstar_opt_", Case, ".mat"), 'rwdstar_opt')

% In preparation for the simulation analysis, we now consider what the
% policy would look like for other combinations of health behaviours. Begin
% by reading in the optimal policy. Note the inclusion of the parameter
% 'TextType', which is intended to ensure variables in the alternative
% policy table will be strings rather than cell arrays. Incidentally ...
% 1) A similar operation can be performed by applying the function STRING. 
% 2) A variable with STRING data type can always be converted into a cell
% array of strings using the function CELLSTR.
Alt_Pol = readtable(strcat(MAT_Path, ...
    "\Baseline_Results\MATLAB_Output\Optimal_Policies\Case_", Case, ...
    "\", "OptPolicy_", Case, ".dat"), 'TextType', 'String');
% 011119: Expand Alt_Pol_Array to contain one new scenario in which people
% who previously did not adhere to either LTPA or healthy eating habits
% become adherent upon diagnosis with T2DM.
Alt_Pol_Array = cell(9,1); pstar_alt_array = cell(9,1);
rwdstar_alt_array = cell(9,1);

% Define what health behaviour characterizes the alternative policy. Since
% there are eight possible combinations, Des_Health_Behav is in [1, 8].
% 011119: We are introducing one new scenario, but will compile this
% outside the loop, since it is derived in a slightly different way.
for DSB = 1:9
    
    Alt_Pol_Temp = Alt_Pol;
    % Duration is defined as a character. Since we'll be carrying out
    % mathematical operations with this variable, we'll need to convert it to
    % numeric data type.
    Alt_Pol_Temp.Duration_1 = str2double(Alt_Pol_Temp.Duration_1);
    % We now want to define each possible level of labour market participation
    % as an integer. In R, this can be easily accomplished using CASE_WHEN.
    % However, I could not find evidence that a similar function exists for
    % MATLAB, and acheiving this using a conditional statement (i.e., IF)
    % turned out to be difficult. I settled upon an alternative approach that
    % calculates the weighted sum of the product of a logical statement and the
    % desired corresponding integer value. Since only one of the logical
    % statements can be true at any time, and since a value of 1 is assigned to
    % conditions that are determined to be true, this is effectively equivalent
    % to conditional assignment of integer values based on the value in
    % Alt_Pol_01.LMP_01.
    Alt_Pol_Temp.LMP_01 = (Alt_Pol_Temp.LMP == "FT").*1 + ...
        (Alt_Pol_Temp.LMP == "PT").*2 + (Alt_Pol_Temp.LMP == "RT").*3;
    % This is a complicated statement. What we are trying to convey is that the
    % alternate policy should be the same as the optimal policy except in a
    % relatively small number of cases. This is because health-related choices
    % for individuals with T2DM are effectively fixed, except at specific
    % times (i.e., Years 0 and 9). Individuals who are healthy are free to make
    % whatever health-related choices they like. Therefore, if the person is
    % healthy OR has been newly-diagnosed with T2DM OR has had T2DM for 9
    % years, then apply the adjustment. Otherwise, leave things as they are.
    % 
    % How does this work? As in the last line of code, we are calculating the
    % sum of two products. Both products include logical statements, only one
    % of which can be true at any given time. If the first statement is true,
    % the alternative policy is captured in the value of an integer that
    % includes:
    % 1) An integer that reflects the individual's current state.
    % 2) A second integer that effectively acts as a "bookmark", narrowing down
    % the range of values the alternative policy can taken by assuming the
    % individual's chosen level of labour market participation would remain
    % unchanged (noting that the first 1826856 values are associated with
    % individuals working full-time, the second 1826856 values are associated
    % with individuals working part-time, etc.).
    % 3) A third integer that simply indicates what health-related behaviour
    % the individual should adopt under the alternative policy.
    
    if DSB ==9 % 011119: Compilation of new alternate policy begins here.
        % 011119: This is the only line that actually deviates from the
        % above description. What we are saying is this:
        % a) If you are healthy, the integer denoting the health-related 
        % behaviour you should adopt is set to 8, which denotes uniform 
        % non-adherence to health-related behaviours.
        % b) Otherwise, do what you would have done anyway (i.e., follow 
        % the optimal policy).
        %         
        % The intent here is to address the question of whether it is
        % possible to significantly reduce indirect costs associated with
        % T2DM after having been diagnosed. This is important, because it
        % speaks to the issue of whether it is more efficient to
        % concentrate our health care resources on T2DM prevention or
        % management.
        Alt_Pol_Temp.Alt_Policy = ...
            (Alt_Pol_Temp.Health_Stat_1 == "Healthy").*...
            (Alt_Pol_Temp.ID_States_Opt_Pol_01 + ...
            (228357*((Alt_Pol_Temp.LMP_01 - 1).*8 + 8 - 1))) ...
            + (1-(Alt_Pol_Temp.Health_Stat_1 == "Healthy")).*...
            Alt_Pol_Temp.Policy;
    else
        Alt_Pol_Temp.Alt_Policy = ...
            (Alt_Pol_Temp.Health_Stat_1 == "Healthy" | ...
            ismember(Alt_Pol_Temp.Duration_1, [0, 9])).*...
            (Alt_Pol_Temp.ID_States_Opt_Pol_01 + ...
            (228357*((Alt_Pol_Temp.LMP_01 - 1).*8 + DSB - 1))) ...
            + (1-(Alt_Pol_Temp.Health_Stat_1 == "Healthy" |...
            ismember(Alt_Pol_Temp.Duration_1, [0, 9]))).*...
            Alt_Pol_Temp.Policy; 
    end

    % Alt_Pol_01.Alt_Policy reflects both the state the individual 
    % currently occupies as well as the action they've chosen to take. We 
    % want to convert this to an integer that reflects only the latter.
    Alt_Pol_Temp.ID_Actions_Alt_Pol_01 = ...
        floor(Alt_Pol_Temp.Alt_Policy/228357) + 1;

    % NOTE 112318: We now save the alternate policy. We can call upon this 
    % later when we need to track employment status within the simulation
    % cohort; however, we do not really need to, since the optimal and
    % alternate policies differ only by health-related choices, not labour
    % decisions around labour-market participation. Since we are storing 
    % this information solely to exclude retirement periods from the
    % efficiency loss calculations, nothing is loss by simply referencing 
    % the optimal policy.
    %     writetable(Alt_Pol_Temp, strcat(MAT_Path, "\Simulation_Results\", ...
    %     "AltPolicy_", Case, "_", "hybrid", ".dat"));

    pstar_alt_temp = TP(Alt_Pol_Temp.Alt_Policy, :);
    rwdstar_alt_temp = rwdstar(Alt_Pol_Temp.Alt_Policy, :);
    
    if DSB ==9
        save(strcat(MAT_Path, "Simulation_Results\", "pstar_alt_", ...
            Case, "_", "hybrid", ".mat"), 'pstar_alt_temp')
        save(strcat(MAT_Path, "Simulation_Results\", "rwdstar_alt_", ...
            Case, "_", "hybrid", ".mat"), 'rwdstar_alt_temp')
    else % 011119: ENDS HERE
        save(strcat(MAT_Path, "Simulation_Results\", "pstar_alt_", ...
            Case, "_", string(DSB), ".mat"), 'pstar_alt_temp')
        save(strcat(MAT_Path, "Simulation_Results\", "rwdstar_alt_", ...
            Case, "_", string(DSB), ".mat"), 'rwdstar_alt_temp')
    end
            
    Alt_Pol_Array{DSB} = Alt_Pol_Temp; 
    pstar_alt_array{DSB} = pstar_alt_temp;
    rwdstar_alt_array{DSB} = rwdstar_alt_temp;
    
    clearvars Alt_Pol_Temp pstar_alt_temp rwdstar_alt_temp
    
end

% NOTE 092518: However, we'll also need to compile a matrices that sum
% probabilities in a cumulative fashion. This is because our optimal
% transition probability matrices are sparse, whereas the original MDPSIM
% function was developed with ordinary (i.e., dense) functions in mind.

% First, we identify all the non-zero elements in the sparse matrix,
% splitting its contents into three components - namely, the row and column
% indices and the values themselves.
% https://www.mathworks.com/matlabcentral/answers/338433-how-to-save-variables-from-function-output-to-workspace
% 110118: This line of code results in the following error, beginning with
% the second loop: "Unable to perform assignment because the indices on the
% left side are not compatible with the size of the right side. I believe
% the source of the error lies in the fact that by the second loop,
% pstar_opt_cum will have been defined as a 228357x228357 sparse double
% (due to additional manipulations carried out below), even while the RHS
% of the expression (i.e., find(pstar_opt)) is a 557360x1 double. This
% issue can be addressed by clearing the value of pstar_opt_cum on every
% loop (which I had actually done in a previous version of this file but
% then relegated to a comment before the first "end" statement).
[pstar_opt_cum(:,1), pstar_opt_cum(:,2), pstar_opt_cum(:,3)] = ...
    find(pstar_opt);
% Next, we need to sort on the first two columns, so that all of the
% starting state values are clumped together. It's noteworthy that there's
% no need to insert any brackets, which speaks to the efficiency of the
% MATLAB syntax.
% https://www.mathworks.com/help/matlab/ref/sortrows.html
pstar_opt_cum = sortrows(pstar_opt_cum,1:2);
% Now we replace the old column of probabilities by a new one, which has
% applied a cumulative summation operation. Note that we must apply the 
% cell2mat function to avoid output being aggregated across columns (it is
% hard to articulate what occurs since I do not fully understand it - but I
% do know the output that is otherwise generated is not suitable for our
% purposes).
pstar_opt_cum2 = [pstar_opt_cum(:,1), pstar_opt_cum(:,2), ...
    cell2mat(accumarray(pstar_opt_cum(:,1), pstar_opt_cum(:,3),[],...
    @(x){cumsum(x)}))];
% Convert to a table so that we can name the columns. This, in turn, will
% enable us to merge multiple data structures.
pstar_opt_cum2 = array2table(pstar_opt_cum2);
pstar_opt_cum2.Properties.VariableNames = {'Start','End','CumProb'};

% Now we contend with a second problem; the maximum value in each row is
% not necessarily 1, due (presumably) to small rounding errors at some
% point in the calculations. To address this, we first identify the largest
% value in every group, where groups refer in this context to starting
% states.
Groups = findgroups(pstar_opt_cum2(:,1));
% The next step is to apply the function, which involves identifying the
% largest value in the third column of PSTAR_CUM2, by group.
% https://www.mathworks.com/help/matlab/ref/splitapply.html
MaxVals = splitapply(@max, pstar_opt_cum2(:,3), Groups);
% Create an index we can subsequently match on when conducting the merge.
MaxVals = [transpose(1:228357), MaxVals];
% Convert MAXVALS to a table (again, it appears we cannot merge simple
% arrays).
MaxVals = array2table(MaxVals);
% Name the columns appropriately. Note that the column 'Start' appears in
% both PSTAR_CUM2 and MAXVALS - which is at it should be, since this is the
% variable we will be merging on.
MaxVals.Properties.VariableNames = {'Start','MaxProb'};
% Now carry out the merging operation.
pstar_opt_cum3 = outerjoin(pstar_opt_cum2, MaxVals, 'Type', 'Left');
% Convert the resulting table into an array, as the command we will be
% executing next doesn't seem to work well with tables.
pstar_opt_cum3 = table2array(pstar_opt_cum3);
% This is an important line of code. Essentially, we compare our original
% probabilities against the maximum probabilities associated with each
% group. If the values match, the cumulative probability should be one or
% close to one; otherwise, it is less than one. We are telling MATLAB to
% assign 1 to each of these cells, thereby addressing cases where the value
% would otherwise be slightly more or less than one.
% https://stackoverflow.com/questions/18741764/matlab-vectorization-within-if-else-if-else-statements
pstar_opt_cum3(pstar_opt_cum3(:,3)==pstar_opt_cum3(:,5),3) = 1;
% There's another significant feature of the line of code above, which is
% how quickly it runs, taking advantage of vectorization.
% The alternative would be to use a for-loop such as the one below - but
% this is MUCH slower and thus less efficient.
% for i = 1:571279
% 
%         if pstar_cum3.CumProb(i)==pstar_cum3.MaxProb(i)
%             pstar_cum3.CumProb(i) = 1;
%         else 
%         end
%     
% end

% Convert the matrix back into a sparse matrix, then save the result for
% use in our simulations.
pstar_opt_cum = sparse(pstar_opt_cum3(:,1), pstar_opt_cum3(:,2), ...
    pstar_opt_cum3(:,3));
save(strcat(MAT_Path, ...
    "Simulation_Results\", "pstar_opt_cum_", Case, ".mat"), ...
    'pstar_opt_cum')

for DSB = 1:9
    
    % Now do exactly the same thing with the alternative policy.
    [pstar_alt_cum_temp(:,1), pstar_alt_cum_temp(:,2), ...
        pstar_alt_cum_temp(:,3)] = find(pstar_alt_array{DSB});
    pstar_alt_cum_temp = sortrows(pstar_alt_cum_temp,1:2);

    pstar_alt_cum_temp2 = [pstar_alt_cum_temp(:,1), ...
        pstar_alt_cum_temp(:,2), ...
        cell2mat(accumarray(pstar_alt_cum_temp(:,1), ...
        pstar_alt_cum_temp(:,3),[],...
        @(x){cumsum(x)}))];
    pstar_alt_cum_temp2 = array2table(pstar_alt_cum_temp2);
    pstar_alt_cum_temp2.Properties.VariableNames = ...
        {'Start','End','CumProb'};

    Groups = findgroups(pstar_alt_cum_temp2(:,1));
    MaxVals = splitapply(@max, pstar_alt_cum_temp2(:,3), Groups);
    MaxVals = [transpose(1:228357), MaxVals];
    MaxVals = array2table(MaxVals);
    MaxVals.Properties.VariableNames = {'Start','MaxProb'};

    pstar_alt_cum_temp3 = ...
        outerjoin(pstar_alt_cum_temp2, MaxVals, 'Type', 'Left');
    pstar_alt_cum_temp3 = table2array(pstar_alt_cum_temp3);
    pstar_alt_cum_temp3(pstar_alt_cum_temp3(:,3) == ...
        pstar_alt_cum_temp3(:,5),3) = 1;
    pstar_alt_cum_temp = sparse(pstar_alt_cum_temp3(:,1), ...
        pstar_alt_cum_temp3(:,2), pstar_alt_cum_temp3(:,3));

    % 011119: Compilation of new alternate policy begins here.
    if DSB == 9
        save(strcat(MAT_Path, ...
            "Simulation_Results\", "pstar_alt_cum_", Case, "_", ...
            "hybrid", ".mat"), 'pstar_alt_cum_temp')
    else % 011119: ENDS HERE
        save(strcat(MAT_Path, ...
            "Simulation_Results\", "pstar_alt_cum_", Case, "_", ...
            string(DSB), ".mat"), 'pstar_alt_cum_temp')      
    end
    
    clearvars pstar_alt_cum_temp pstar_alt_cum_temp2 ...
        pstar_alt_cum_temp3 Groups MaxVals
    
end

%% 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');

writetable(Val_Fun_01, strcat(MAT_Path, ...
    "\Baseline_Results\MATLAB_Output\Value_Functions\Case_", Case, "\", ...
    "ValFun_", Case, ".dat"));

clearvars -except alpha beta gamma R_alpha R_beta R_gamma ...
    BASE Ref_Path MAT_Path RefActions RefStates

        end
    end
end
% toc