###MITRIE NOAMER 1400 PC1
#####
	#this scripts collates and emulates the calculation of the AD1400 PCs in the mitrie data set mitrie_new_proxy_pcs_1400_v01.nc
	#the results can be replicated up to sign and scale (4 9s correlation) for 20 of 30 series
	#emulated are all the 70-series NOAMER PCs and all the SWM PCs
	#not replicated are the 56-series subset. 

#Fig. 1. Proxy principal components: the rst principal component of the North American ITRDB network
#of Mann et al., 1998. (1) Using the normalisation as in Mann et al. 1998, (2) as (1), but using full variance
#for normalisation rather than detrended variance, (3) normalised and centred on the whole series, (4)
#centred only (5) as archived by MBH1998. 21-year running means.


#`mbh98' is the collection of proxies used by Mann et al., (1998) which extend back to AD1400. `mbhx' and `mbh98x' are variations of `mbh' and
#`mbh98' respectively in which the proxy principal components have been recalculated.


#### MITRIE CODING
	#mbh: using the standardisation of mbh: centred and standardised on detrended-variance of last 79 years.
	#mbhx: centred and standardised on variance of last 79 years
	#mbhl: centred and standardised on detrended-variance of last 125 years
	#cen: centred on the whole period.
	#std: centred and normalised on the whole period.


######
 #COLLATE MBH1400 PC1s AS SAVED BY MITRIE
#######

	##################
	#### N. B. #####
	#the zip file has been previous downloaded and unzipped into a directory "C:/climate/data/mitrie"
	############
	#################

	library(ncdf)
	url<- "c:/climate/data/mitrie"
	loc<-file.path(url, "mitrie_new_proxy_pcs_1400_v01.nc"   ) #[1] "c:/climate/data/mitrie/mitrie_new_reconstructions_1400_v01.nc"
	fred<-open.ncdf(loc)
	mbh.pc1400<-array(NA,dim=c(581,30))
	for (i in 1:30) {
		v1 <- fred$var[[i]]
		mbh.pc1400[,i] <- get.var.ncdf( fred, v1 ) # by default, reads ALL the data
		}

 	col0<-c("black","red","blue")


###PC COLLECTION #this is MITIRE nomenclature # see coding described above
# [1] "mppc01_itrdb_namer_1400_mbh_00"  "mppc02_itrdb_namer_1400_mbh_00"  "mppc01_itrdb_namer_1400_mbhl_00"
# [4] "mppc02_itrdb_namer_1400_mbhl_00" "mppc01_itrdb_namer_1400_mbhx_00" "mppc02_itrdb_namer_1400_mbhx_00"
# [7] "mppc01_itrdb_namer_1400_std_00"  "mppc02_itrdb_namer_1400_std_00"  "mppc01_itrdb_namer_1400_cen_00" 
#[10] "mppc02_itrdb_namer_1400_cen_00"  "mppc01_stahle_swm_1400_mbh_00"   "mppc01_stahle_swm_1400_mbhl_00" 
#[13] "mppc01_stahle_swm_1400_mbhx_00"  "mppc01_stahle_swm_1400_std_00"   "mppc01_stahle_swm_1400_cen_00"  
#[16] "mppc01_itrdb_namer_1400_mbh_01"  "mppc02_itrdb_namer_1400_mbh_01"  "mppc01_itrdb_namer_1400_mbhl_01"
#[19] "mppc02_itrdb_namer_1400_mbhl_01" "mppc01_itrdb_namer_1400_mbhx_01" "mppc02_itrdb_namer_1400_mbhx_01"
#[22] "mppc01_itrdb_namer_1400_std_01"  "mppc02_itrdb_namer_1400_std_01"  "mppc01_itrdb_namer_1400_cen_01" 
#[25] "mppc02_itrdb_namer_1400_cen_01"  "mppc01_stahle_swm_1400_mbh_01"   "mppc01_stahle_swm_1400_mbhl_01" 
#[28] "mppc01_stahle_swm_1400_mbhx_01"  "mppc01_stahle_swm_1400_std_01"   "mppc01_stahle_swm_1400_cen_01"  


