Commit d57fe4d3 authored by Leonie Pick's avatar Leonie Pick

Test of event selection parameters

parent 8b98f721
This diff is collapsed.
......@@ -127,17 +127,17 @@ def Search_TargetEvents(HMC, HMC11y, HMC5d, dHMC, HTime, DTime, grid, Save):
if HMC_diff <= -thresHMC:
IndexMin2.append(i)
pl.Selection(HTime,DTime,HMC,Index_thres1,StormIndices,IndexMin,IndexMin1,IndexMin2,Save)
#pl.Selection(HTime,DTime,HMC,Index_thres1,StormIndices,IndexMin,IndexMin1,IndexMin2,Save)
return IndexMin2
###
###
def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Training, Time, Date, YearStart, Save):
def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Training, Time, Date, YearStart, Plot, Save):
YearsIndex = np.where(Time[:,0] >= YearStart)[0]
TrClass = Training[:,0]
TrTimeIndex = Training[:,1]
TrFound = Training[:,2]
#TrFound = Training[:,2]
CIRs = np.where(np.logical_or(TrClass == 0, TrClass == 2))[0]
CMEs = np.where(TrClass==1)[0]
......@@ -149,54 +149,68 @@ def Get_TargetEvents(HMC, HMC11y, HMC5d, dHMC, Kp_all, KpHours_all, Training, Ti
## Parameters
percentile = np.arange(29,29.1,1) #percentile = np.arange(8,9,1)
#percentile = np.arange(5,55,5)
power = np.arange(1,1.1,1) #power = np.arange(0.9,1.0,0.1)
mini_range = np.arange(26,26.1,1) #mini_range = np.arange(78,79,1)
#mini_range = np.arange(2,58,6)
thresHMC = np.arange(7,7.1,1) #thresHMC = np.arange(10,11,1)
#thresHMC = np.arange(0,20,2)
grid = np.stack(np.meshgrid(percentile,power,mini_range,thresHMC),-1).reshape(-1,4)
TrFound = np.zeros((len(grid),len(Training)))
StormsWin = dict()
FoundStorms = dict()
#print(grid)
threshold = 0
for i in range(len(grid)):
Storms = Search_TargetEvents(HMC[YearsIndex], HMC11y[YearsIndex], HMC5d[YearsIndex], dHMC[YearsIndex],
Time[YearsIndex,:], Date[YearsIndex], grid[i,:], Save)
Time[YearsIndex,:], Date[YearsIndex], grid[i,:], Plot)
Found_CIRs = np.where(np.in1d(TrTimeIndex[np.logical_or(TrClass==0,TrClass==2)],Storms+YearsIndex[0]))[0]
Found_CMEs = np.where(np.in1d(TrTimeIndex[TrClass==1],Storms+YearsIndex[0]))[0]
#Found_CIRs = np.where(np.in1d(TrTimeIndex[np.logical_or(TrClass==0,TrClass==2)],Storms+YearsIndex[0]))[0]
#Found_CMEs = np.where(np.in1d(TrTimeIndex[TrClass==1],Storms+YearsIndex[0]))[0]
FoundStorms['CIRs_'+str(int(i))] = np.where(np.in1d(TrTimeIndex[np.logical_or(TrClass==0,TrClass==2)],Storms+YearsIndex[0]))[0]
FoundStorms['CMEs_'+str(int(i))] = np.where(np.in1d(TrTimeIndex[TrClass==1],Storms+YearsIndex[0]))[0]
#IndexWin = i
#StormsWin = Storms+YearsIndex[0]
#CIRIndices = CIRs[Found_CIRs]
#CMEIndices = CMEs[Found_CMEs]
if ((len(Found_CIRs) >= np.floor((len(CIRs)-len(CIRs)*accepted_loss_CIRs/100))) and (len(Found_CMEs) >= np.floor((len(CMEs)-len(CMEs)*accepted_loss_CMEs/100)))):
criterium = (len(Found_CIRs)+len(Found_CMEs))+(len(YearsIndex)-len(Storms))
if criterium > threshold:
IndexWin = i
StormsWin = Storms+YearsIndex[0]
CIRIndices = CIRs[Found_CIRs]
CMEIndices = CMEs[Found_CMEs]
threshold = criterium
StormsWin[str(int(i))] = Storms+YearsIndex[0]
CIRIndices = CIRs[FoundStorms['CIRs_'+str(int(i))]]
CMEIndices = CMEs[FoundStorms['CMEs_'+str(int(i))]]
#if ((len(Found_CIRs) >= np.floor((len(CIRs)-len(CIRs)*accepted_loss_CIRs/100))) and (len(Found_CMEs) >= np.floor((len(CMEs)-len(CMEs)*accepted_loss_CMEs/100)))):
#criterium = (len(Found_CIRs)+len(Found_CMEs))+(len(YearsIndex)-len(Storms))
#if criterium > threshold:
#IndexWin = i
#StormsWin = Storms+YearsIndex[0]
#CIRIndices = CIRs[Found_CIRs]
#CMEIndices = CMEs[Found_CMEs]
#threshold = criterium
#if Save == True:
#np.save('./Dump/Out/development/Selection_P'+str(int(grid[i,0])),StormsWin)
TrFound[CIRIndices] = 1
TrFound[CMEIndices] = 1
Training[:,2] = TrFound
#print(grid[IndexWin,:])
#if Save == True:
#np.save('./Dump/Events',StormsWin)
TrFound[i,CIRIndices] = 1
TrFound[i,CMEIndices] = 1
#Training[:,2] = TrFound
##print(grid[IndexWin,:])
IndexWin = 0
SelectStr = str(int(IndexWin))
Training[:,2] = TrFound[IndexWin,:]
###
TargetResults = {'No. targets':len(StormsWin),'Fraction [%]':len(StormsWin)*100/len(HMC[YearsIndex])}
TrainingResults = [[len(CIRs), len(CMEs), len(CIRs)+len(CMEs)],[len(Found_CIRs),len(Found_CMEs),len(Found_CIRs)+len(Found_CMEs)]]
SelectionParams = {'HMC percentile':grid[IndexWin,0],'Scaling power':grid[IndexWin,1], 'Separation [h]':grid[IndexWin,2], 'Min. HMC drop [nT]':grid[IndexWin,3]}
display(pd.DataFrame(data=SelectionParams,index=['Selection parameters']))
display(pd.DataFrame(data=TargetResults,index=['Selection result']))
TargetResults = {'No. targets': len(StormsWin[SelectStr]),'Fraction [%]':len(StormsWin[SelectStr])*100/len(HMC[YearsIndex])}
TrainingResults = [[len(CIRs), len(CMEs), len(CIRs)+len(CMEs)],[len(FoundStorms['CIRs_'+SelectStr]),len(FoundStorms['CMEs_'+SelectStr]),len(FoundStorms['CIRs_'+SelectStr])+len(FoundStorms['CMEs_'+SelectStr])]]
SelectionParams = {'HMC percentile':grid[:,0],'Scaling power':grid[:,1], 'Separation [h]':grid[:,2], 'Min. HMC drop [nT]':grid[:,3]}
display(pd.DataFrame(data=SelectionParams))
display(pd.DataFrame(data=TargetResults,index=['Selection']))
display(pd.DataFrame(data=TrainingResults,columns=['No. CIRs','No. CMEs', 'Total'], index=['Training set', 'in Target set']))
###
pl.IndexDist(Time,YearsIndex,StormsWin,Kp_all,KpHours_all,HMC,Save)
#pl.IndexDist(Time,YearsIndex,StormsWin,Kp_all,KpHours_all,HMC,Plot)
return StormsWin
return StormsWin, TrFound, FoundStorms
###
###
def Get_F0(i,HMC,HMC_filt1):
......@@ -401,7 +415,7 @@ def Get_Features(Time, Storms, Training, HMC, HMC11y, HMC1y, dHMC, B_MLT, ASY, S
print('Caution!')
if Save == True:
np.savez('./Dump/Out/Features',Features=FeatureMatrix,FeatureIndices=PeakIndices)
np.savez('./Dump/Out/development/Features_SC0',Features=FeatureMatrix,FeatureIndices=PeakIndices)
return PeakIndices, FeatureMatrix
###
......@@ -526,7 +540,7 @@ def Get_Pipelines(whichEstimator,nClasses):
return BPipe,BEstimateParams,OPipe,SelectParams,OEstimateParams
###
###
def Select_Features(Data,Target,N2,K2,Pipe,SelectParams,EstimateParams,scoring,scoring_name,n_classes,estimators,est,ax,Save):
def Select_Features(Data,Target,N2,K2,Pipe,SelectParams,EstimateParams,scoring,scoring_name,n_classes,estimators,est,ax,Save,SaveName):
nFeatures = Data.shape[1]
Scores = np.zeros((N2,nFeatures,2))
......@@ -572,7 +586,7 @@ def Select_Features(Data,Target,N2,K2,Pipe,SelectParams,EstimateParams,scoring,s
ranking = np.argsort(transformer.ranking_)
## Plot test results
ax = pl.Ranking(Scores_mean, Scores_std, Data, ranking, FeaturesSelect, estimators, scoring_name, est, ax, Save)
ax = pl.Ranking(Scores_mean, Scores_std, Data, ranking, FeaturesSelect, estimators, scoring_name, est, ax, Save, SaveName)
return EstimateParams, ax
###
......@@ -618,7 +632,7 @@ def Optimize_Model(Data,Target,N1,K1,model,params,nModels,scoring,scoring_name):
return CV1_counts, CV1_means, Param_Comb
###
###
def Find_Model(Data,Target,N1,K1,N2,K2,model,params,scoring,scoring_name,n_classes,Save):
def Find_Model(Data,Target,N1,K1,N2,K2,model,params,scoring,scoring_name,n_classes,Save,SaveName):
nModels = 1
for p in params.keys(): nModels *= len(params[p])
......@@ -662,7 +676,7 @@ def Find_Model(Data,Target,N1,K1,N2,K2,model,params,scoring,scoring_name,n_class
display(pd.DataFrame(data=Hyperparams,index=['Best model'],columns=list(ModelB.keys())))
###
pl.ModelCounts(nModels,SumCounts,CV2_means,N1,K1,N2,K2,Param_Comb,BestModelTestB,scoring_name,Save)
pl.ModelCounts(nModels,SumCounts,CV2_means,N1,K1,N2,K2,Param_Comb,BestModelTestB,scoring_name,Save,SaveName)
return Param_Comb[BestModelTestB]
###
......@@ -697,7 +711,7 @@ def Get_MyScores(CM_means, n_classes):
return C, Scores
###
###
def Assess_Model(Data,Target,N2,K2,modelTest,n_classes,Plot,Save):
def Assess_Model(Data,Target,N2,K2,modelTest,n_classes,Plot,Save,SaveName):
## ROC, Precision-Recall curves, confusion matrix
curve_i = np.linspace(0,1,100)
......@@ -764,8 +778,8 @@ def Assess_Model(Data,Target,N2,K2,modelTest,n_classes,Plot,Save):
# Plot mean confusion matrix and curves
if Plot == True:
pl.CM(CM_means[:,0],C,n_classes,Save)
pl.Curves(N2,K2,Curves,curve_i,Model_Mean[0:2,0],Model_Std[0:2,0],C,Save)
pl.CM(CM_means[:,0],C,n_classes,Save,SaveName)
pl.Curves(N2,K2,Curves,curve_i,Model_Mean[0:2,0],Model_Std[0:2,0],C,Save,SaveName)
return Model_Mean, Model_Std
###
......
......@@ -68,8 +68,8 @@ def Selection(Time,Date,HMC,Index_thres1,StormIndices,IndexMin,IndexMin1,IndexMi
ax[i].tick_params(axis='x',which='minor',direction='in')
if Save == True:
fig.savefig('./Dump/Fig/EventSelection.pdf',format='pdf',dpi=200,transparent=True)
#fig.savefig('./Dump/Fig/EventSelection.png',format='png',dpi=200,transparent=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()
###
###
......@@ -190,8 +190,8 @@ def IndexDist(Time,YearsIndex,Storms,Kp_all,KpHours_all,HMC,Save):
#ax3.yaxis.label.set_color('maroon')
if Save == True:
fig.savefig('./Dump/Fig/HMCDist.pdf',format='pdf',dpi=200,transparent=True)
#fig.savefig('./Dump/Fig/HMCDist.png',format='png',dpi=200,transparent=True)
fig.savefig('./Dump/Fig/development/HMCDist.pdf',format='pdf',dpi=200,transparent=True)
#fig.savefig('./Dump/Fig/development/HMCDist.png',format='png',dpi=200,transparent=True)
plt.show()
###
###
......@@ -243,12 +243,12 @@ def Diagnostics(n_features,n_classes,NData, Save):
k += 1
if Save == True:
fig.savefig('./Dump/Fig/Features.pdf',format='pdf',dpi=200,transparent=True)
fig.savefig('./Dump/Fig/development/Features_SC0.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):
def Ranking(Scores, Std, Data, Ranking, nfeaturesFinal, estimators, scoring_name, est, ax, Save, SaveName):
nEstimators = len(estimators)
......@@ -287,14 +287,14 @@ def Ranking(Scores, Std, Data, Ranking, nfeaturesFinal, estimators, scoring_name
if est == nEstimators-1:
if Save == True:
fig.savefig('./Dump/Fig/Ranking.pdf',format='pdf',dpi=300,transparent=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):
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)
......@@ -333,29 +333,30 @@ def ModelCounts(nModels,SumCounts,CV2_means,N1,K1,N2,K2,Param_Comb,BestID,scorin
k += 1
if Save == True:
fig.savefig('./Dump/Fig/ModelStats.pdf',format='pdf',dpi=300,transparent=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):
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 = ['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)
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))
##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)
......@@ -381,13 +382,13 @@ def CM(CM,Condition,n_classes,Save):
ax.grid(False)
if Save == True:
fig.savefig('./Dump/Fig/ConfusionCol.pdf',format='pdf',dpi=300,transparent=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):
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)
......@@ -444,7 +445,7 @@ def Curves(N2,K2,Curves,curve_i,Model_Mean,Model_Std,C,Save):
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/Curves.pdf',format='pdf',dpi=300,transparent=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()
......@@ -634,12 +635,12 @@ def Musical(Time,HTime,Storms,Kp_all,KpHours_all,SN,SNYears,HMC,HMC_filt,Trainin
n, bins, patches = axs2.hist(x=HTime[Storms,4],bins=np.arange(1900,2017,1),orientation='horizontal',histtype='stepfilled',color='silver',alpha=1)
if Save == True:
fig.savefig('./Dump/Fig/StormsSetup.pdf',format='pdf',dpi=200)#,transparent=True)
#fig.savefig('./Dump/Fig/StormsSetup.png',format='png',dpi=200)#,transparent=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):
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)
......@@ -788,13 +789,13 @@ def MusicalClassified(Time,HTime,Storms,Class_Proba,Kp_all,KpHours_all,SN,SNYear
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/StormsClassified.pdf',format='pdf',dpi=200)#,transparent=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):
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]
......@@ -861,7 +862,7 @@ def Stats(HTime,Storms,SN,SNYears,HMC,Proba,Class,Save):
axID += 1
if Save == True:
fig.savefig('./Dump/Fig/Statistics.pdf',format='pdf',dpi=200,bbox_inches='tight',pad_inches=0.01)
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()
###
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment