##MCINTYRE AND MCKITRICK [EE 2005] FUNCTIONS

#TO EXTEND SERIES BY PERSISTENCE
	extend.persist<-function(tree) {
	extend.persist<-tree
	for (j in 1:ncol(tree) )  {
		test<-is.na(tree[,j])
		end1<-max ( c(1:nrow(tree)) [!test])
		test2<-( c(1:nrow(tree))>end1) & test
		extend.persist[test2,j]<-tree[end1,j]
		}
	extend.persist
	}
	
#JONES GRIDCELL ID FROM LAT AND LONG
	jones<-function(lat,long) {
	i<-18-floor(lat/5);
	j<- 36+ceiling(long/5);
	jones<-72*(i-1) + j
	jones
	}


#GAUSSIAN FILTER
	#see http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf

	gaussian.filter.weights<-function(N,year) {
		#N number of points; year - is number of years
	sigma<-year/6
	i<-( (-(N-1)/2):((N-1)/2)   )  /sigma
	w<- ((2*pi)^-0.5) * exp(-0.5* i^2)
	gaussian.filter.weights<-w/sum(w)
	gaussian.filter.weights
	}

	truncated.gauss.weights<-function(year) {
	a<-gaussian.filter.weights(2*year+1,year)
	temp<-(a>0.05*max(a))
	a<-a[temp]
	a<-a/sum(a)
	truncated.gauss.weights<-a
	truncated.gauss.weights
	}

  filter.combine.pad<-function(x,a,M=NA) {
	temp<-!is.na(x)  
   	year<-c(tsp(x)[1]:tsp(x)[2])[temp]
	w<-x[temp]
	N<-length(w)
	if(is.na(M)) M<-trunc (length(a)/2)+1
	w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M))
   	y<-filter(w,a,method="convolution",sides=2)
	y<-y[(M+1):(length(y)-M)]
	z<-ts (rep(NA,length(x)),start=tsp(x)[1],end=tsp(x)[2])
	z[temp]<-y
   	filter.combine.pad<-ts.union(x,z)
   	#filter.combine<-y[,2]
	filter.combine.pad
   	}

#MANNOMATIC DATA TRANSFORMATINO
   sd.detrend<-function(x) {t<-c(1:length(x));fm<-lm(x~t);sd.detrend<-sd(fm$residuals);sd.detrend}

   mannomatic<-function(x,M=78) {
	N<-length(x)
	xstd<- (x- mean( x[(N-M):N]))/sd(x[(N-M):N])
	sdprox<-sd.detrend(xstd[(N-M):N]) 
	mannomatic<- xstd/sdprox
	mannomatic
	}