#######
  ##COLLATE EMULATIONS OF MBH1400 PCs
########

      #SET UP PLACEHOLDER ARRAYS
	emulation<-array(NA,dim=c(581,30))
	eigen.NOAMER<-array(NA, dim=c(70,30))
	eigen.SWM<-array(NA,dim=c(6,30))


###########
# NOAMER NETWORK #SERIES 1-10; 16-25 #SERIES 1-10 are 56-series network; #16-25 are 70-series network
###########


  #LOAD NOAMER NETWORK, EXTEND BY PERSISTENCE AND TRUNCATE TO 1980 
	#url.source<-"ftp://ftp.agu.org/apend/gl/2004GL021750"
	url.source<-"c:/climate/scripts/2004GL021750"
	loc<-file.path(url.source,"2004GL021750-NOAMER.1400.txt")
	tree<-read.table(loc,header=TRUE,sep="\t") #collated AD1400 network
	tree<-read.table("c:/climate/data/mann/UVA/NOAMER.txt",header=TRUE)#593 213

	tree<-ts(tree[,2:ncol(tree)],start=tree[1,1])
	temp<-!is.na(tree[1,])
	tree<-tree[,temp]
	tree<-ts(tree[1:581,],start=1400)
	temp81<-!is.na(tree[581,]);sum(temp81)#56

	tree<-ts( tree[1:581,] ,start=1400)	
	tree<-extend.persist(tree)  
	dim(tree) #581 70

  ##COLS 1-2 mbh_00	
	tree.mannomatic<-mannomatic(tree)
	pca.mannomatic<-prcomp(tree.mannomatic[,temp81],center=FALSE)
	emulation[,1:2]<-pca.mannomatic$x[,1:2]
	c(cor(emulation[,1],mbh.pc1400[,1]), cor(emulation[,2],mbh.pc1400[,2]))# [1] -0.999995  0.999994
	eigen.NOAMER[temp81,1:2]<-pca.mannomatic$rotation[,1:2]

   ##cols 16-17 mbh_01
	pca.mannomatic<-prcomp(tree.mannomatic,center=FALSE)
	emulation[,16:17]<-pca.mannomatic$x[,1:2]
	c(cor(emulation[,16],mbh.pc1400[,16]), cor(emulation[,17],mbh.pc1400[,17]))# [1] -0.9999953 -0.9999964
	eigen.NOAMER[,16:17]<-pca.mannomatic$rotation[,1:2]

   ##cols 3-4:  #mbhl: centred and standardised on detrended-variance of last 125 years
	tree.mannomatic<-mannomatic(tree,M=124) 
	pca<- prcomp(tree.mannomatic[,temp81],center=FALSE) 
	emulation[,3:4]<-pca$x[,1:2]
	c(cor(emulation[,3],mbh.pc1400[,3]), cor(emulation[,4],mbh.pc1400[,4]))# [1] -0.9999960  0.9999941
	eigen.NOAMER[temp81,3:4]<-pca$rotation[,1:2]

   ##cols 18-19:  #mbhl: centred and standardised on detrended-variance of last 125 years
	pca<- prcomp(tree.mannomatic,center=FALSE) 
	emulation[,18:19]<-pca$x[,1:2]
	c(cor(emulation[,18],mbh.pc1400[,18]), cor(emulation[,19],mbh.pc1400[,19]))# [1]  -0.9999962 -0.9999963
	eigen.NOAMER[,18:19]<-pca$rotation[,1:2]

  ##cols 5-6 #mbhx: centred and standardised on variance of last 79 years
	tree.mannomatic<-mannomatic(tree,sdmethod="undetrended") 
	pca<-prcomp(tree.mannomatic[,temp81],center=FALSE) 
	emulation[,5:6]<-pca$x[,1:2]
	c(cor(emulation[,5],mbh.pc1400[,5]), cor(emulation[,6],mbh.pc1400[,6]))# [1]-0.9999952  0.9999930
	eigen.NOAMER[temp81,5:6]<-pca$rotation[,1:2]

  ##cols 20-21:  # #mbhx: centred and standardised on variance of last 79 years
	pca<- prcomp(tree.mannomatic,center=FALSE) 
	emulation[,20:21]<-pca$x[,1:2]
	c(cor(emulation[,20],mbh.pc1400[,20]), cor(emulation[,21],mbh.pc1400[,21]))# [1]  -0.9999955 -0.9999963
	eigen.NOAMER[,20:21]<-pca$rotation[,1:2]

  ###cols 7-8 std: centred and normalised on the whole period.
	tree.scaled<-scale(tree)
	pca<- prcomp(tree.scaled[,temp81]) 
	emulation[,7:8]<-pca$x[,1:2]
  	c(cor(emulation[,7],mbh.pc1400[,7]), cor(emulation[,8],mbh.pc1400[,8]))# [1]  0.999996 -0.999993
	eigen.NOAMER[temp81,7:8]<-pca$rotation[,1:2]
	
  ###cols 22-23 std: centred and normalised on the whole period.
	pca<- prcomp(tree.scaled) 
	emulation[,22:23]<-pca$x[,1:2]
	c(cor(emulation[,22],mbh.pc1400[,22]), cor(emulation[,23],mbh.pc1400[,23]))# [1] -0.9999965 -0.9999937
	eigen.NOAMER[,22:23]<-pca$rotation[,1:2]

  ###cols 9-10 COVARIANCE #cen: centred on the whole period.
	tree.centered<-scale(tree,scale=FALSE)
	pca<-prcomp(tree.centered[,temp81]) 
	emulation[,9:10]<-pca$x[,1:2]
 	c(cor(emulation[,9],mbh.pc1400[,9]), cor(emulation[,10],mbh.pc1400[,10]))# [1]   -0.9999971  0.9999956
	eigen.NOAMER[temp81,9:10]<-pca$rotation[,1:2]

  #cols 24-25 std: centered on whole period full ser 
  	pca<-prcomp(tree.centered) 
	emulation[,24:25]<-pca$x[,1:2]
 	c(cor(emulation[,24],mbh.pc1400[,24]), cor(emulation[,25],mbh.pc1400[,25]))# [1]  -0.9999985  0.9999967
	eigen.NOAMER[,24:25]<-pca$rotation[,1:2]

