% This script loads and plots score dependent policy experiments and generates Figures 3 and 4

%% Setup

experi.name = 'Debt-to-Earnings Limit (in %)';
experi.short = 'di';

edulong = {'no college','college'};

% Output variable names and rows in output file
outPlot = {'bankruptcies','debt-to-income','interest','welfare'}';
outPlot(:,2) = num2cell([5,13,10,54]);

% Initialize Figures
Figure3 = figure('Visible', 'off');
Figure3col = figure('Visible', 'off');
Figure4 = figure('Visible', 'off');
plotOrder = [1 2 5 6 3 4 7 8];  % order plots for Figure 4
Figure4.Position = [100, 100, 1200, 600]; % Make Figure 4 wider

%% Import data from text file.

%Get names of all subfolders that contain experiments

ee = 0;
for edugroup = {'hs','col'}
    ee = ee+1;


    % foldercount = sum(vertcat(files.isdir)==1);
    policyFolders = {};
    for i = 1:length(experimentFolders)
        if  contains(experimentFolders(i).name,experi.short) && contains(experimentFolders(i).name,'_score') && contains(experimentFolders(i).name,edugroup) && ~contains(experimentFolders(i).name,'robust')
            policyFolders{end+1,1} = fullfile(runDir,experimentFolders(i).name);
        end

        if contains(experimentFolders(i).name,'benchmark') && contains(experimentFolders(i).name,edugroup)  && ~contains(experimentFolders(i).name,'dump')
            benchmark = load(fullfile(runDir,experimentFolders(i).name,'NewOrderPlus.txt'));
        end
    end

    % generate input-output matrix
    neworders = zeros(length(benchmark),length(policyFolders));
    experiments = cell(3,length(policyFolders));

    for ii = 1:length(policyFolders)
        % retrieve outcomes
        neworders(:,ii) = load([policyFolders{ii},'/NewOrderPlus.txt']);

        % retrieve experiment details from folder name.
        % if name contains di_11_s22: d/i limit = 11, score limit = 22
        % type of experiment
        experiments(1,ii) = {experi.short};
        % policy cap
        cap_pat = experi.short + "_" + digitsPattern + "_"; %where does cap sit?
        experiments{2,ii} = str2double(extract(extract(policyFolders{ii},cap_pat),digitsPattern));
        % type score threshold
        experiments{3,ii} = str2double(extract(extract(policyFolders{ii},"_s" + digitsPattern),digitsPattern));
    end

    %% PREPROCESSING

    % devide type scores / 100 to be in [0,1]
    experiments(3,:) = num2cell([experiments{3,:}]./100);

    % sort experiments AND OUTCOMES such that policy caps and type scores are ascending
    [experiments, reorderindex] = sortrows(experiments',[1 2 3]);
    experiments = experiments';
    neworders = neworders(:,reorderindex);

    % calculate annual D/I caps
    experiments(2,strcmp(experiments(1,:),experi.short)) = num2cell([experiments{2,strcmp(experiments(1,:),experi.short)}].*3);

    % calculate CEV with new outputs but place at old position
    welfareID = outPlot{strcmp(outPlot,'welfare'),2};
    neworders(welfareID:welfareID+1,:) = (benchmark(welfareID:welfareID+1)./neworders(welfareID:welfareID+1,:) - 1).*100;

    % calculate filing annually
    filingID = outPlot{strcmp(outPlot,'bankruptcies'),2};
    neworders(filingID:filingID+1,:) = neworders(filingID:filingID+1,:)./3.*100;

    benchmark(filingID:filingID+1,:) = benchmark(filingID:filingID+1,:)./3.*100;

    % calculate interest rates annually
    interestID = outPlot{strcmp(outPlot,'interest'),2};
    neworders(interestID:interestID+1,:) = ((1+neworders(interestID:interestID+1,:)).^(1/3)-1).*100;
    benchmark(interestID:interestID+1,:) = ((1+benchmark(interestID:interestID+1,:)).^(1/3)-1).*100;

    % calculate debt-to-income annually
    dtiID = outPlot{strcmp(outPlot,'debt-to-income'),2};
    neworders(dtiID:dtiID+1,:) = neworders(dtiID:dtiID+1,:).*3.*100;
    benchmark(dtiID:dtiID+1,:) = benchmark(dtiID:dtiID+1,:).*3.*100;

    %% PLOT POLICY EXPERIMENTS

    for j=1:length(outPlot)

        % Plot D/I limit change

        % which type score cut-off was used?
        modecap = mode([experiments{3,strcmp(experiments(1,:),experi.short)}]);
        colplot = strcmp(experiments(1,:),experi.short) & [experiments{3,:}]==modecap; %pick caps with type score cutoff used the most

        % define variable to be plotted
        plotvariable = [neworders(outPlot{j,2},colplot) ; %rational
            neworders(outPlot{j,2}+1,colplot)];  %behavioral

        if strcmp(outPlot{j,1},'welfare')
            plotName = [outPlot{j,1},' (CEV)'];
            xbench=[];
            ybenchreal=[];
            ybenchgoofy=[];
        else
            plotName = outPlot{j,1};
            xbench=linspace(min([experiments{2,colplot}]),max([experiments{2,colplot}]),10);
            ybenchreal = linspace(benchmark(outPlot{j,2}),benchmark(outPlot{j,2}),10);
            ybenchgoofy = linspace(benchmark(outPlot{j,2}+1),benchmark(outPlot{j,2}+1),10);
        end
        if strcmp(edugroup,'hs')
            figure(Figure3)
        else
            figure(Figure3col)
        end
        subplot(2,2,j)
        plot([experiments{2,colplot}],plotvariable(1,:),'Color',coloruse(1,:), 'LineWidth', 2)
        hold on
        plot(xbench,ybenchreal,'Color',coloruse(1,:),'LineStyle',':', 'LineWidth', 2);
        plot([experiments{2,colplot}],plotvariable(2,:),'Color',coloruse(2,:),'LineStyle','--', 'LineWidth', 2);
        plot(xbench,ybenchgoofy,'Color',coloruse(2,:),'LineStyle','-.', 'LineWidth', 2);
        if j == 1
            legend('Rational','BM_{Rational}','Behavioral','BM_{Behavioral}','Location','E')
            %             legend boxoff
        end

        % if j >= 3 %only xaxis when at bottom in paper
        xlabel(experi.name)
        % end
        ylabel('in %')
        title([plotName,' -- ',edulong{ee}])



        % Plot type score threshold change

        colplot = [experiments{2,:}] == mode([experiments{2,strcmp(experiments(1,:),experi.short)}])...
            & strcmp(experiments(1,:),experi.short); %pick most common d/i cap

        % define variable to be plotted
        plotvariable = [neworders(outPlot{j,2},colplot) ; %rational
            neworders(outPlot{j,2}+1.*(outPlot{j,2}~=8 && outPlot{j,2}~=16),colplot)];  %behavioral


        if strcmp(outPlot{j,1},'welfare')
            plotName = [outPlot{j,1},' (CEV)'];
            xbench=[];
            ybenchreal=[];
            ybenchgoofy=[];
        else
            plotName = outPlot{j,1};
            xbench=[min([experiments{3,colplot}]),max([experiments{3,colplot}])];
            ybenchreal = [benchmark(outPlot{j,2}),benchmark(outPlot{j,2})];
            ybenchgoofy = [benchmark(outPlot{j,2}+1),benchmark(outPlot{j,2}+1)];
        end

        figure(Figure4)
        subplot(2,4,plotOrder((ee-1)*length(outPlot)+j))
        plot([experiments{3,colplot}],plotvariable(1,:),'Color',coloruse(1,:), 'LineWidth', 2)
        hold on
        plot(xbench,ybenchreal,'Color',coloruse(1,:),'LineStyle',':', 'LineWidth', 2);
        plot([experiments{3,colplot}],plotvariable(2,:),'Color',coloruse(2,:),'LineStyle','--', 'LineWidth', 2);
        plot(xbench,ybenchgoofy,'Color',coloruse(2,:),'LineStyle','-.', 'LineWidth', 2);
        if j == 1 && ee == 1
            legend('Rational','BM_{Rational}','Behavioral','BM_{Behavioral}','Location','W')
            %             legend boxoff
        end

        % if j >= 3 %only xaxis when at bottom in paper
        xlabel('Type Score Threshold')
        % end

        ylabel('(in %)')
        title([plotName,' -- ',edulong{ee}])

        if false
            print('-dpng',[outFolder,matlab.lang.makeValidName(extractBefore(experi.name,"(in %)")),...
                '_change_score_',edulong{ee},'_',outName{j}])
        end
    end
end

saveas(Figure3,fullfile(outFolder,'Figure3.png'))
saveas(Figure4,fullfile(outFolder,'Figure4.png'))
close all

disp('Success!')


