from Imports import *
###
def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin,IndexMin1,IndexMin2,Save):
    
    fig, ax = plt.subplots(3,1,facecolor='w', edgecolor='k')
    fig.set_size_inches([6.69, 8.86], forward=True)
    fig.subplots_adjust(bottom=0.06,top=0.99,left=0.11,right=0.99,hspace=0.3, wspace=0.0)
    ax = ax.ravel()

    gs0 = gridspec.GridSpec(4,1)
    ax[0] = plt.subplot(gs0[2:4,0])
    ax[1] = plt.subplot(gs0[1,0])
    ax[2] = plt.subplot(gs0[0,0])

    starts = [1951.665811,1951.000057,1944.000057]; ends = [1951.748002,1952.000057,1955.000057]
    
    for i in range(3):

        start = np.where(Time[:,4]==starts[i])[0][0]
        end = np.where(Time[:,4]==ends[i])[0][0]

        ax[i].plot(Time[start:end,4],HMC[start:end],color='gray')
        ax[i].plot(Time[start:end,4],Index_thres1[start:end],linestyle='--',color='midnightblue',label=r'$Hl_t\in$ ['+str(np.around(max(Index_thres1[start:end]),2))+' nT ,'+str(np.around(min(Index_thres1[start:end]),2))+' nT]')
        ax[i].scatter(Time[StormIndices,4],HMC[StormIndices],color='midnightblue',s=5,zorder=3)
        
        ax[i].set_xlim([Time[start,4],Time[end,4]])
        ax[i].set_ylabel('HMC [nT]')
        ax[i].legend(ncol=1,frameon=False,fontsize=10,loc=3)    

        if i == 0:
            
            ax[i].plot(Time[start:end,4],HMC[start:end]-100,color='gray')
            ax[i].plot(Time[start:end,4],HMC[start:end]-200,color='gray')
            ax[i].plot(Time[start:end,4],HMC[start:end]-300,color='gray')
            ax[i].scatter(Time[np.asarray(IndexMin),4],HMC[np.asarray(IndexMin)]-100,color='midnightblue',s=10,zorder=4)
            ax[i].scatter(Time[np.asarray(IndexMin1,dtype=int),4],HMC[np.asarray(IndexMin1,dtype=int)]-200,color='midnightblue',s=10,zorder=5)
            ax[i].scatter(Time[np.asarray(IndexMin2,dtype=int),4],HMC[np.asarray(IndexMin2,dtype=int)]-300,color='maroon',s=10,zorder=6)
 
            ax[i].text(0.02,0.94,'Step 1',transform=ax[i].transAxes);ax[i].text(0.02,0.745,'Step 2a',transform=ax[i].transAxes)
            ax[i].text(0.02,0.535,'Step 2b',transform=ax[i].transAxes);ax[i].text(0.02,0.325,'Step 3',transform=ax[i].transAxes)
            ax[i].set_yticks([0,-100,-200,-300])
            
            x_unique, x_indices = np.unique(Time[start:end,2],return_index=True)
            ax[i].set_xticklabels(Time[start:end,2][x_indices][::2])
            ax[i].set_xlabel('Days since '+str(pd.to_datetime(Date[start])))
            ax[i].set_ylim([np.min(HMC[start:end]-300)-20,10])
            ax[i].set_xticks(Time[start:end,4][x_indices],minor=True)
            ax[i].set_xticks(Time[start:end,4][x_indices][::2])

        else:
            if i == 1:
                x_unique, x_indices = np.unique(Time[start:end,1],return_index=True)
                ax[i].set_xticklabels(Time[start:end,1][x_indices])
                ax[i].set_xlabel('Months since '+str(pd.to_datetime(Date[start])))
            else:
                x_unique, x_indices = np.unique(Time[start:end,0],return_index=True)
                ax[i].set_xticklabels(Time[start:end,0][x_indices])
                ax[i].set_xlabel('Years since '+str(pd.to_datetime(Date[start])))

            ax[i].text(0.02,0.3,'Step 1',transform=ax[i].transAxes)
            ax[i].set_yticks([0])
            
            ax[i].set_ylim([np.min(HMC[start:end])-20,10])
            ax[i].axvspan(starts[i-1], ends[i-1], facecolor='silver')
            ax[i].set_xticks(Time[start:end,4][x_indices])

        ax[i].tick_params(axis='x',direction='in')
        ax[i].tick_params(axis='x',which='minor',direction='in')
      
    if Save == True:
        fig.savefig('./Dump/Fig/development/EventSelection.pdf',format='pdf',dpi=200,transparent=True)
        fig.savefig('./Dump/Fig/development/EventSelection.png',format='png',dpi=200,transparent=True)
    plt.show()