################
##STAHLE SWM NETWORK  - 11-15: 2 series; 26-30 - 6 series
##################
	loc<-"c:/climate/data/mann/UVA/STAHLE.SWM.txt"
	tree<-read.table(loc,header=TRUE,sep="\t") #starts 1375
	tree<-ts(tree[(1400:1980)-1374,2:ncol(tree)],start=1400)
		dim(tree)
	temp<-!is.na(tree[1,])
	tree<-tree[,temp];dim(tree)
	temp81<-!is.na(tree[581,]);sum(temp81) #2
	id<-dimnames(tree)[[2]][!temp81] ;id<-strip(id,2)
		# "swmxdfew01" "swmxdfew11" "swmxdflw01" "swmxdflw11"
	round(apply(tree[,!temp81],2,hockeystat),2)
	tree<-extend.persist(tree)
	tree.mannomatic<-mannomatic(tree) 
	
   ##COLS 11 mbh_00	
	pca.mannomatic<-prcomp(tree.mannomatic[,temp81],center=FALSE)
	emulation[,11]<-pca.mannomatic$x[,1]
	cor(emulation[,11],mbh.pc1400[,11])# [1]  -0.9999862
	eigen.SWM[temp81,11]<-pca.mannomatic$rotation[,1]

   ##cols 26 mbh_01
	pca.mannomatic<-prcomp(tree.mannomatic,center=FALSE)
	emulation[,26]<-pca.mannomatic$x[,1]
	cor(emulation[,26],mbh.pc1400[,26])# [-0.9999912
	eigen.SWM[,26]<-pca.mannomatic$rotation[,1]

   ##cols 12:  #mbhl: centred and standardised on detrended-variance of last 125 years
	tree.mannomatic<-mannomatic(tree,M=124) 
	pca<- prcomp(tree.mannomatic[,temp81],center=FALSE) 
	emulation[,12]<-pca$x[,1]
	cor(emulation[,12],mbh.pc1400[,12])#  -0.999986
	eigen.SWM[temp81,12]<-pca$rotation[,1]

   ##cols 27:  #mbhl: centred and standardised on detrended-variance of last 125 years
	pca<- prcomp(tree.mannomatic,center=FALSE) 
	emulation[,27]<-pca$x[,1]
	cor(emulation[,27],mbh.pc1400[,27])# [1] -0.9999902
	eigen.SWM[,27]<-pca$rotation[,1]
	
  ##cols 13 #mbhx: centred and standardised on variance of last 79 years
	tree.mannomatic<-mannomatic(tree,sdmethod="undetrended") 
	pca<-prcomp(tree.mannomatic[,temp81],center=FALSE) 
	emulation[,13]<-pca$x[,1]
	cor(emulation[,13],mbh.pc1400[,13])# [1] -0.9999862
	eigen.SWM[temp81,13]<-pca$rotation[,1]

  ##cols 28:  # #mbhx: centred and standardised on variance of last 79 years
	pca<- prcomp(tree.mannomatic,center=FALSE) 
	emulation[,28]<-pca$x[,1]
	cor(emulation[,28],mbh.pc1400[,28])# [1]  -0.9999915
	eigen.SWM[,28]<-pca$rotation[,1]

  ###cols 14 std: centred and normalised on the whole period.
	tree.scaled<-scale(tree)
	pca<- prcomp(tree.scaled[,temp81]) 
	emulation[,14]<-pca$x[,1]
  	cor(emulation[,14],mbh.pc1400[,14])# [1]  -0.9999875
	eigen.SWM[temp81,14]<-pca$rotation[,1]
	
  ###cols 29 std: centred and normalised on the whole period.
	pca<- prcomp(tree.scaled) 
	emulation[,29]<-pca$x[,1]
	cor(emulation[,29],mbh.pc1400[,29])# [1] -0.9999935
	eigen.SWM[,29]<-pca$rotation[,1]

  ###cols 15 COVARIANCE #cen: centred on the whole period.
	tree.centered<-scale(tree,scale=FALSE)
	pca<-prcomp(tree.centered[,temp81]) 
	emulation[,15]<-pca$x[,1]
 	cor(emulation[,15],mbh.pc1400[,15])# [1] -0.9999875
	eigen.SWM[,15]<-pca$rotation[,1]

  #cols 30 std: centered on whole period full ser 
  	pca<-prcomp(tree.centered) 
	emulation[,30]<-pca$x[,1]
 	cor(emulation[,30],mbh.pc1400[,30])# [1]-0.9999938
	eigen.SWM[,30]<-pca$rotation[,1]


