% Copyright 2020 Enrico Camporeale
% Licensed under a FreeBSD License (see below)

% evaluate the score of classification models
% on a given time period

clear all,
close all

station_list={'frn','ott','iqa'};
%station_list={'iqa'};

%root_dir = cd;
%models_dir=[root_dir '/models/'];
%data_dir=[root_dir '/data/'];

station_number=0;

pwd_dir=pwd;

fig1=figure('Position',[1250 100 440 855]);
fig2=figure('Position',[1250 100 440 855]);
fig3=figure('Position',[1250 100 440 855]);
fig4=figure('Position',[1250 100 440 855]);

fig5=figure('Position',[390 530 1000 420]);


root = ['~/geospace_work/'];

markers={'o','+','x'};
for station = station_list
 station_number = station_number +1;
 cd([root char(station)])
 disp(['Processing ' char(station)])
 
 load([char(station) '_thresholds.mat'])

 filename = [char(station) '_train_test_sets.mat']
 load(filename)
 
% no geospace  
  load([char(station) '_trained_model_classification_black_box.mat'])

  Mdl_nogeospace=Mdl;

  X=table2array(T_test(:,4:end-1));
  
%  X=X*C;
  
 for class=1:5
  t = double(table2array(T_test(:,end))/60>=threshold_list(class));
 
  ypred_val=(predict(Mdl_nogeospace{class},X));
  pod_val(class) = sum(round(ypred_val)==1 & t==1)/sum(t==1);
  pofd_val(class) = 1-sum(round(ypred_val)==0 & t==0)/sum(t==0);
  tss(class) = pod_val(class)-pofd_val(class);

  TP = sum(round(ypred_val)==1 & t==1);
  TN = sum(round(ypred_val)==0 & t==0);
  FP = sum(round(ypred_val)==1 & t==0);
  FN = sum(round(ypred_val)==0 & t==1);
    
  hss(class) = 2*(TP.*TN - FN.*FP)./((TP+FN).*(FN+TN) + (TP+FP).*(FP+TN));
  
   
  end


  figure(fig1)
  subplot(3,1,station_number)
  plot(threshold_list,pod_val,'-o','LineWidth',2);  ylabel('Prob.  detection'),hold on

  figure(fig2)
  subplot(3,1,station_number)
  plot(threshold_list,pofd_val,'-o','LineWidth',2);  ylabel('Prob.  false detection'),hold on

  figure(fig3)
  subplot(3,1,station_number)
  plot(threshold_list,tss,'-o','LineWidth',2);  ylabel('True Skill Statistic'),hold on

  figure(fig4)
  subplot(3,1,station_number)
  plot(threshold_list,hss,'-o','LineWidth',2);  ylabel('Heidke Skill Score'), hold on

% surrogate  

  load([char(station) '_trained_model_classification_grey_box.mat'])

  Mdl_surr=Mdl;
  X=table2array(T_test(:,2:end-1));
  
%    X=X*C;

 
  
 for class=1:5
  t = double(table2array(T_test(:,end))/60>=threshold_list(class));
 
  ypred_val=(predict(Mdl_surr{class},X));
  pod_val(class) = sum(round(ypred_val)==1 & t==1)/sum(t==1);
  pofd_val(class) = 1-sum(round(ypred_val)==0 & t==0)/sum(t==0);
  tss(class) = pod_val(class)-pofd_val(class);

  TP = sum(round(ypred_val)==1 & t==1);
  TN = sum(round(ypred_val)==0 & t==0);
  FP = sum(round(ypred_val)==1 & t==0);
  FN = sum(round(ypred_val)==0 & t==1);
    
  hss(class) = 2*(TP.*TN - FN.*FP)./((TP+FN).*(FN+TN) + (TP+FP).*(FP+TN));

  end

  tss_gb=tss;hss_gb=hss;

  figure(fig1)
  subplot(3,1,station_number)
  plot(threshold_list,pod_val,'-o','LineWidth',2);  ylabel('Prob.  detection'),hold on

  figure(fig2)
  subplot(3,1,station_number)
  plot(threshold_list,pofd_val,'-o','LineWidth',2);  ylabel('Prob.  false detection'),hold on

  figure(fig3)
  subplot(3,1,station_number)
  plot(threshold_list,tss,'-o','LineWidth',2);  ylabel('True Skill Statistic'),hold on

  figure(fig4)
  subplot(3,1,station_number)
  plot(threshold_list,hss,'-o','LineWidth',2);  ylabel('Heidke Skill Score'), hold on