###
###
def IndexDist(Time,YearsIndex,Storms,Kp_all,KpHours_all,HMC,Save):
    
    TimeDecYear = Time[:,4]

    # Overlap of storm times with Kp
    Kp_indices = np.where(np.logical_and(KpHours_all >= min(TimeDecYear[YearsIndex]),KpHours_all < max(TimeDecYear[YearsIndex])))[0]
    Kp_red = Kp_all[Kp_indices]; KpHours_red = KpHours_all[Kp_indices]
    Kp_step = 0.3333333
    Kp_range = np.around(np.arange(0,9,Kp_step),3)
    Storms1932 = np.where(TimeDecYear[Storms] >= 1932.000171)[0]
    HMC_step = 7.78
    HMC_range = np.arange(10,480+HMC_step,HMC_step)
    
    Index_range = HMC_range #Kp_range
    step = HMC_step #Kp_step
    
    years = np.arange(1932,2016,1)
    Percentages = np.zeros((len(years),len(Index_range)))
    Hist = np.zeros((len(years),len(Index_range)-1))
    KpValuesAll = []
    
    for l in range(len(years)):
        if l == len(years)-2:
            Kp_indices = np.where(np.logical_and(KpHours_red >= years[l],KpHours_red <= years[l]+1))[0]
            StormsBlock = np.where(np.logical_and(TimeDecYear[Storms] >= years[l],TimeDecYear[Storms] <= years[l]+1))[0]
        else:
            Kp_indices = np.where(np.logical_and(KpHours_red >= years[l],KpHours_red < years[l]+1))[0]
            StormsBlock = np.where(np.logical_and(TimeDecYear[Storms] >= years[l],TimeDecYear[Storms] < years[l]+1))[0]
        #Kp_year = Kp_red[Kp_indices]; KpHours_year = KpHours_red[Kp_indices]
        Kp = Kp_red[Kp_indices]; KpHours = KpHours_red[Kp_indices]

        Kp_hours = np.zeros((len(Kp)*3,2))
        indices = np.where(np.in1d(TimeDecYear,KpHours))[0] 
        j = 0
        for i in range(len(Kp)):
            Kp_hours[j:j+3,0] = TimeDecYear[indices[i]-1:indices[i]+2]
            Kp_hours[j:j+3,1] = Kp[i]
            j+=3

        KpStorms = np.where(np.in1d(np.around(Kp_hours[:,0],6),TimeDecYear[Storms][Storms1932]))[0]
        #KpValues = Kp_hours[KpStorms,1]
        KpValues = np.zeros(len(KpStorms))
        for k in range(len(KpStorms)): 
            if (KpStorms[k]+13 in range(len(Kp_hours)) and KpStorms[k]-13 in range(len(Kp_hours))):
                KpValues[k] = np.nanmax(Kp_hours[KpStorms[k]-13:KpStorms[k]+14,1])
            else:
                KpValues[k] = np.nanmax(Kp_hours[KpStorms[k],1])

        HMCValues = -HMC[Storms[StormsBlock]]  

        #KpValuesAll.extend(KpValues.tolist())
        #Hist[l,:], bin_edges = np.histogram(KpValues, bins = Index_range)

        KpValuesAll.extend(HMCValues.tolist())
        Hist[l,:], bin_edges = np.histogram(HMCValues, bins = Index_range)

        Hist[l,:] *= 100/len(Storms1932)
        
    SumHist = np.sum(Hist,axis=0)
    CumSum = np.cumsum(SumHist)
    CumSum2 = np.cumsum(np.insert(SumHist[0:len(SumHist)-1],0,0))
    bound25 = np.percentile(np.asarray(KpValuesAll).flatten(),25)
    bound75 = np.percentile(np.asarray(KpValuesAll).flatten(),75)
    bound50 = np.percentile(np.asarray(KpValuesAll).flatten(),50)
    #print(bound25,bound50, bound75)

    fig = plt.figure()
    fig.set_size_inches([4, 8], forward=True)
    fig.subplots_adjust(left=0.1,right=0.99)
    gs = gridspec.GridSpec(2,2,height_ratios=[3,1],width_ratios=[20,1],hspace=0.06,wspace=0.04,left=0.2,right=0.85,bottom=0.06,top=0.99)
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[2])
    #ax3 = ax2.twinx()

    cmap = plt.get_cmap('Greys')
    im=ax1.imshow(Hist,cmap=cmap,vmin=np.min(Hist),vmax=np.max(Hist),aspect='auto')
    ax1.set_xticks(np.arange(0,len(Index_range),3)-0.5)
    ax1.set_xticks(np.arange(0,len(Index_range),1)-0.5,minor=True)
    ax1.set_xticklabels([])
    ax1.set_xlim([-0.5,26.5])
    ax1.invert_yaxis()
    ax1.set_ylabel('Years')
    ax1.set_yticks(np.arange(0,len(years)-1,1)-0.5,minor=True)
    ax1.set_yticks(np.arange(0,len(years)-1,4)-0.5)
    ax1.set_yticklabels(years[0:len(years)-1][::4])
    ax1.axvline((bound25-Index_range[0])/step-0.5,color='black',linestyle='--')
    ax1.axvline((bound75-Index_range[0])/step-0.5,color='black',linestyle='--')
    cax=plt.subplot(gs[1])
    bar=plt.colorbar(im,cax=cax,orientation='vertical')
    bar.set_label('Occurrence [%]')

    ax2.bar(np.arange(0,len(Index_range)-1,1),SumHist,width=1.0,align='edge',color='gray',zorder=1)
    ax2.axvline((bound25-Index_range[0])/step,color='black',linestyle='--',label='Q1')
    ax2.axvline((bound75-Index_range[0])/step,color='black',linestyle='--',label='Q3')
    ax2.set_xticks(np.arange(0,28,3))
    ax2.set_xticks(np.arange(0,28,1),minor=True)
    ax2.set_xlim([0,27])
    ax2.legend(loc=4,ncol=1,frameon=False,fontsize=10)
    ax2.set_ylabel('Occurrence [%]')

    #ax2.set_xticklabels(np.arange(0,10,1))
    #ax2.set_xlabel('Kp')
    #ax2.set_ylim([0,10])

    ax2.set_xticklabels(np.around(Index_range[::3]/100,1))
    ax2.set_xlabel('-HMC/100 [nT]')
    ax2.set_ylim([0,20])

    #ax3.fill_between(np.arange(0,len(Kp_range)-1,1)+0.5,CumSum,CumSum2,color='maroon')
    #ax3.set_ylim([0,105])
    #ax3.set_ylabel('Cum. Occurrence [%]')
    #ax3.spines['right'].set_color('maroon')
    #ax3.tick_params(axis='y', colors='maroon')
    #ax3.yaxis.label.set_color('maroon')
    
    if Save == True:
        fig.savefig('./Dump/Fig/development/KpDist_test.pdf',format='pdf',dpi=200,transparent=True)
        fig.savefig('./Dump/Fig/development/KpDist_test.png',format='png',dpi=200,transparent=True)
    plt.show()
###
###
def Diagnostics(n_features,n_classes,NData, Save):

    NCIRs = NData[0]; NCMEs = NData[1]

    fig, ax = plt.subplots(3,4,facecolor='w',edgecolor='k',sharex=True)
    fig.set_size_inches([14,7], forward=True)
    fig.subplots_adjust(bottom=0.08,top=0.95,left=0.05,right=0.98,wspace=0.2,hspace=0.2)
    ax = ax.ravel()
    
    #fNames = [r'$\Delta$HMC', r'dHMC/dt,r', r'$\Delta$dHMC/dt', r'r$\left(\hat{Z}_6,\hat{Z}_{18}\right)$', r'Range r$\left(\hat{Z}_{6\mathrm{ lag}}\mathrm{, }\hat{Z}_{18}\right)$', r'MAD$\left(d\hat{Z}_6/dt\right)$', r'arg max $Z_{\mathrm{MLT}}$', r'max$\left(ASY_{max}\right)$', 'FWHM $ASY_{\mathrm{DD}}$', 'Solar cycle phase', 'Recurrence', 'Ideal (synthetic)']
    fNames = [r'$k$=1',r'$k$=9',r'$k$=11',r'$k$=2',r'$k$=10',r'$k$=3',r'$k$=6',r'$k$=7',r'$k$=8',r'$k$=4',r'$k$=5',r'$k$=12']
    All = sum(len(NData[j]) for j in range(n_classes))

    if n_classes == 2: 
        colors=['midnightblue','maroon']; labels=['C/SIRs','ICMEs']
    else: 
        colors=['midnightblue','maroon','green']; labels=['C/SIRs','CMEs','SIRs']
    
    k = 0
    for i in [0,3,5,9,10,6,7,8,1,4,2,11]: #range(n_features):
 
        for j in range(n_classes):
            hist,bins = np.histogram(NData[j][:,i],np.arange(0,1.01,0.01))
            ax[k].plot(np.arange(0,1.0,0.01)+0.005,hist*100/All,color=colors[j],linewidth=0.75,label=labels[j])

            xmin = np.percentile(NData[j][:,i],25); xmax = np.percentile(NData[j][:,i],75)
            ax[k].axvspan(xmin,xmax,0,1,color=colors[j],alpha=0.2,zorder=0)
            ax[k].axvline(np.nanmedian(NData[j][:,i]),0,1,color=colors[j],linestyle='--')
        
        if k == 0:
            # dummy
            ax[k].axvspan(2,3,0,1,color='black',alpha=0.1,zorder=0,label='IQR')
            ax[k].axvline(2,0,1,color='black',linestyle='--',label='Q2')
            ax[k].legend(ncol=1,loc=1,fontsize=18,frameon=False)

        if k == 8:
            ax[k].set_ylabel('Occurrence [%]', fontsize=18)
            ax[k].set_xlabel('Standardized feature value', fontsize=18)
        #if i == n_features-1: ax[i].text(0.26,1.02,'Synthetic',fontsize=8,transform=ax[i].transAxes)
        #else: 
        #ax[k].text(0,1.02,str(k+1)+': '+fNames[i],fontsize=18,transform=ax[k].transAxes)
        ax[k].text(0,1.02,fNames[i],fontsize=18,transform=ax[k].transAxes)
        
        ax[k].set_xlim([0,1]); ax[k].set_xticks([0,0.25,0.5,0.75,1.0])
        ax[k].tick_params(axis = 'y', which='both',direction = 'in',labelsize=16)
        ax[k].tick_params(axis = 'x', which='both',direction = 'inout',labelsize=16)
        
        k += 1
    
    if Save == True:
        fig.savefig('./Dump/Fig/development/Features.pdf',format='pdf',dpi=200,transparent=True)
        #fig.savefig('./Dump/Fig/Features.png',format='png',dpi=200,transparent=True)
    plt.show()
