##RECONCILE WA
	

###SUMMARY

#1. Load the WA archived proxy and reconstruction from prior WA run-through.
#Then do 1-dimensional matrix algebra and compare result to recon. This yields exact match. 

#2. Repeat with 2-D algebra in 1-D case.  These all reconcile.

#3. Add in PC2 and re-do. Very close match - under 0.1 deg C for nearly all years.

#4. Do AD1450 with one PC first. Close match to AD1400 1-D

###FUNCTIONS
  mannian=function(proxy,weightsprx=rep(1,ncol(proxy)),ipc) {  #this does recon given proxy and weights
	N=nrow(proxy)
	m0=apply(proxy[(N-78):N,],2,mean);sd0=apply(proxy[(N-78):N,],2,sd)
	Y=scale(proxy,center=m0,scale=sd0)
	X=as.matrix(pc[1:79,1:ipc]) # 0.02310907
	sd.obs=apply(X,2,sd)
	P=diag(weightsprx)^2
	rho=t(as.matrix(cor(pc[1:79,1:ipc],Y[(N-78):N,]))); #16 19
	index=(N-78):N
	if(ipc==1){
	   sd.est =  sqrt(diag(  t(rho)%*%P %*% t(Y[index,])%*% Y[index,] %*% P %*% rho) )  /sqrt(78) ;
	   sd.adj=as.matrix( sd.obs[1:ipc]/sd.est[1:ipc])
	   Uhat = Y %*% P %*% rho %*% diag(sd.adj)
	   tau= ts(Uhat %*% weights.nh[1:ipc,],start=tsp(proxy)[1])
	   beta= c(P %*% rho %*% diag(sd.adj) %*% weights.nh[1:ipc,])
	          }
	if(ipc>1) {
	  Cuu= t(U)%*% U; Cuy= t(U)%*%Y[index,]
	  norm.obs= sqrt(diag( t(U)%*% U));	# 0.7293962 0.9303059
	  norm.est = sqrt(diag( Cuu %*% solve( Cuy %*% P%*% t(Cuy)) %*% Cuy%*% P%*% t(Y[index,]) %*% Y[index,] %*%P%*%  t(Cuy) %*% solve( Cuy %*% P%*% t(Cuy)) %*% Cuu ));
	  Uhat= Y%*%P%*% t(Cuy) %*% solve( Cuy %*% P%*% t(Cuy)) %*% Cuu %*% diag(norm.obs/norm.est)
	  tau=  ts( Uhat%*%as.matrix(weights.nh[1:ipc]), start=tsp(proxy)[1])
	  beta=c( t(Cuy) %*% solve( Cuy %*% P%*% t(Cuy)) %*% Cuu %*% diag(norm.obs/norm.est) %*% as.matrix(weights.nh[1:ipc]))
	 	}	
	mannian=list(tau,beta,Uhat)
	names(mannian)=c("series","weights","RPC")
	mannian}


###
#CREATE WAHL OBJECT
######################
	source("http://www.climateaudit.org/scripts/wahl/mann_masterfile.windows.txt") 
	    #takes a little while: will say Read nn items 	
  	    #creates large R-list object wahlfit.tab 
		#this is from prior run of WA results and is large R-list

     #confirm that you've got the right object
	W=wahlfit[[1]];
	W[[2]]$TmatD[1:4] 
		#1] 154.02858  75.17558  63.15213  61.86446
	attributes(W[[2]])
	  #$names
	  #[1] "TmatU"                   "TmatV"                   "TmatD"                   "TmatD.truncated"        
	  #[5] "ctrs.TmatU.ret.anom.cal" "stds.TmatU.ret.anom.cal" "TmatU.ret.anom.cal"      "TmatV.truncated"        
	  # These are WA objects
	
 	nhtemp=c(rep(TRUE,707),rep(FALSE,1082-707))
        weights.nh=diag(W[[2]]$TmatD[1:16]) %*% t(W[[2]]$TmatV[nhtemp,1:16])%*%diag(W[[1]]$stds.Tmat.anomaly[nhtemp])%*% as.matrix(W[[1]]$grid.cell.stds[nhtemp])/sum(W[[1]]$cos.lat.instrumental[nhtemp])
	c(weights.nh)[1:5]
	  # [1]  1.81831789 -0.48908006  0.08932557  0.53319471  0.08284691

######################
##1.  AD1400 1-D ALGEBRA
#################
	W=wahlfit[[1]];#names(W)
	weights.nh=diag(W[[2]]$TmatD[1:16]) %*% t(W[[2]]$TmatV[nhtemp,1:16])%*% diag(W[[1]]$stds.Tmat.anomaly[nhtemp])%*%as.matrix(W[[1]]$grid.cell.stds[nhtemp])/sum(W[[1]]$cos.lat.instrumental[nhtemp])
	proxy=ts( wahlfit$"1400"[[1]]$Ymat.anomaly,start=1400);dim(proxy)
	   # 581  22
	recon=ts(c(wahlfit$"1400"[[3]]$Tmat.anomaly.precal.fitted.NH.mean.corrected,wahlfit$"1400"[[3]]$Tmat.anomaly.cal.fitted.NH.mean.corrected),start=1400);length(recon)
	   #583
	   #this is WA version
	pc=wahlfit$"1400"[[2]]$TmatU[1:79,1:16] ;pc[1:5,1] #WA pc 
	   #-0.04026887 -0.09598517 -0.17558648 -0.08449567 -0.03783115

	ipc=1
	N=nrow(proxy)
	m0=apply(proxy[(N-78):N,],2,mean);sd0=apply(proxy[(N-78):N,],2,sd)
	Y=scale(proxy,center=m0,scale=sd0)
	X=as.matrix(pc[1:79,1:ipc]) # 0.02310907
	sd.obs=apply(X,2,sd)
	weightsprx=rep(1,ncol(Y))
	P=diag(weightsprx)^2
	rho=t(as.matrix(cor(pc[1:79,1:ipc],Y[(N-78):N,]))); #16 19
	index=(N-78):N
	   sd.est =  sqrt(diag(  t(rho)%*%P %*% t(Y[index,])%*% Y[index,] %*% P %*% rho) )  /sqrt(78) ;
	   sd.adj=as.matrix( sd.obs[1:ipc]/sd.est[1:ipc])
	   Uhat = Y %*% P %*% rho %*% diag(sd.adj)
	   tau= ts(Uhat %*% weights.nh[1:ipc,],start=tsp(proxy)[1])
	   beta= P %*% rho %*% diag(sd.adj) %*% weights.nh[1:ipc,]
     	max(tau-recon) # 1.165734e-15
	   #this shows match to WA result

	test=mannian(proxy,weightsprx=rep(1,ncol(proxy)),ipc=1)
	max(test$series-recon) #  1.165734e-15
		#two results match

	par(mar=c(3,4,3,1))
	plot(1400:1980,test$series,type="l",main="WA Reconstruction - AD1400 Step",xlab="",ylab="deg C")
	abline(h=0)