###COMPARE TO BENCHMARK
	stat<-rep(NA,30)
	for (i in 1:30){ stat[i]<-cor(emulation[,i],mbh.pc1400[,i])}
		# [1] -0.9999951  0.9999939 -0.9999960  0.9999941 -0.9999952  0.9999930  0.9999961 -0.9999931 -0.9999971
		#[10]  0.9999956 -0.9999862 -0.9999860 -0.9999862 -0.9999875 -0.9999875 -0.9999953 -0.9999964 -0.9999962
		#[19] -0.9999963 -0.9999955 -0.9999963 -0.9999965 -0.9999937 -0.9999985  0.9999967 -0.9999912 -0.9999902
		#[28] -0.9999915 -0.9999935 -0.9999938


##NOAMER EIGENVECTOR COEFFICIENTS
	apply(eigen.NOAMER[,c(1:10,16:25)],2,sum,na.rm=T)
	#	 5.16945204  6.60617788  5.21884208 -1.69359497  5.31797049  6.44227423  8.36867447  0.09653672 -4.08336239
	#[10]  0.87856786  5.16197122  5.14073506  5.23599562  5.16053705  5.31599727  5.00335364  7.45450166  0.66516379
	#[19]  6.02031330  2.85946920		
		#eigenvector sums are positive

#SWM EIGENVECTOR COEFFICIENTS
	apply(eigen.SWM[,c(11:15,26:30)],2,sum,na.rm=T)
	#[1] 1.340917 1.348494 1.340307 1.414214 4.236397 2.410663 2.423051 2.412154 2.241449 2.207887
		#all are positive