###
###
def Ranking(Scores, Std, Data, Ranking, nfeaturesFinal, estimators, scoring_name, est, ax, Save, SaveName):
    
    nEstimators = len(estimators)

    if est == 0:
        fig, ax = plt.subplots(1,nEstimators,facecolor='w',edgecolor='k',sharey=True)
        fig.set_size_inches([8,4],forward=True)
        fig.subplots_adjust(bottom=0.14,top=0.85,left=0.11,right=0.99,wspace=0.05,hspace=0.0)
        if nEstimators > 1:
            ax = ax.ravel()
        else:
            ax = [ax]

    ax[est].errorbar(np.arange(1,Data.shape[1]+1,1),Scores[:,0],color='midnightblue',yerr=2*Std[:,0],linewidth=1.5,elinewidth=1.5,label='Validation set')
    ax[est].errorbar(np.arange(1,Data.shape[1]+1,1),Scores[:,1],color='gray',yerr=2*Std[:,1],linewidth=1.5,elinewidth=1.5,label='Training set')
    ax[est].set_xticks(np.arange(1,Data.shape[1]+1,1))
    ax[est].set_xlim([0.5,11.5])
    ax[est].scatter(np.argmax(Scores[:,0])+1,np.max(Scores[:,0]),s=50,color='maroon',zorder=3,label='Best')
    #ax[est].scatter(np.argmax(Scores[:,1])+1,np.max(Scores[:,1]),s=30,color='gray',zorder=3,label='Best (Training)')
    #ax[est].scatter(nfeaturesFinal,Scores[nfeaturesFinal-1,0],s=20,color='midnightblue',zorder=2,label='Best (Validation)')
    #ax[est].text(0.0,1.15,estimators[est],transform=ax[est].transAxes,fontsize=12,fontweight='bold')
    ax[est].tick_params(axis ='x',which='both',direction='inout',labelsize=16)
    ax[est].tick_params(axis ='y',which='both',direction='inout',labelsize=16)

    axu=ax[est].twiny()
    axu.set_xticks(np.arange(1,Data.shape[1]+1,1))
    axu.set_xticklabels(Ranking)
    axu.set_xlim([0.5,11.5])
    axu.tick_params(axis ='x',which='both',direction='inout',labelsize=16)

    if est == 0:
        ax[est].set_xlabel('Number of features',fontsize=18)
        ax[est].legend(loc=4,ncol=1,frameon=False,fontsize=18,scatterpoints=1)
        ax[est].set_ylabel(scoring_name,fontsize=18)
        axu.set_xlabel(r'Feature IDs ( $\rightarrow$ cumulative)',fontsize=18)
    
    if est == nEstimators-1:
      
        if Save == True:
            fig.savefig('./Dump/Fig/development/Ranking_'+SaveName+'.pdf',format='pdf',dpi=300,transparent=True)
            #fig.savefig('./Dump/Fig/Ranking.png',format='png',dpi=300,transparent=True)
        plt.show()
        
    return ax
###
###
def ModelCounts(nModels,SumCounts,CV2_means,N1,K1,N2,K2,Param_Comb,BestID,scoring_name,Save,SaveName):

    fig, ax = plt.subplots(1,1)
    fig.set_size_inches([7,5],forward=True)
    fig.subplots_adjust(bottom=0.15,top=0.98,left=0.1,right=0.85)
    ax.bar(np.arange(0,nModels,1)-0.25,SumCounts[:,0]*100/(N1*K1*N2*K2),width=0.5,color='darkgray')
    ax.set_ylim([0,max(SumCounts[:,0])*100/(N1*K1*N2*K2)])
    ax.yaxis.grid(False)
    ax.set_ylabel('Count out of '+str(N1*K1*N2*K2)+' runs [%]', fontsize=18)
    #ax.set_title('Model selection statistics')
    ax.spines['right'].set_color('darkgray')
    ax.tick_params(axis='y', colors='darkgray')
    ax.yaxis.label.set_color('darkgray')
    ax.set_xlabel('Model ID', fontsize=18)
    ax.set_xticks(np.arange(0,nModels,1))
    ax.set_xlim(-0.5,nModels-1+0.5)
    ax.tick_params(axis ='x',which='both',direction='inout',labelsize=16)
    ax.tick_params(axis ='y',which='both',direction='inout',labelsize=16)

    ax1 = ax.twinx()
    ax1.set_ylim([min(CV2_means[np.nonzero(CV2_means[:,0]),0].squeeze())-0.01,max(CV2_means[:,0])])
    ax1.bar(np.arange(0,nModels,1)+0.25,CV2_means[:,0],width=0.5,color='black')
    ax1.bar(BestID+0.25,CV2_means[BestID,0],width=0.5,color='maroon',edgecolor='black',linewidth=2)
    ax1.yaxis.grid(False)
    ax1.set_ylabel(scoring_name, fontsize=18)
    ax1.tick_params(axis ='y',which='both',direction='inout',labelsize=16)

    k = 0
    ParamKeys = list(Param_Comb[0].keys())
    nKeys = len(ParamKeys)
    for i in range(nModels):
        if i == BestID:
            for j in range(nKeys):
                if ParamKeys[j] != 'select__n_features_to_select':
                    ax1.text(i+0.25,0.52+k*0.012,r'O'+str(k+1)+' ='+str(Param_Comb[i][ParamKeys[j]]),rotation=90,fontsize=16, color='white', horizontalalignment = 'center')
                    #print(Param_Comb[i][ParamKeys[j]])
                    k += 1

    if Save == True:
        fig.savefig('./Dump/Fig/development/ModelStats_'+SaveName+'.pdf',format='pdf',dpi=300,transparent=True)
        #fig.savefig('./Dump/Fig/ModelStats.png',format='png',dpi=300,transparent=True)
    plt.show()
