%Michael Angus 2021
%Verification plots. 
%Start with rank histograms.
%follows on from load_all.m and ver_rd.m both in verification.
EMOS_tps_drynan=EMOS_tps;
EMOS_tps_drynan(EMOS_tps==0)=NaN;
Fcst_tps_drynan=Fcst_tps;
Fcst_tps_drynan(Fcst_tps==0)=NaN;
UQM_tps_drynan=UQM_tps;
UQM_tps_drynan(UQM_tps==0)=NaN;
IMD_tps_drynan=IMD_tps;
IMD_tps_drynan(IMD_tps==0)=NaN;
IMD_tps_all_drynan=IMD_tps_all;
IMD_tps_all_drynan(IMD_tps_all==0)=NaN;
IMERG_tps_drynan=IMERG_tps;
IMERG_tps_drynan(IMERG_tps==0)=NaN;
%%
%For the rankings, need to find a cell structure of actual rainfall values
%at each location. Then from there, rank with respect to those values. 
%also strip out if any NaN values.
for loc_rainvals=1:9
%forecast
A=squeeze(Fcst_tps(loc_rainvals,:,IMD_tps(loc_rainvals,:)>0));
A(24,:)=IMD_tps(loc_rainvals,IMD_tps(loc_rainvals,:)>0);
B=A(:,~isnan(mean(A)));
ranked_IMD_Fcst_tps{loc_rainvals}=sort(B,1,'descend');
ranked_IMD_Fcst_rain{loc_rainvals}=B(24,:);
%UQM
A2=squeeze(UQM_tps(loc_rainvals,:,IMD_tps(loc_rainvals,:)>0));
A2(24,:)=IMD_tps(loc_rainvals,IMD_tps(loc_rainvals,:)>0);
B2=A2(:,~isnan(mean(A2)));
ranked_IMD_UQM_tps{loc_rainvals}=sort(B2,1,'descend');
ranked_IMD_UQM_rain{loc_rainvals}=B2(24,:);
%EMOS
A3=squeeze(EMOS_tps(loc_rainvals,:,IMD_tps(loc_rainvals,:)>0));
A3(24,:)=IMD_tps(loc_rainvals,IMD_tps(loc_rainvals,:)>0);
B3=A3(:,~isnan(mean(A3)));
ranked_IMD_EMOS_tps{loc_rainvals}=sort(B3,1,'descend');
ranked_IMD_EMOS_rain{loc_rainvals}=B3(24,:);
end
%%
for stat=1:9
    for time=1:size(ranked_IMD_Fcst_tps{stat},2)
        fr_ens_mems=squeeze(ranked_IMD_Fcst_tps{stat}(:,time));
        UQM_ens_mems=squeeze(ranked_IMD_UQM_tps{stat}(:,time));
        EMOS_ens_mems=squeeze(ranked_IMD_EMOS_tps{stat}(:,time));
        try
        Ens_ranker_fcst{stat}(time)=find(fr_ens_mems==ranked_IMD_Fcst_rain{stat}(time));
        catch %theres some draws, so just randomize between the equals.
        X=find(fr_ens_mems==ranked_IMD_Fcst_rain{stat}(time));
        Ens_ranker_fcst{stat}(time)=randi([X(1),X(end)]);
        end
        try
        Ens_ranker_UQM{stat}(time)=find(UQM_ens_mems==ranked_IMD_UQM_rain{stat}(time));
        catch %theres some draws, so just randomize between the equals.
        Y=find(UQM_ens_mems==ranked_IMD_UQM_rain{stat}(time));
        Ens_ranker_UQM{stat}(time)=randi([Y(1),Y(end)]);
        end
        try
        Ens_ranker_EMOS{stat}(time)=find(EMOS_ens_mems==ranked_IMD_EMOS_rain{stat}(time));
        catch %theres some draws, so just randomize between the equals.
        Z=find(EMOS_ens_mems==ranked_IMD_EMOS_rain{stat}(time));
        Ens_ranker_EMOS{stat}(time)=randi([Z(1),Z(end)]);
        end
    end
end
%% make the plot
loc_opts={'Mumbai','Rajasthan','Kerala','Shimla','Delhi',...
    'Hyderabad','Patna','Bhubaneswar','Meghalaya'};
figure('Renderer', 'painters', 'Position', [40 40 1200 600])
Orig_Ens_count=zeros(9,24);
UQM_Ens_count=zeros(9,24);
EMOS_Ens_count=zeros(9,24);

for locs=1:9
    for hist_place=1:24
    Orig_Ens_count(locs,hist_place)=length(find(Ens_ranker_fcst{locs}==hist_place))/length(Ens_ranker_fcst{locs});
    UQM_Ens_count(locs,hist_place)=length(find(Ens_ranker_UQM{locs}==hist_place))/length(Ens_ranker_UQM{locs});
    EMOS_Ens_count(locs,hist_place)=length(find(Ens_ranker_EMOS{locs}==hist_place))/length(Ens_ranker_EMOS{locs});
    end
    bar_plot_comb=cat(2,Orig_Ens_count(locs,:)',UQM_Ens_count(locs,:)',EMOS_Ens_count(locs,:)');
    subplot(3,3,locs)
    h_bar=bar(bar_plot_comb);
    h_bar(1).EdgeColor='none';
    h_bar(2).EdgeColor='none';
    h_bar(3).EdgeColor='none';
    h_bar(1).FaceColor=[0 66 37]/255;
    h_bar(2).FaceColor=[242 9 242]/255;
    h_bar(3).FaceColor=[0 204 204]/255;
    grid on
    if ismember(locs,[1 3 4])
    ylim([0 .7])
    else
    ylim([0 .5])
    end
    title(loc_opts(locs));
    set(gca,'xtick',0:2:24);
    set(gca,'ytick',0:.1:.7);
    if locs==1
    axP = get(gca,'Position');
    legend('Original Ensemble','QM ensemble','EMOS ensemble',...
        'location','northwestoutside')
    legend boxoff   
    set(gca, 'Position', axP)
    end
end
    % add some general axes
    ht = text(-70,0.65,...
        'Frequency of days');
    set(ht,'Rotation',90)
    set(ht,'FontSize',14)
    ht2=text(-33,-0.2,'Rank of IMD to ensemble members');
    set(ht2,'Fontsize',14)
saveas(gcf,'rank_hist_IMD_EMOS_DGQM','png');