% white-box

 for class=1:5
  t = double(table2array(T_test_geo(:,end))/60>=threshold_list(class));
 
  ypred_val=double((T_test_geo.("maxdBdt_geo"))/60>=threshold_list(class));
  pod_val(class) = sum(round(ypred_val)==1 & t==1)/sum(t==1);
  pofd_val(class) = 1-sum(round(ypred_val)==0 & t==0)/sum(t==0);
  tss(class) = pod_val(class)-pofd_val(class);

  TP = sum(round(ypred_val)==1 & t==1);
  TN = sum(round(ypred_val)==0 & t==0);
  FP = sum(round(ypred_val)==1 & t==0);
  FN = sum(round(ypred_val)==0 & t==1);
    
  hss(class) = 2*(TP.*TN - FN.*FP)./((TP+FN).*(FN+TN) + (TP+FP).*(FP+TN));

  
  end

   tss_wb=tss;hss_wb=hss;

    figure(fig1)
  subplot(3,1,station_number)
  plot(threshold_list,pod_val,'-o','LineWidth',2);  ylabel('Prob.  detection'),hold on

  title(char(upper(station)+" station"),'fontsize',12)
  set(gca,'fontsize',14)
  xlim([threshold_list(1)*0.9 threshold_list(end)*1.1])
  xlabel('Threshold (nT/s)')
  legend('black-box','grey-box','white-box');legend('boxoff')

  figure(fig2)
  subplot(3,1,station_number)
  plot(threshold_list,pofd_val,'-o','LineWidth',2);  ylabel('Prob.  false detection'),hold on

   title(char(upper(station)+" station"),'fontsize',12)
  set(gca,'fontsize',14)
  xlim([threshold_list(1)*0.9 threshold_list(end)*1.1])
xlabel('Threshold (nT/s)')
legend('black-box','grey-box','white-box');legend('boxoff')

  figure(fig3)
  subplot(3,1,station_number)
  plot(threshold_list,tss,'-o','LineWidth',2);  ylabel('True Skill Statistic'),hold on

   title(char(upper(station)+" station"),'fontsize',12)
  set(gca,'fontsize',14)
  xlim([threshold_list(1)*0.9 threshold_list(end)*1.1])
xlabel('Threshold (nT/s)')
legend('black-box','grey-box','white-box');legend('boxoff')

  figure(fig4)
  subplot(3,1,station_number)
  plot(threshold_list,hss,'-o','LineWidth',2);  ylabel('Heidke Skill Score'), hold on
 title(char(upper(station)+" station"),'fontsize',12)
  set(gca,'fontsize',14)
  xlim([threshold_list(1)*0.9 threshold_list(end)*1.1])
xlabel('Threshold (nT/s)')
legend('black-box','grey-box','white-box');legend('boxoff')
  
 

  figure(fig5)
  subplot(1,2,1)
  plot(tss_wb,tss_gb,char(markers(station_number)),'markersize',8,'linewidth',2),hold on
  set(gca,'fontsize',16)
  xlabel('Geospace only');ylabel('Geospace + ML')
  title('True Skill Statistic')
  axis([0 0.8 0 0.8]), axis square, grid on

  subplot(1,2,2)
  plot(hss_wb,hss_gb,char(markers(station_number)),'markersize',8,'linewidth',2),hold on
  set(gca,'fontsize',16)
  xlabel('Geospace only');ylabel('Geospace + ML')
  title('Heidke Skill Score')
  axis([0 0.6 0 0.6]), axis square, grid on
  
  
  end 

  figure(fig5)
  subplot(1,2,1)
  plot([0 1],[0 1],'k','linewidth',2)
  legend('FRN','OTT','IQA');legend('boxoff')
  subplot(1,2,2)
  plot([0 1],[0 1],'k','linewidth',2)
  legend('FRN','OTT','IQA');legend('boxoff')

cd(pwd_dir)  

%%%%%%%%%%% COPYRIGHT NOTE %%%%%%%%%%%%%%%%%%

% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

% THIS SOFTWARE IS PROVIDED BY ENRICO CAMPOREALE "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS 
% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE 
% GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 
% STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