###
###
def CM(CM,Condition,n_classes,Save,SaveName):

    CM = np.reshape(CM,(n_classes,n_classes))
    #CM_names = ['TNR','FPR','FNR','TPR']
    #CM_names = ['NPV','FDR','FOR','PPV']
    CM_names = ['TN','FP,','FN','TP']
    CM_names = np.asarray(CM_names).reshape((2,2))

    classes = range(n_classes)
    class_names = ['0: C/SIR', '1: ICME']

    #print(Condition)
    print(CM)

    ##Normalization
    #for i in range(n_classes):
        ##CM[i,:] /= Condition[i,0]
        #CM[:,i] /= Condition[i,1]
        #class_names.append(np.str(i))

    fig, ax = plt.subplots(1,1,facecolor='w',edgecolor='k',sharey=True)
    fig.set_size_inches([4,4], forward=True)
    fig.subplots_adjust(bottom=0.01,top=0.99,left=0.2,right=0.9,wspace=0.0,hspace=0.0)

    cmap=plt.cm.Blues
    Im = ax.imshow(CM, interpolation='nearest', cmap=cmap)
    CBar = plt.colorbar(Im,shrink=0.66)

    thresh = CM.max()/2.
    for i, j in itertools.product(range(CM.shape[0]), range(CM.shape[1])):
        ax.text(j, i, CM_names[i,j], horizontalalignment='center', verticalalignment = 'bottom', color='white' if CM[i, j] > thresh else "black", fontsize=16)
        ax.text(j, i, format(CM[i, j], '.3f'), horizontalalignment='center', verticalalignment = 'top', color='white' if CM[i, j] > thresh else "black", fontsize=16)

    #ax.set_title('Confusion matrix')
    ax.set_xticks(classes); ax.set_yticks(classes)
    ax.set_xticklabels(class_names); ax.set_yticklabels(class_names,rotation=90)
    ax.set_xlabel('Predicted class',fontsize=18); ax.set_ylabel('True class',fontsize=18)
    ax.tick_params(axis ='x',which='both',direction='inout',labelsize=16)
    ax.tick_params(axis ='y',which='both',direction='inout',labelsize=16)
    CBar.set_label('Normalized count',fontsize=18)
    CBar.ax.tick_params(axis ='y',which='both',direction='inout',labelsize=16)
    ax.grid(False)

    if Save == True:
        fig.savefig('./Dump/Fig/development/Confusion_'+SaveName+'.pdf',format='pdf',dpi=300,transparent=True)
        #fig.savefig('./Dump/Fig/ConfusionCol.png',format='png',dpi=300,transparent=True)
    
    plt.show()
###
###
def Curves(N2,K2,Curves,curve_i,Model_Mean,Model_Std,C,Save,SaveName):

    fig, axs = plt.subplots(1,2,facecolor='w',edgecolor='k',sharey=True)
    fig.set_size_inches([14,6], forward=True)
    fig.subplots_adjust(bottom=0.1,top=0.95,left=0.1,right=0.9,wspace=0.08,hspace=0.0)
    axs = axs.ravel()
    axs[0].axis('equal')
    axs[1].axis('equal')

    precision = Curves[:,:,0,:]
    tpr = Curves[:,:,1,:]

    for i in range(N2):
        for ii in range(K2):

            #tpr[i,ii,0] = 0.0; tpr[i,ii,-1] = 1.0
            #precision[i,ii,0] = 1.0; precision[i,ii,-1] = 0.0
            if i== 0 and ii == 0:
                axs[0].plot(curve_i,tpr[i,ii,:],color='gray',label='1 fold',linewidth=0.5,zorder=0)
            else:
                axs[0].plot(curve_i,tpr[i,ii,:],color='gray',linewidth=0.5,zorder=0)
            axs[1].plot(curve_i,precision[i,ii,:],color='gray',linewidth=0.5,zorder=0)

        ROC_inner = np.mean(tpr[i,:,:],axis=0)
        PR_inner = np.mean(precision[i,:,:],axis=0)
        if i == 0:
            axs[0].plot(curve_i,ROC_inner,color='black',label='Split mean (4 folds)')
        else:
            axs[0].plot(curve_i,ROC_inner,color='black',zorder=1)
        axs[1].plot(curve_i,PR_inner,color='black',zorder=1)
    
    ROC_inner = np.nanmean(tpr,axis=1); ROC_outer = np.nanmean(ROC_inner,axis=0)
    PR_inner = np.nanmean(precision,axis=1); PR_outer = np.nanmean(PR_inner,axis=0)

    axs[0].plot(curve_i,ROC_outer,color='maroon',label=r'AUC = '+str(np.around(Model_Mean[1],3))+'$\pm$'+str(np.around(Model_Std[1],5)),zorder=2)
    axs[1].plot(curve_i,PR_outer,color='maroon',label=r'AUC = '+str(np.around(Model_Mean[0],3))+'$\pm$'+str(np.around(Model_Std[0],5)),zorder=2)

    P = C[1,0]; PP = C[1,1]; N = C[0,0]; PN = C[0,1]; POP = sum(C[:,0])
    #P,PP,N,PN = C
    axs[0].plot([0,1],[0,1],linestyle='--',color='midnightblue',label='Chance')
    #axs[0].axvspan(0,1,0,1,color='silver',alpha=0.5,zorder=0,lw=0)
    axs[1].plot([0,1],[P/(P+N),P/(P+N)],linestyle='--',color='midnightblue',label='Prevalence = P/POP')
    #axs[1].axvspan(0,1,P/(P+N),1,color='silver',alpha=0.5,zorder=0,lw=0)
    axs[0].set_xlabel(r'FPR = FP/N', fontsize=18); axs[1].set_xlabel(r'TPR = TP/P', fontsize=18)
    axs[0].set_ylabel(r'TPR = TP/P', fontsize=18); axs[1].set_ylabel(r'PPV = TP/PP', fontsize=18)
    axs[0].legend(loc=0, frameon=False, fontsize=16); axs[1].legend(loc=0, frameon=False, fontsize=16)
    axs[0].set_title('ROC curve',fontsize=18)
    axs[0].tick_params(axis ='x',which='both',direction='inout',labelsize=16)
    axs[0].tick_params(axis ='y',which='both',direction='inout',labelsize=16)
    
    axs[1].set_title('Precision-Recall curve',fontsize=18)
    axs[1].yaxis.set_label_position('right'); axs[1].yaxis.set_ticks_position('right')
    axs[1].tick_params(axis ='x',which='both',direction='inout',labelsize=16)
    axs[1].tick_params(axis ='y',which='both',direction='inout',labelsize=16)
    axs[0].set_xlim([0,1]);axs[0].set_ylim([0,1.03]);axs[1].set_xlim([0,1]);axs[1].set_ylim([0,1.03])

    if Save == True:
        fig.savefig('./Dump/Fig/development/Curves_'+SaveName+'.pdf',format='pdf',dpi=300,transparent=True)
        #fig.savefig('./Dump/Fig/Curves.png',format='png',dpi=300,transparent=True)

    plt.show()