#4. FUNCTION TO RETURN VERIFICATION PERIOD STATISTICS
	#calibration and verificaiton are both vectors of length 2 with start and end period
	#estimator and observed are time series which may have different start and end periods

  verification.stats<-function(estimator,observed,calibration,verification)  {
        combine<-ts.union(estimator,observed)
	index.cal<- (calibration[1]-tsp(combine)[1]+1):(calibration[2]-tsp(combine)[1]+1)       
	index.ver<-(verification[1]-tsp(combine)[1]+1):(verification[2]-tsp(combine)[1]+1)
	xmean<-mean ( combine[index.cal,"observed"],na.rm=TRUE  )
       	mx<-mean(  combine[index.ver,"observed"],na.rm=TRUE  ) # 
	test.ver<-cov ( combine[index.ver,],use=use0)
	test<-cov(combine[index.cal,],use=use0)
	R2.cal<-(test[1,2]*test[2,1])/(test[1,1]*test[2,2])
	R2.ver<-(test.ver[1,2]*test.ver[2,1])/(test.ver[1,1]*test.ver[2,2])
	RE.cal<-1-sum( (combine[index.cal,"estimator"]- combine[index.cal,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE)
	RE.ver<-1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	RE.risk<- - sum( combine[index.ver,"estimator"]^2,na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE)
	RE.bias<- 2*length(index.ver) * mean(combine[index.ver,"estimator"],na.rm=TRUE)*mean(combine[index.ver,"observed"],na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE)
	RE.covar<- 2* (length(index.ver)-1)* cov(combine[index.ver,],use="pairwise.complete.obs")[1,2]/sum( (xmean- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	CE<- 1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mx- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	test<-sign(diff(combine[index.ver,],lag=1))
	temp<-	( test[,1]*test[,2] ==1)&!is.na(test[,1])&!is.na(test[,2])
	sign.test<-sum(temp)/length(index.ver)
	test1<-scale(combine[index.ver,],scale=FALSE)
	test1<-test1[!is.na(test1[,1])&!is.na(test1[,2]),]
	temp<-sign(test1)
	temp<-(temp[,1]*temp[,2]==1)
	test2<-test1[,1]*test1[,2]
	if(sum(!temp) >0) prod.test<-  ( mean(test2[temp])+mean(test2[!temp]))/  sqrt( var(test2[temp])/sum(temp)  + var(test2[!temp])/sum(!temp) ) else prod.test<-NA
	verification.stats<-data.frame(RE.cal,RE.ver,R2.cal,R2.ver,CE,sign.test,prod.test,RE.risk,RE.bias,RE.covar)
	names(verification.stats)<-c( "RE.cal","RE.ver","R2.cal","R2.ver","CE","sign.test","prod.test","RE.risk","RE.bias","RE.covar")
	verification.stats
	}


     hockeystat<-function(test,N=502) { hockeystat<-  ( mean(test[N:(N+79)],na.rm=TRUE)- mean(test[1:(N+79)],na.rm=TRUE) )/sd(test[1:(N+79)],na.rm=TRUE); hockeystat }
	#mean of "blade" - overall mean/sd 

     hockeystat3<-function(test) { hockeystat3<-  ( mean(test[502:581],na.rm=TRUE)- mean(test[1:51],na.rm=TRUE) )/sd(test[1:581],na.rm=TRUE); hockeystat3 }
	#compare 20th to early 15th
 
     hockeystat4<-function(test) { hockeystat4<-  ( mean(test[532:581],na.rm=TRUE)- mean(test[1:51],na.rm=TRUE) )/sd(test[1:581],na.rm=TRUE); hockeystat4 }
	#compare last 50 years to first 50 years


#  UTILITY
	#Create a directory and subdirectories if they do not exist.
	createsubdirs<-function(dirpath1) {
	if( ! file.exists(dirpath1) ) {
		if( ! file.exists(dirname(dirpath1)) ) {
			# If parent directory is missing, request its creation.
			createsubdirs( dirname(dirpath1) )
		}
		# the parent directory exists but not this dir
		dir.create( dirpath1 )
	}
	}

 	use0<-"pairwise.complete.obs" #abbreviation used in correlation options

##MBH TEMPERATURE RECONSTRUCTION ALGORITHMS

  #NHBETA
    ##used in E&E 2005
	#this is to be used with config3
	#series - this is a table listing the number of PC series by 

   NHbeta<-function(proxy,case,series=series.MBH,directory=directory.MBH,method3="test",method4="scale.NH",
	bristlecone.tag=FALSE,NOAMER.CENSORED=FALSE,flag="new") {
	proxy0<-proxy

	nhscale<-array(  nhscale,  dim= c(length(nhscale),1) ) #dim 707 in config3
	np<-length(select) #number of steps
	mm<-period[1]-period[length(period)]

	G<-array(c(rep(NA,16*112)),dim=c(16,112))
	D<- array ( rep(NA,np*112)  , dim=c(112,np)  ) #returned as second item in list
	recon<-rep(NA,581) #reconstruction stepwise
	recontemp<-rep(NA,581) #reconstruction step
	etemp<-rep(NA,581)  #sparse reconstruction step
	e<-rep(NA,581)  #sparse reconstruction stepwise
	ftemp<-rep(NA,581) #other sparse not used presently
	f<-rep(NA,581) #not used
	scale.factor<-rep(NA,11) 
	H<-rep(list(NA),11)  #reconstructed RPCs

      #stepwise calculations
	for (k in 1:11) {
	   #arrange proxy
		proxy<-proxy0
		#base: read in principal component series; 
		#tags for special studies on NOAM.CENSORED and PC4 effect
		for (m in 1:6)  {
			if (series[m,k]>0) {if ((m==3) & (k==11) & NOAMER.CENSORED) load(file.path(pcstore,case,paste(regionid[m],directory[m,k],"CENSORED.tab",sep="."))) else load(file.path(pcstore,case,paste(regionid[m],directory[m,k],"tab",sep=".")))} #this proides for using CENSORED version
			if (series[m,k]>0) { if (m==3 & bristlecone.tag) proxy[,(mannid[m]-1+(1:(series[m,k]-1)))]<-pca[,c(1:3,5:series[m,k])] else proxy[,(mannid[m]-1+(1:series[m,k]))]<-pca[,(1:series[m,k])] }						
			}
		m1<-apply(proxy[503:581,],2,mean,na.rm=TRUE);## means
		s1<-apply(proxy[503:581,],2,sd,na.rm=TRUE);##standard deviations
		for (j in 1:112) {
			temp<-!is.na(proxy[,j]);
			proxy[temp,j]<-(proxy[temp,j]-m1[j])/s1[j]
			}
		selection<-select[[k]]
		roster<-!is.na(proxy[period[k+1],]) #length 112
		
		U<-as.matrix(pc[1:79,selection])
		Ut<-t(U)
  		B<-proxy[503:581,roster]
		G[selection,roster]<- solve(Ut %*% U ) %*% Ut %*% B;  #this is array of regression coefficients in calibration
		P<-diag(weights1[roster]^0.5)  #this is used in weighted regression
		y<-proxy [period[k+1]:(period[1]-1),roster] #stepwise proxy rosters

	#calculate weighting factors	
		C <- P %*%  ( t(G[selection,roster] %*% P) %*% solve ((G[selection,roster] %*% P) %*% t((G[selection,roster] %*%P)))   )

	#calculate reconstructed temperature principal components #this is somewhat redundant
		H[[k]]<-y %*% C 

	#calculate weighting factor for NH temperature construction
		if (length(selection)>1 ) C1<- C %*% diag( lambda[1:16][selection] )    %*% t(eof) [selection,]    else {
			C1<- (C * lambda[1:16][selection])    %*% t(eof) [selection,] }    #dimension 112x1082 #allows for vector form
		C1<-C1[,!nil] #dim 112x1078#
		C2<-C1[,nhtemp.nil] #this picks out NH cells #dim 112x703
		C3<-C2 %*% nhscale #length 112
		D[roster,k]<-C3

	#re-scale RPCs #
		if (method3=="scale.RPC") {
			scale0<-apply(as.matrix(H[[k]] [(1902-period.m[k]+1):(1980-period.m[k]+1),]),2,sd)
			scale0<-scale0/apply(as.matrix(pc[1:79,selection]),2,sd)
			H[[k]]<- scale(H[[k]],center=FALSE,scale=scale0) #this matches sd as per scaling
			if (length(selection)>1 ) V<- H[[k]] %*% diag( lambda[1:16][selection] )    %*% t(eof) [selection,]    else {
				V<- (H[[k]] * lambda[1:16][selection])    %*% t(eof) [selection,] }    #dimension 112x1082 #allows for vector form
			V<-V[,!nil] #dim 112x1078#
			V<-V[,nhtemp.nil] #this picks out NH cells #dim 112x703
			V<-V %*% nhscale #length 201 1
			D[roster,k]<-C3
			recontemp[period[k+1]:(period[1]-1)]<-V  #dim 150x112  * 112x1
			} else recontemp[period[k+1]:(period[1]-1)]<-y %*% C3  #dim 150x112  * 112x1
	
	#re-scale NH reconstruction #this is used in EE to get scale to fit
		if (method4=="scale.NH") {
			scale.factor[k]<-sd(recontemp[(1902-1399):(1980-1399)])
			if (flag=="old") recontemp[period[k+1]:(period[1]-1)]<-recontemp[period[k+1]:(period[1]-1)]*scale.factor[k]/NH.scale.factor else
				recontemp[period[k+1]:(period[1]-1)]<-recontemp[period[k+1]:(period[1]-1)]*NH.scale.factor/scale.factor[k]
			}

	recon[period[k+1]:(period[k]-1)]<-recontemp[period[k+1]:(period[k]-1)]

	#sparse calculations
		E3<-C1[,mask[[2]][!nil]]%*% nhscale[mask[[2]]]
		etemp[period[k+1]:(period[1]-1)]<-y %*% E3  #dim 150x112  * 112x1
		e[period[k+1]:(period[k]-1)]<-etemp[period[k+1]:(period[k]-1)] #this is 219 cell reconstruction
		}
	NHbeta<-list(recon,D,recontemp,scale.factor,e,etemp,f,ftemp,H)
	names(NHbeta)<-c("reconstruction","weights","AD1400","scale","mask.219","mask.219.AD1400","mask.11","mask.11.AD1400","RPC")
	NHbeta
	} 


   ##NHDELTA 
     #this is current (May 2005) version and is reconciled to Ammann and Wahl

     NHdelta<-function(case,proxy.method="MBH98",series=series.MBH,directory=directory.MBH,method3="scale.RPC",method4="default",
	case.NOAMER="MBH98",weight.method="MBH98",sd.rescaling=rep(1,1082),steps=c(1:11),gaspe.method="default") {

	nhscale<-array(  nhscale,  dim= c(length(nhscale),1) ) #dim 707 in config3
	np<-length(select) #number of steps
	mm<-period[1]-period[length(period)]
	D<- array ( rep(NA,np*112)  , dim=c(112,np)  ) #returned as second item in list
	recon<-rep(NA,581) #reconstruction stepwise
	recontemp<-rep(NA,581) #reconstruction step
	sparse.fitted<-rep(NA,581)
	scale.factor<-rep(NA,11) 
	H<-rep(list(NA),11)  #reconstructed RPCs
 
       # stepwise calculations
	for (k in steps[1]:steps[length(steps)]) {
	   #arrange proxy
		proxy<-loadproxy(step0=k,case,proxy.method,mbhdata,series,directory,case.NOAMER,gaspe.method)
		if(proxy.method=="WA") {proxy<-rbind (array ( rep(NA, ncol(proxy)*(581-nrow(proxy))),dim=c( 581-nrow(proxy),ncol(proxy))),proxy);
				proxy<-cbind (proxy,  array ( rep (NA, 581*(112-ncol(proxy))),dim=c(581,(112-ncol(proxy))))) }
		selection<-select[[k]]
		roster<-!is.na(proxy[period[k+1],]) #length 112
		U<-as.matrix(pc[1:79,selection])
		Ut<-t(U)
  		B<-proxy[503:581,roster]
		if(proxy.method=="MBH98") G<-array(c(rep(NA,16*112)),dim=c(16,112)) else G<-array(c(rep(NA,16*ncol(proxy))),dim=c(16,ncol(proxy))) 
		G[selection,roster]<- solve(Ut %*% U ) %*% Ut %*% B;  #this is array of regression coefficients in calibration
		if(weight.method=="MBH98") P<-diag(weights1[roster]^0.5)  else P<-diag (rep(1,sum(roster))) # #this is used in weighted regression
		y<-proxy [period[k+1]:(period[1]-1),roster] #stepwise proxy rosters

	#calculate weighting factors	
		C <- P %*%  ( t(G[selection,roster] %*% P) %*% solve ((G[selection,roster] %*% P) %*% t((G[selection,roster] %*%P)))   )

	#calculate reconstructed temperature principal components #this is somewhat redundant
		H[[k]]<-y %*% C 

	#calculate weighting factor for NH temperature construction
		if (length(selection)>1 ) C1<- C %*% diag( lambda[1:16][selection] )    %*% t(eof) [selection,]    else {
			C1<- (C * lambda[1:16][selection])    %*% t(eof) [selection,] }    #dimension 112x1082 #allows for vector form
		C1<-C1[,!nil] #dim 112x1078#
		##D[roster,k]<-C1 [,nhtemp.nil]%*% (mann.scale[nhtemp.nil]*sd.rescaling[nhtemp.nil])/sum(mu[nhtemp.nil]) #length 112

	#do weighting of RPC
		if (method3=="scale.RPC") {
			scale.est<-apply(as.matrix(H[[k]] [(1902-period.m[k+1]+1):(1980-period.m[k+1]+1),]),2,sd) #adjusted May 2005
			scale.adj<-scale.est/apply(as.matrix(pc[1:79,selection]),2,sd)
			H[[k]]<- scale(H[[k]],center=FALSE,scale=scale.adj) #this matches sd as per scaling
			if (length(selection)>1 ) V<- H[[k]] %*% diag( lambda[1:16][selection] )    %*% t(eof) [selection,]    else {
				V<- (H[[k]] * lambda[1:16][selection])    %*% t(eof) [selection,] }    #dimension 112x1082 #allows for vector form
			V<-V[,!nil] #dim 112x1078#
			#V<-V[,nhtemp.nil] #this picks out NH cells #dim 112x703
			##D[roster,k]<-C3
			recontemp[period[k+1]:(period[1]-1)]<-V[,nhtemp.nil] %*% (mann.scale[nhtemp.nil]*sd.rescaling[nhtemp.nil])/sum(mu[nhtemp.nil]) #edited May 2005
			sparse.fitted[period[k+1]:(period[1]-1)]<- V[,mask.logical] %*% (mann.scale[mask.logical]*sd.rescaling[mask.logical])/sum(mu[mask.logical])		
			} else recontemp[period[k+1]:(period[1]-1)]<-y %*% D[roster,k]  #dim 150x112  * 112x1
	
	#now do reconstruction using weighting factors directly
		if (method4=="scale.NH") {
			scale.factor[k]<-sd(recontemp[(1902-1399):(1980-1399)])
			recontemp[period[k+1]:(period[1]-1)]<-recontemp[period[k+1]:(period[1]-1)]*NH.scale.factor/scale.factor[k]
			}

		recon[period[k+1]:(period[k]-1)]<-recontemp[period[k+1]:(period[k]-1)]
		}
	tags<-list(case,case.NOAMER,gaspe.method,proxy.method,series,directory,method3,method4,weight.method,sd.rescaling,steps)
	NHdelta<-list(ts(recon,start=1400,end=1980),D,ts(recontemp,start=1400,end=1980),scale.factor,H,ts(sparse.fitted,start=1400,end=1980),proxy,tags)
	names(NHdelta)<-c("reconstruction","weights","AD1400","scale","RPC","sparse.AD1400","proxy","tags")
	NHdelta
	}