# "R" SCRIPT TO DO RE STATISTIC THRESHOLD ANALYSES -- BASED ON MONTE CARLO ANALYSES USING FULL-SPECTRUM RED NOISE PSEUDO
# PROXIES (ALL PROXIES) FOR MBH/WA AND WA SCENARIOS 5 AND 6, INPUT INTO WA EMULATION OF MBH RECONSTRUCTION ALGORITHM.


# FILE_PATH:  In command line below, specify file path as directed.

# FILE NAMES have a specific structure that is set according to the calling and looping done in this script.
# ...Either specify the Validation.stats output files with the names as set herein, or give them different names that
#      can still conform to the script as used here.
# ...Do the same for the Verif_RE_tests script files that are called for Wahl-Ammann Scenarios 5 and 6.


rm(list=ls())


FILE_PATH<-paste("SPECIFY_FULL_PATH_TO_APPROPRIATE_Validation.stats_AND_Verif_RE_tests_FILES_WITHIN_QUOTATION_MARKS")


#************************** SECTION TO DO MAIN MBH/WA RECONSTRUCTIONS ***************************


Network<-array(dim=c(1,12))

Network[1,1]<-1000
Network[1,2]<-1400
Network[1,3]<-1450
Network[1,4]<-1500
Network[1,5]<-1600
Network[1,6]<-1700
Network[1,7]<-1730
Network[1,8]<-1750
Network[1,9]<-1760
Network[1,10]<-1780
Network[1,11]<-1800
Network[1,12]<-1820


WA_verifRE<-array(dim=c(1,12))

WA_verifRE[1,1]<-0.2215955 #full-time-period centering AND no CO2 correction
WA_verifRE[1,2]<-0.4811857
WA_verifRE[1,3]<-0.4391771
WA_verifRE[1,4]<-0.4651978
WA_verifRE[1,5]<-0.4626561
WA_verifRE[1,6]<-0.5045778
WA_verifRE[1,7]<-0.5462276
WA_verifRE[1,8]<-0.6058250
WA_verifRE[1,9]<-0.5431693
WA_verifRE[1,10]<-0.5902400
WA_verifRE[1,11]<-0.6025641
WA_verifRE[1,12]<-0.6169262


#WA_calRE_1400<-0.3948435
#mean_ratio_WA_calRE_to_WA_verif_RE<-1.196  # for 1400-1820 segments



MC_reps<-1000


Validation.stats_cols<-2



print(paste("******************* MAIN MBH/WA RECONSTRUCTIONS *********************"))
print(paste(""))
print(paste(""))



for (h in 1:12) {


if (h==1) {
print(paste("****** MBH/WA ",Network[1,h],"-NETWORK__no CO2 correction/full-period-centered proxy PCs ******",sep=""))}
if (h!=1) {
print(paste("******************* MBH/WA ",Network[1,h],"-NETWORK *********************",sep=""))}
print(paste(""))
print(paste(""))


if (h==1) {
	Validation.stats<-scan(paste(FILE_PATH,"Validation.stats_MBH_RE_simulations",Network[1,h],"_fullCentering_noCO2.dat",sep=""))}

if (h!=1) {
	Validation.stats<-scan(paste(FILE_PATH,"Validation.stats_MBH_RE_simulations",Network[1,h],".dat",sep=""))}

	Validation.stats<-matrix(Validation.stats,nrow=MC_reps,ncol=Validation.stats_cols, byrow=TRUE)


for (hh in c(0.5,0.6666667,0.75,0.9,1.0,1.25,1.5)) {

ratio_cal_verif<-hh


	test_array<-array(0,dim=c(MC_reps,2))

		for (k in 1:MC_reps) {


		if (Validation.stats[k,1]>0) {

			if (Validation.stats[k,1]/Validation.stats[k,2]>=ratio_cal_verif) {

				if (Validation.stats[k,2]>=WA_verifRE[1,h]) {

				test_array[k,1]<-1 } #close third if test

			 } #close second if test

		} #close first if test


		test_array[k,2]<-Validation.stats[k,2]

		if (Validation.stats[k,1]<=0) {

			test_array[k,2]<-(-9999) } else {

			if (Validation.stats[k,1]/Validation.stats[k,2]<ratio_cal_verif) {

			test_array[k,2]<-(-9999) } #close fifth if test

		} #close fourth if test

		} #close k loop


	z<-sum(test_array[,1])

	test_statistic<-(1-(z/MC_reps))


	print(paste("test_statistic",Network[1,h]," calib/verif RE ratio ",hh," = ",test_statistic,sep=""))

	quantiles<-quantile(test_array[,2],c(.5,.8,.9,.95,.96,.97,.98,.985,.99))

	print(paste("quantiles",Network[1,h]," calib/verif RE ratio ",hh,sep=""))

	print(quantiles)
	print(paste(""))
	print(paste(""))


	rm(z)

	rm(test_statistic)

	rm(quantiles)


} #close hh loop

	rm(Validation.stats)


} #close h loop


print(paste(""))
print(paste(""))



#******************************** SECTION TO DO WA SCENARIO 5a/b/c/d RECONSTRUCTIONS ***********************


print(paste("********************* WA SCENARIO 5a/b/c/d RECONSTRUCTIONS *********************"))



source(paste(FILE_PATH,"WA-MBH_SCENARIOS_5a-d_Verif_RE_tests.txt",sep=""))



source(paste(FILE_PATH,"WA-MBH_SCENARIOS_5d_Verif_RE_tests_ALTERNATE_FOR_NEG-Verif-RE.txt",sep=""))



print(paste(""))
print(paste(""))



#******************************** SECTION TO DO WA SCENARIO 6a/b/c RECONSTRUCTIONS ***********************



print(paste("*********************** WA SCENARIO 6a/b/c RECONSTRUCTIONS *********************"))



source(paste(FILE_PATH,"WA-MBH_SCENARIOS_6a-c_Verif_RE_tests.txt",sep=""))



source(paste(FILE_PATH,"WA-MBH_SCENARIOS_6a-c_Verif_RE_tests_ALTERNATE_FOR_NEG-Verif-RE.txt",sep=""))



print(paste(""))
print(paste(""))



#**********************************************************************************************************



# END