###
###
def Musical(Time,HTime,Storms,Kp_all,KpHours_all,SN,SNYears,HMC,HMC_filt,Training,Save):

    fig, axs = plt.subplots(1,4,facecolor='w', edgecolor='k',sharey=True)
    fig.set_size_inches([6.69, 8.86], forward=True)
    fig.subplots_adjust(bottom=0.06,top=0.94,left=0.115,right=0.99,hspace=0.00,wspace=0.03)

    gs0 = gridspec.GridSpec(1,5)
    axs1 = plt.subplot(gs0[0,0:4])
    axs2 = plt.subplot(gs0[0,4])#,sharey=axs1)
    axs3 = axs2.twiny()

    axs1.set_facecolor('white')
    axs2.set_facecolor('white')

    axs1.grid(False)

    #############################################

    SolarRot = 27.27
    Years = np.unique(HTime[:,0])
    Days = len(Years)*365+28
    Lines = int(np.ceil(Days/SolarRot))

    #############################################

    axs1.set_xlim([0,SolarRot])
    axs1.set_xticks(np.arange(0,28,1),minor=True)
    axs1.set_xticks(np.arange(0,28,5))
    axs1.set_xlabel('Day of Solar Rotation')

    # Plot y-ticks (years)
    y_y = np.zeros(len(Years))
    for i in range(len(Years)):
        y_y_diff = dt.datetime(year=int(Years[i]),month=1,day=1,hour=0,minute=0,second=0) - dt.datetime(1900,1,1,0,0,0)
        y_y[i] = np.divide(y_y_diff.total_seconds()/(60*60*24),SolarRot)

    axs1.set_yticks(y_y,minor=True) # every year
    axs1.set_yticks(y_y[::5])       # every 5ths year
    axs1.set_yticklabels(Years[::5])
    axs1.set_ylabel('Year')
    axs1.tick_params(axis='x',which='both',bottom=False,top=True,labelbottom=False,labeltop=True)
    axs1.xaxis.set_label_position('top') 

    plt.setp(axs2.get_yticklabels(), visible=False)
    axs2.set_yticks(Years,minor=True)
    axs2.set_yticks(Years[::5])
    axs2.tick_params(axis='y',direction='in')
    axs2.tick_params(axis='y',which='minor',direction='in')
    axs2.set_yticklabels([])
    axs2.set_xlim([0,125])
    axs2.set_xticks([0,25,50,75,100,125],minor=True)
    axs2.set_xticks([0,50,100])
    axs2.set_ylim([HTime[Storms[0],4],HTime[Storms[-1],4]])
    axs2.set_xlabel('# Storms')
    axs2.spines['bottom'].set_color('black')
    #axs2.spines['top'].set_color('midnightblue')
    
    axs3.set_ylim([HTime[Storms[0],4],HTime[Storms[-1],4]])
    axs3.set_xlim([0,1])
    axs3.set_xticks([0,0.25,0.5,0.75,1],minor=True)
    axs3.set_xticks([0,0.5,1])
    axs3.set_xticklabels(['0','0.5','1'])
    axs3.set_xticks([0,0.25,0.5,0.75,1],minor=True)
    axs3.set_xlabel('Normalized SC')
    
    #axs3.spines['bottom'].set_color('black')
    axs3.spines['top'].set_color('dimgray')

    ####################################################################

    # Plot Kp events
    Kp_indices = np.where(np.logical_and(Kp_all <= 4.666, Kp_all >= 3.0))[0]
    Kp = Kp_all[Kp_indices]; KpHours = KpHours_all[Kp_indices]
    y_Kp = np.zeros(len(Kp)); x_Kp = np.zeros(len(Kp))
    KpTime = np.where(np.in1d(HTime[:,4],KpHours))[0]
    for i in range(len(KpTime)):
        Kp_diff = dt.datetime(year=int(HTime[KpTime[i],0]),month=int(HTime[KpTime[i],1]),day=int(HTime[KpTime[i],2]),hour=int(HTime[KpTime[i],3]),minute=30,second=30)-dt.datetime(1900,1,1,0,0,0)
        y_Kp[i] = np.divide(Kp_diff.total_seconds()/(60*60*24),SolarRot)
        x_Kp[i] = np.mod(Kp_diff.total_seconds()/(60*60*24),SolarRot)
    #axs1.barh(y_Kp,3/24,left=x_Kp,color='black',edgecolor='black',align='edge',height=1/Lines)
    #axs1.scatter(x_Kp,y_Kp,s=5,marker='s',color='dimgray',alpha=0.7,zorder=1)

    ## Dummy for legend
    #axs1.scatter(0,-50,s=20,marker='s',color='black',alpha=0.7,label=r'3 $\leq$ Kp $\leq$ 4+$^+$')
    axs1.plot([0,1],[-50,-25],color='dimgray',linewidth=1.5,label=r'SN')
    axs1.plot([0,1],[-50,-25],color='dimgray',linewidth=1.5,linestyle='--',label=r'HMC$_{11y}$')
    axs1.scatter(0,-50,s=20,marker='s',color='dimgray',label=r'Event')

    # Plot storm events from catalogs (papers)
    cm = plt.cm.get_cmap('seismic')

    KnownE = np.where(Training[:,2]==1)[0]
    diff_t = (pd.to_datetime(Time[Training[KnownE,1]])-dt.datetime(1900,1,1,0,0,0))
    y_t = np.zeros(len(Training[KnownE,1])); x_t = np.zeros(len(Training[KnownE,1])) 
    for i in range(len(Training[KnownE,1])):
        y_t[i] = np.divide(diff_t[i].total_seconds()/(60*60*24),SolarRot)
        x_t[i] = np.mod(diff_t[i].total_seconds()/(60*60*24),SolarRot)

    z = Training[KnownE,0]
    sc = axs1.scatter(x_t, y_t, c=z, vmin=0.0, vmax=1.0, s=5, cmap=cm, zorder = 3)
    cbar_axis = fig.add_axes([0.45, 0.032, 0.3, 0.01]) 
    cbar = plt.colorbar(sc,cax=cbar_axis,use_gridspec=False, orientation ='horizontal',ticks=[0,0.5,1])
    #cbar.set_label('CME probability',fontsize=8)
    cbar.ax.xaxis.set_ticks(np.arange(0,1.1,0.1), minor=True)
    cbar.ax.set_xticklabels(['0:C/SIR','ICME probability','1:ICME'],fontsize=11)
    cbar.ax.tick_params(axis='x',which='both',direction='in')

    #CIRs = np.where(Training[KnownE,0] == 0)[0]
    #CMEs = np.where(Training[KnownE,1] == 1)[0]

    # Plots class occurrences
    nip = 0
    for i in [1,0]:
        zi = np.where(z==i)[0]
        ni, bins = np.histogram(HTime[Training[KnownE[zi],1],4],np.arange(1900,2017,1))
        axs2.barh(np.arange(1900,2016,1), ni, left = nip, height=0.8,color=cm(float(i)),zorder=5)
        nip += ni
            
    #Turner_CIR = np.where(np.logical_and(Training[:,2] == b'Turner',Target['type'] ==0))[0]
    #Jian_CIR =  np.where(np.logical_and(Target['source'] == b'Jian',Target['type'] ==0))[0]
    #Shen_CIR = np.where(np.logical_and(Target['source'] == b'Shen',Target['type'] ==0))[0]
    #Turner_CME = np.where(np.logical_and(Target['source'] == b'Turner',Target['type'] ==1))[0]
    #Jian_CME =  np.where(np.logical_and(Target['source'] == b'Jian',Target['type'] ==1))[0]
    #Shen_CME = np.where(np.logical_and(Target['source'] == b'Shen',Target['type'] ==1))[0]

    #axs1.scatter(x_t[Turner_CIR],y_t[Turner_CIR],s=20,marker='x',color='dodgerblue')
    #axs1.scatter(x_t[Jian_CIR],y_t[Jian_CIR],s=20,marker='x',color='orange')
    #axs1.scatter(x_t[Shen_CIR],y_t[Shen_CIR],s=20,marker='x',color='limegreen')
    #axs1.scatter(x_t[Turner_CME],y_t[Turner_CME],s=20,marker='d',color='dodgerblue')
    #axs1.scatter(x_t[Jian_CME],y_t[Jian_CME],s=20,marker='d',color='orange')
    #axs1.scatter(x_t[Shen_CME],y_t[Shen_CME],s=20,marker='d',color='limegreen')
   
    ## Dummy for legend
    #axs1.scatter(0,-50,s=20,color='dodgerblue',label='Turner') 
    #axs1.scatter(0,-50,s=20,color='orange',label='Jian')       
    #axs1.scatter(0,-50,s=20,color='limegreen',label='Shen')    
    #axs1.scatter(0,-50,s=20,marker='d',color='black',label='CME') 
    #axs1.scatter(0,-50,s=20,marker='x',color='black',label='CIR') 
    #axs1.scatter(0,-50,s=20,color='black',label='Targets') 
    #axs1.scatter(0,-50,s=20,marker='s',color='gray',label=r'Kp $\geq$ 4$^+$')
    #axs1.plot([0,1],[-50,-25],color='blue',linewidth=1.5,label=r'HMC$_{11}$,p=1')
    #axs1.plot([0,1],[-50,-25],color='black',linewidth=1.5,label=r'SN')

    # Plot storm events (my picks)
    #z = np.zeros(len(Storms))+0.5
    diff = (pd.to_datetime(Time[Storms])-dt.datetime(1900,1,1,0,0,0))
    y = np.zeros(len(Storms)); x = np.zeros(len(Storms))
    for i in range(len(Storms)):
        y[i] = np.divide(diff[i].total_seconds()/(60*60*24),SolarRot)
        x[i] = np.mod(diff[i].total_seconds()/(60*60*24),SolarRot)

    #sc = axs1.scatter(x, y, c=z, vmin=0.0, vmax=1.0, s=3, cmap=cm, zorder=2)
    sc = axs1.scatter(x, y, s=3, color='gray', zorder=2)

    handles, labels = axs1.get_legend_handles_labels()
    order = [0,1]
    axs1.legend([handles[idx] for idx in order],[labels[idx] for idx in order], ncol=5,frameon=True,fontsize=11,loc='center left',bbox_to_anchor=[-0.1,-0.037])
    axs1.set_ylim([np.nanmin(y),np.nanmax(y)])
    
    #handles, labels = axs1.get_legend_handles_labels()
    #order = [4,3,2,7,6,5,8,0,1]
    #axs1.legend([handles[idx] for idx in order],[labels[idx] for idx in order], ncol=5,frameon=True,fontsize=9,loc='center left',bbox_to_anchor=[-0.09,-0.037])
    #axs1.set_ylim([np.nanmin(y),np.nanmax(y)])

    ################################################
    y_0 = np.arange(0,Lines)*SolarRot
    for i in range(len(y_0)):
        if i%13 == 0:
            axs1.plot(np.linspace(0,SolarRot,100),np.linspace(y_0[i],y_0[i+1],100)/SolarRot,linewidth=0.3,color='gray',zorder=0)
    ################################################

    # Plot yearly sunspot number
    SN_start = np.where(SNYears == 1930.5)[0][0]
    SN_end = np.where(SNYears == 2016.5)[0][0]
    #axs3.plot((SN[SN_start:SN_end,1]-np.min(SN[SN_start:SN_end,1]))/np.ptp(SN[SN_start:SN_end,1]),SN[SN_start:SN_end,0],color='black',linewidth=1.5, label='SN')
    axs3.plot((SN[SN_start:SN_end]/max(SN[SN_start:SN_end]))**1,SNYears[SN_start:SN_end],color='dimgray',linewidth=1.5, label='SN')

    # Plot HMC threshold scaling
    axs3.plot((HMC_filt/min(HMC_filt))**1,HTime[:,4],color='dimgray',linestyle='--',linewidth=1.5, label='HMC')
    
    # Plot histogram
    n, bins, patches = axs2.hist(x=HTime[Storms,4],bins=np.arange(1900,2017,1),orientation='horizontal',histtype='stepfilled',color='silver',alpha=1)
    #print(n)
   
    if Save == True:
        fig.savefig('./Dump/Fig/development/StormsSetup.pdf',format='pdf',dpi=200)#,transparent=True)
        #fig.savefig('./Dump/Fig/development/StormsSetup.png',format='png',dpi=200)#,transparent=True)
    plt.show()
###
###
def MusicalClassified(Time,HTime,Storms,Class_Proba,Kp_all,KpHours_all,SN,SNYears,HMC,HMC_filt,Save,SaveName):

    fig, axs = plt.subplots(1,4,facecolor='w', edgecolor='k',sharey=True)
    fig.set_size_inches([6.69, 8.86], forward=True)
    fig.subplots_adjust(bottom=0.06,top=0.94,left=0.115,right=0.99,hspace=0.00,wspace=0.03)

    gs0 = gridspec.GridSpec(1,5)
    axs1 = plt.subplot(gs0[0,0:4])
    axs2 = plt.subplot(gs0[0,4])#,sharey=axs1)
    axs3 = axs2.twiny()

    axs1.set_facecolor('silver')
    axs2.set_facecolor('white')

    #############################################

    SolarRot = 27.27
    Years = np.unique(HTime[:,0])
    Days = len(Years)*365+28
    Lines = int(np.ceil(Days/SolarRot))

    #############################################

    axs1.set_xlim([0,SolarRot])
    axs1.set_xticks(np.arange(0,28,1),minor=True)
    axs1.set_xticks(np.arange(0,28,5))
    axs1.set_xlabel('Day of Solar Rotation')

    # Plot y-ticks (years)
    y_y = np.zeros(len(Years))
    for i in range(len(Years)):
        y_y_diff = dt.datetime(year=int(Years[i]),month=1,day=1,hour=0,minute=0,second=0) - dt.datetime(1900,1,1,0,0,0)
        y_y[i] = np.divide(y_y_diff.total_seconds()/(60*60*24),SolarRot)

    axs1.set_yticks(y_y,minor=True) # every year
    axs1.set_yticks(y_y[::5])       # every 5ths year
    axs1.set_yticklabels(Years[::5])
    axs1.set_ylabel('Year')
    axs1.tick_params(axis='x',which='both',bottom=False,top=True,labelbottom=False,labeltop=True)
    axs1.xaxis.set_label_position('top') 

    plt.setp(axs2.get_yticklabels(), visible=False)
    axs2.set_yticks(Years,minor=True)
    axs2.set_yticks(Years[::5])
    axs2.tick_params(axis='y',direction='in')
    axs2.tick_params(axis='y',which='minor',direction='in')
    axs2.set_yticklabels([])
    axs2.set_xlim([0,1])
    axs2.set_xticks([0,0.25,0.5,0.75,1],minor=True)
    axs2.set_xticks([0,0.5,1])
    axs2.set_xticklabels(['0','0.5','1'])
    axs2.set_ylim([HTime[Storms[0],4],HTime[Storms[-1],4]])
    axs2.set_xlabel('Occurrence [%]')
    #axs2.spines['bottom'].set_color('khaki')
    ##axs2.spines['top'].set_color('midnightblue')
    
    axs3.set_ylim([HTime[Storms[0],4],HTime[Storms[-1],4]])
    axs3.set_xlim([0,1])
    axs3.set_xticks([0,0.25,0.5,0.75,1],minor=True)
    axs3.set_xticks([0,0.5,1])
    axs3.set_xticklabels(['0','0.5','1'])
    axs3.set_xticks([0,0.25,0.5,0.75,1],minor=True)
    axs3.set_xlabel('Normalized SC')
    #axs3.spines['bottom'].set_color('khaki')
    #axs3.spines['top'].set_color('midnightblue')

    ####################################################################

    # Plot Kp events
    Kp_indices = np.where(np.logical_and(Kp_all <= 4.666, Kp_all >= 3.0))[0]
    Kp = Kp_all[Kp_indices]; KpHours = KpHours_all[Kp_indices]
    y_Kp = np.zeros(len(Kp)); x_Kp = np.zeros(len(Kp))
    KpTime = np.where(np.in1d(HTime[:,4],KpHours))[0]
    for i in range(len(KpTime)):
        Kp_diff = dt.datetime(year=int(HTime[KpTime[i],0]),month=int(HTime[KpTime[i],1]),day=int(HTime[KpTime[i],2]),hour=int(HTime[KpTime[i],3]),minute=30,second=30)-dt.datetime(1900,1,1,0,0,0)
        y_Kp[i] = np.divide(Kp_diff.total_seconds()/(60*60*24),SolarRot)
        x_Kp[i] = np.mod(Kp_diff.total_seconds()/(60*60*24),SolarRot)
    #axs1.barh(y_Kp,3/24,left=x_Kp,color='black',edgecolor='black',align='edge',height=1/Lines)
    #axs1.scatter(x_Kp,y_Kp,s=5,marker='s',color='black',alpha=0.7)
   
    ## Dummy for legend
    #axs1.scatter(0,-50,s=20,marker='s',color='black',alpha=0.7,label=r'Kp $\geq$ 4$^+$')
    axs1.plot([0,1],[-50,-25],color='dimgray',linewidth=1.5,label=r'SN', zorder=5)
    axs1.plot([0,1],[-50,-25],color='dimgray',linewidth=1.5,linestyle='--',label=r'HMC$_{11y}$', zorder=5)
    axs1.scatter(0,-50,s=20,marker='s',color='black',label=r'Event')

    # Plot storm events (my picks)
    CIRs_pred = np.where(Class_Proba[:,0] > 0.5)[0]
    CMEs_pred = np.where(Class_Proba[:,1] > 0.5)[0]
    CIRs_med = np.median(Class_Proba[CIRs_pred][:,1])
    CMEs_med = np.median(Class_Proba[CMEs_pred][:,1])

    z = Class_Proba[:,1]
    #print(pd.to_datetime(Time[Storms]))
    diff = (pd.to_datetime(Time[Storms])-dt.datetime(1900,1,1,0,0,0))
    y = np.zeros(len(Storms)); x = np.zeros(len(Storms))
    for i in range(len(Storms)):
        y[i] = np.divide(diff[i].total_seconds()/(60*60*24),SolarRot)
        x[i] = np.mod(diff[i].total_seconds()/(60*60*24),SolarRot)

    cm = plt.cm.get_cmap('seismic')
    sc = axs1.scatter(x,y,c=z, vmin=0.0, vmax=1.0, s=3, cmap=cm)
    cbar_axis = fig.add_axes([0.45, 0.032, 0.3, 0.01]) 
    cbar = plt.colorbar(sc,cax=cbar_axis,use_gridspec=False, orientation ='horizontal',ticks=[0,0.5,1])
    #cbar.set_label('CME probability',fontsize=8)
    cbar.ax.xaxis.set_ticks(np.arange(0,1.1,0.1), minor=True)
    cbar.ax.set_xticklabels(['0:C/SIR','ICME probability','1:ICME'],fontsize=11)
    cbar.ax.tick_params(axis='x',which='both',direction='in')
    cbar_axis.annotate('',(CIRs_med,0.5),(CIRs_med,2.5), arrowprops=dict(arrowstyle='simple',color=cm(CIRs_med),lw=0.5))
    cbar_axis.annotate('',(CMEs_med,0.5),(CMEs_med,2.5), arrowprops=dict(arrowstyle='simple',color=cm(CMEs_med),lw=0.5))
    
    handles, labels = axs1.get_legend_handles_labels()
    order = [0,1]
    axs1.legend([handles[idx] for idx in order],[labels[idx] for idx in order], ncol=5,frameon=True,fontsize=11,loc='center left',bbox_to_anchor=[-0.1,-0.037])
    axs1.set_ylim([np.nanmin(y),np.nanmax(y)])

    ################################################
    y_0 = np.arange(0,Lines)*SolarRot
    for i in range(len(y_0)):
        if i%13 == 0:
            axs1.plot(np.linspace(0,SolarRot,100),np.linspace(y_0[i],y_0[i+1],100)/SolarRot,linewidth=0.3,color='gray')
    ################################################

    # Plots class occurrences
    n, bins = np.histogram(HTime[Storms,4],np.arange(1900.5,2017.5,1))
    start = 1; step=0.1; loops = 1./step; nip = 0
    for i in range(int(loops)):
        zi = np.where(np.logical_and(z>start-step,z<=start))[0]
        ni, bins = np.histogram(HTime[Storms[zi],4],np.arange(1900.5,2017.5,1))
        axs2.barh(np.arange(1900.5,2016.5,1), (ni*100/n)/100, left = nip, height=0.8,color=cm(1-(i*step)-step/2)) #left = (nip*100/n)/100
        start -= step
        nip += (ni*100/n)/100
    #nCIR, bins = np.histogram(HTime['dyear'][Storms[CIRs_pred,1]],np.arange(1900,2017,1))
    #nCME, bins = np.histogram(HTime['dyear'][Storms[CMEs_pred,1]],np.arange(1900,2017,1))
    #axs2.barh(np.arange(1900,2016,1), (nCIR*100/n)/100, left = (nCME*100/n)/100, height=0.8,color='blue',alpha=0.5)
    #axs2.barh(np.arange(1900,2016,1), (nCME*100/n)/100, height=0.8,color='red',alpha=0.5)

    boundaries = [1936.0,1942.0,1946.0,1952.0,1956.0,1962.0,1967.0,1973.0,1978.0,1984.0,1988.0,1994.0,1998.0,2004.0,2011.0,2015.0]
    for i in range(len(boundaries)):
        axs3.axhline(xmin=0,xmax=1, y=boundaries[i], linewidth=1.5, color='black',zorder=5)

    # Plot yearly sunspot number
    SN_start = np.where(SNYears == 1930.5)[0][0]
    SN_end = np.where(SNYears == 2016.5)[0][0]
    #axs3.plot((SN[SN_start:SN_end,1]-np.min(SN[SN_start:SN_end,1]))/np.ptp(SN[SN_start:SN_end,1]),SN[SN_start:SN_end,0],color='black',linewidth=1.5, label='SN')
    axs3.plot((SN[SN_start:SN_end]/max(SN[SN_start:SN_end]))**1,SNYears[SN_start:SN_end],color='white',linewidth=1.5, label='SN')
    axs3.plot((HMC_filt/min(HMC_filt))**1,HTime[:,4],color='white',linestyle='--',linewidth=1.5, label='HMC')
    
    if Save == True:
        fig.savefig('./Dump/Fig/development/StormsClassified_'+SaveName+'.pdf',format='pdf',dpi=200)#,transparent=True)
        #fig.savefig('./Dump/Fig/StormsClassified.png',format='png',dpi=200)#,transparent=True)
    
    plt.show()
