##LONG FUNCTION TO PRODUCE RPC FROM PROXIES USING MANN ROSTERS
RPC<-function(proxy) {
mm<-581

##RUNCATE AND NORMALIZE EACH PROXY SERIES BASIS 1902-1980
proxy<-ts(proxy[1:mm,],start=1400,end=(1399+mm))
m1<-apply(proxy[503:581,],2,mean,na.rm=TRUE);## means
s1<-apply(proxy[503:581,],2,sd,na.rm=TRUE);##standard deviations
for (k in 1:112) {
temp<-!is.na(proxy[,k]);
proxy[temp,k]<-(proxy[temp,k]-m1[k])/s1[k]
}

period.m<-c(1981,1970,1820,1800,1780,1760,1750,1730,1700,1600,1500,1450,1400)
period<-period.m-1399
select<-list(c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),c(1:5,7,9,11,14:16),
c(1:5,7,9,11,15),c(1:3,5,7,9,11,15),c(1,2,5,11,15),c(1,2,5,11,15),
c(1,2,11,15),c(1,2),c(1,2),c(1))
np<-length(select)

H<-array(c(rep(NA,mm*16)),dim=c(mm,16))
G<-array(c(rep(NA,16*112)),dim=c(16,112))
C<- array ( rep(NA,np*16*112)  , dim=c(np,16,112)  ) #not used in this version

#NOW 1970-1980
k<-1
selection<-select[[k]]
for (j in 571:mm) {
temp<-!is.na(proxy[j,]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X4+X5+X7+X9+X11+X14+X15+X16-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}


#FIRST 1820-1970
k<-2
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X4+X5+X7+X9+X11+X14+X15+X16-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}


#NEXT 1800-1820. 102 proxies, 11 EVs
k<-3
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X4+X5+X7+X9+X11+X14+X15+X16-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}


##NEXT 1780-1800. 97 proxies,11 EVs
k<-4
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X4+X5+X7+X9+X11+X14+X15+X16-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}


#NOW CHANGE EV BASIS
#FIRST TRANCH 1760-1779 
k<-5
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X4+X5+X7+X9+X11+X15-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#NOW FOR STRANGE 1750-1759. USE 7 and 9 instead of 6 and 8
k<-6
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X3+X5+X7+X9+X11+X15-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#NOW CHANGE EV BASIS AND PROXIES FOR 1730-1749
k<-7
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X5+X11+X15-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#JUST CHANGE PROXIES FOR 1700-1729
k<-8
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X5+X11+X15-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#CHANGE SELECTION AND PROXIES FOR 1600-1699
k<-9
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2+X11+X15-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#CHANGE SELECTION AND PROXIES FOR 1500-1599
k<-10
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#
#JUST CHANGE PROXIES FOR 1450-1499
k<-11
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1+X2-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

#CHANGE SELECTION AND PROXIES FOR 1400-1449
k<-12
selection<-select[[k]]
temp<-!is.na(proxy[period[k+1],]) #length 112
G<-Gcalc(selection,temp,proxy)
Gt<-t(G)
z<-data.frame(Gt)
Y<-c(1:ncol(G))
z<-cbind(Y,z)
for (j in period[k+1]:(period[k]-1)) {
z[,1]<-proxy[j,]
fm<-lm(Y~X1-1,data=z[temp,],weights=weights1[temp])
H[j,selection]<-fm$coefficients
}

RPC<-list(H,C)
return(RPC)
}  #end function


Gcalc<-function(selection,roster,proxy)  {
G<-array (c (rep(NA,16*112)), dim=c(16,112))
U<-pc[,selection]
Ut<-t(U)
G[selection,roster]<- solve(Ut %*% U ) %*% Ut %*% proxy[503:581,roster];  
Gcalc<-G
Gcalc
}

#this redoes with all 16 EVs using rpc.function.3A

nh.alt<-function(rpc,eof,lambda,select)  {
#RECONSTRUCT w-estimate
#the MBH periods are 11 EV 1780-1980; 9 from 1760-1779; 8 from 1750-1759; 5 from 1700 to 1749; 4 from 1600 to 1699; 2 from 1450 to 1599; 1 from 1400 to 1449.
#rpc<-new[1:581,]
#count<-apply(!is.na(rpc),1,sum)
#cbind(1400:1980,count[1:581])
w<-array  (c(rep(NA,581*1082))  ,  dim=c(581,1082)  )
period.m<-c(1981,1970,1820,1800,1780,1760,1750,1730,1700,1600,1500,1450,1400)
period<-period.m-1399
np<-length(select)
for (iter in 1:(np-1))  {
w[period[iter+1]:(period[iter]-1),]<- rpc[period[iter+1]:(period[iter]-1),select[[iter]]]  %*% diag( lambda[1:16][select[[iter]] ] )    %*% t(eof) [select[[iter]],]     
}
for (iter in np:np)  {
if (length(select[[np]])==1) w[period[iter+1]:(period[iter]-1),]<- as.matrix (rpc[period[iter+1]:(period[iter]-1),select[[iter]]])  %*%   as.matrix(lambda[1])     %*%    t( as.matrix(  t(eof) [select[[iter]],])) else w[period[iter+1]:(period[iter]-1),]<- rpc[period[iter+1]:(period[iter]-1),select[[iter]]]  %*% diag( lambda[1:16][select[[iter]] ] )    %*% t(eof) [select[[iter]],]     
}
w<-ts(w,start=1400,end=1980)
dim(w) #[1]  581 1082

#EXPANSION AND NH AVERAGE
load("c:/climate/data/mann/gridpoints.tab")
nhtemp<-(grid2[,4]>0)
load("c:/climate/data/mann/nil.tab")
load("c:/climate/data/mann/mann.scale.tab")
mann.scale<-mann.scale[!nil]/100
w<-w[,!nil]
nhtemp<-nhtemp[!nil]
w<-scale(w,center=FALSE,scale=1/mann.scale)
nh.alt<-apply(w[,nhtemp],1,mean)
nh.alt
}