###
###
def Stats(HTime,Storms,SN,SNYears,HMC,Proba,Class,Save,SaveName):

    SNYearsRed = np.where(np.logical_and(SNYears-0.5>=1930,SNYears-0.5<=2015))[0]
    SN = SN[SNYearsRed]
    SNYears = SNYears[SNYearsRed]-0.5
    
    #HMCB = [-40,-60,-80,-110,-160,-478]
    HMCB = [-10,-40,-50,-65,-90,-500]
    Storms1932 = np.where(HTime[Storms,4] >=1930.000171)[0]
    
    SNMin = np.where(SN <= 70)[0]
    SNMax = np.where(SN > 70)[0]
    SNYearMin = SNYears[SNMin]
    SNYearMax = SNYears[SNMax]

    #print(len(SNYearMin))
    #print(len(SNYearMax))
    #print(SNYearMin)
    #print(SNYearMax)

    fig, axs = plt.subplots(2,5,facecolor='w',edgecolor='k')
    fig.set_size_inches([6.69, 8.86], forward=True)
    fig.subplots_adjust(bottom=0.1,top=0.9,left=0.115,right=0.99,hspace=-0.8,wspace=-0.06)
    axs = axs.ravel()

    axID = 0
    lables=['C/SIR','ICME']
    colors=['midnightblue','maroon']
    #title =[r'-HMC$\in$[40,60)','[60,80)','[80,110)','[110,160)','[160,478]']
    title =[r'-HMC$\in$[16,40)','[40,50)','[50,65)','[65,90)','[90,478]']
    for G in range(5):
        Gindex = np.where(np.logical_and(HMC[Storms[Storms1932]]>HMCB[G+1],HMC[Storms[Storms1932]]<=HMCB[G]))[0]
        Gyears = HTime[Storms[Storms1932[Gindex]],0]
        Gmin = Gindex[np.where(np.in1d(Gyears,SNYearMin))[0]]
        Gmax = Gindex[np.where(np.in1d(Gyears,SNYearMax))[0]]

        MinCIR = np.where(Class[Storms1932[Gmin]] == 0)[0]
        MinCME = np.where(Class[Storms1932[Gmin]] == 1)[0]
        MaxCIR = np.where(Class[Storms1932[Gmax]] == 0)[0]
        MaxCME = np.where(Class[Storms1932[Gmax]] == 1)[0]
        Allmin = len(Gmin)
        Allmax = len(Gmax)

        fracsMin = [len(MinCIR)*100/len(Gmin),len(MinCME)*100/len(Gmin)]
        fracsMax = [len(MaxCIR)*100/len(Gmax),len(MaxCME)*100/len(Gmax)]

        if axID == 0:
            axs[axID].text(-0.05,1.05,title[axID],fontsize=12,transform=axs[axID].transAxes)
        else:
            axs[axID].text(0.3,1.05,title[axID],fontsize=12,transform=axs[axID].transAxes)
        axs[axID].text(0.4,0.92,str(len(Gmin)),fontsize=11,transform=axs[axID].transAxes)
        axs[axID].pie(fracsMin,colors=colors,autopct='%.0f%%',textprops={'size':'smaller','color':'white'})
        
        axs[axID+5].text(0.4,0.92,str(len(Gmax)),fontsize=11,transform=axs[axID+5].transAxes)
        axs[axID+5].pie(fracsMax,colors=colors,autopct='%.0f%%',textprops={'size':'smaller','color':'white'})

        if axID == 0:
            axs[axID].text(-0.05,0.75,'Solar Min.',fontsize=11,rotation=90,transform=axs[axID].transAxes)
            axs[axID+5].plot([0,0.24],[0,0],color='midnightblue',transform=axs[axID+5].transAxes)
            axs[axID+5].plot([0.6,0.84],[0,0],color='maroon',transform=axs[axID+5].transAxes)
            axs[axID+5].text(0.27,-0.025,'C/SIR',fontsize=11,transform=axs[axID+5].transAxes)
            axs[axID+5].text(0.87,-0.025,'ICME',fontsize=11,transform=axs[axID+5].transAxes)
            axs[axID+5].text(-0.05,0.75,'Solar Max.',fontsize=11,rotation=90,transform=axs[axID+5].transAxes)
        
        axID += 1
    
    if Save == True:
        fig.savefig('./Dump/Fig/development/Statistics_'+SaveName+'.pdf',format='pdf',dpi=200,bbox_inches='tight',pad_inches=0.01)
        #fig.savefig('./Dump/Fig/Statistics.png',format='png',dpi=200,bbox_inches='tight',pad_inches=0.01)
    plt.show()
###
###