read.control<-function(control.file){ #controlfile is of the kind "mannoriginal_original.txt" obj<-list() #Read Reconstruction Identifier skipline<-1 obj$caselabel<-scan(control.file,skip=skipline,nlines=1,what="") #Read Years of Calibration skipline<-skipline+2 yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj$years.calibration<-yy[1]:yy[2] #Read Years of Verification skipline<-skipline+2 yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj$years.verification<-yy[1]:yy[2] #Read Years of Instrumental Data skipline<-skipline+2 yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj$years.instrumental<-yy[1]:yy[2] #Read annual-standardized Instrumental File skipline<-skipline+2 obj$file.instrumental<-scan(control.file,skip=skipline,nlines=1,what="") #Read lat-lon Grid of Instrumental Data skipline<-skipline+2 obj$grid.instrumental<-scan(control.file,skip=skipline,nlines=1,what="") #Read cos.weight for each calibration grid cell skipline<-skipline+2 obj$file.cos.lat.instrumental<-scan(control.file,skip=skipline,nlines=1,what="") #Read cos.weight for each verification grid cell skipline<-skipline+2 obj$file.cos.lat.verification<-scan(control.file,skip=skipline,nlines=1,what="") #Read Mean and StdDev for each Calibration grid cell skipline<-skipline+2 obj$file.means.stdev.instrumental<-scan(control.file,skip=skipline,nlines=1,what="") #Read number of Calibration grid cells skipline<-skipline+2 obj$instrumental.gridcells<-scan(control.file,skip=skipline,nlines=1,what=numeric()) #Boolean for instrumental data in hundreds of a degree (TRUE) or in degrees (FALSE) skipline<-skipline+2 obj$hundredths<-scan(control.file,skip=skipline,nlines=1,what=logical()) #Boolean for scaling of reconstructed PCs (see Scenario 6) skipline<-skipline+2 obj$rpc<-scan(control.file,skip=skipline,nlines=1,what=logical()) #Read File of Verification data skipline<-skipline+2 obj$file.verification<-scan(control.file,skip=skipline,nlines=1,what="") #Read lat-lon information of verification grid skipline<-skipline+2 obj$grid.verification<-scan(control.file,skip=skipline,nlines=1,what="") #Read number of Verification Grid Cells skipline<-skipline+2 obj$verification.gridcells<-scan(control.file,skip=skipline,nlines=1,what=numeric()) #Read number of Proxy Periods to be used in Reconstruction skipline<-skipline+2 obj$proxy.periods<-scan(control.file,skip=skipline,nlines=1,what=numeric()) #Read all sets of Reconstruction Periods with : skipline<-skipline+2 for(i in 1:obj$proxy.periods){ #Read Range of Proxy Data in this period yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj[[paste("years.proxy",i,sep="")]]<-yy[1]:yy[2] skipline<-skipline+2 #Read Number of individual predictor series in this period yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj[[paste("py.proxy",i,sep="")]]<-yy skipline<-skipline+2 #Read List of PCs to be retained from reconstruction yy<-scan(control.file,skip=skipline,nlines=1,what=numeric()) obj[[paste("RetPCs.proxy",i,sep="")]]<-yy skipline<-skipline+2 #Read Label of this period yy<-scan(control.file,skip=skipline,nlines=1,what=character()) obj[[paste("label.proxy",i,sep="")]]<-yy skipline<-skipline+2 #Read name of Proxy file for this period yy<-scan(control.file,skip=skipline,nlines=1,what="") obj[[paste("file.proxy",i,sep="")]]<-yy skipline<-skipline+2 } return(obj)} format.datamatrices<-function(controldata,i){ #controldata is the object produced by read.control #i is the index of the proxy period under exam obj<-list() obj$Nc<-length(controldata$years.calibration) obj$Nv<-length(controldata$years.verification) obj$Nt<-length(controldata$years.instrumental) obj$Mt<-controldata$instrumental.gridcells obj$Mv<-controldata$verification.gridcells obj$cos.lat.instrumental<-scan(controldata$file.cos.lat.instrumental) obj$cos.lat.verification<-scan(controldata$file.cos.lat.verification) grid.cell.mean.stds.vectors<-scan(controldata$file.means.stdev.instrumental) grid.cell.mean.stds.vectors<-t(matrix(grid.cell.mean.stds.vectors,nrow=4, ncol=obj$Mt)) grid.cell.mean<-grid.cell.mean.stds.vectors[,1] grid.cell.stds<-grid.cell.mean.stds.vectors[,2] if(controldata$hundredths){ grid.cell.mean<-grid.cell.mean/100 grid.cell.stds<-grid.cell.stds/100 } obj$grid.cell.mean<-grid.cell.mean obj$grid.cell.stds<-grid.cell.stds #2 Tmat.anomaly<-scan(controldata$file.instrumental) Tmat.anomaly<-t(matrix(Tmat.anomaly, nrow=obj$Nt, ncol =obj$Mt)) Tmat.anomaly<-Tmat.anomaly/obj$cos.lat.instrumental obj$ctrs.Tmat.anomaly<-apply(Tmat.anomaly[,match(controldata$years.calibration, controldata$years.instrumental)],1,mean) obj$stds.Tmat.anomaly<-sqrt(apply(Tmat.anomaly[,match(controldata$years.calibration, controldata$years.instrumental)],1,var)) Tmat.anomaly<-(Tmat.anomaly-obj$ctrs.Tmat.anomaly)/obj$stds.Tmat.anomaly Tmat.anomaly<-Tmat.anomaly*obj$cos.lat.instrumental Tmat.anomaly<-t(Tmat.anomaly) obj$Tmat.anomaly<-Tmat.anomaly Tmat.anomaly.verif<-scan(controldata$file.verification) Tmat.anomaly.verif<-matrix(Tmat.anomaly.verif, nrow=obj$Nv, ncol =obj$Mv) if(controldata$hundredths)Tmat.anomaly.verif<-Tmat.anomaly.verif/100 obj$Tmat.anomaly.verif<-Tmat.anomaly.verif obj$YyearsA<-controldata[[paste("years.proxy",i,sep="")]] file.proxy<-controldata[[paste("file.proxy",i,sep="")]] obj$label.proxy<-controldata[[paste("label.proxy",i,sep="")]] Ny<-length(obj$YyearsA) obj$PyA<-controldata[[paste("py.proxy",i,sep="")]] YmatA<-scan(file.proxy) YmatA<-t(matrix(YmatA,nrow=obj$PyA,ncol=Ny)) ctrsY<-apply(YmatA[match(controldata$years.calibration, obj$YyearsA),],2,mean) stdsY<-sqrt(apply(YmatA[match(controldata$years.calibration, obj$YyearsA),],2,var)) obj$Ymat.anomaly<-t((t(YmatA)-ctrsY)/stdsY) obj$Ymat.calib.anomaly<-obj$Ymat.anomaly[match(controldata$years.calibration, obj$YyearsA),] return(obj)} get.svd<-function(controldata,i,datamatrices,logfile){ #controldata is the object produced by read.control #datamatrices is the object produced by format.datamatrices #logfile can be specified as a /path/filename string #i indexes the proxy period obj<-list() RetPCs<-controldata[[paste("RetPCs.proxy",i,sep="")]] maxPC<-max(RetPCs) write("Retained Number of PCs:",file=logfile,append=TRUE) write(RetPCs,ncol=length(RetPCs),file=logfile,append=TRUE) Tmat.svd<-svd(datamatrices$Tmat.anomaly) obj$TmatU<-Tmat.svd$u #each column represents the PC time series of each EOF obj$TmatV<-Tmat.svd$v #each column represents one EOF (spatial) pattern obj$TmatD<-Tmat.svd$d #each element is the "resolved variance of the corresponding EOF" #Tmat.anomaly= UDV' obj$TmatD.truncated<-diag(obj$TmatD)[RetPCs,RetPCs,drop=FALSE] TmatU.ret.anom.cal<-obj$TmatU[match(controldata$years.calibration, controldata$years.instrumental),RetPCs,drop=FALSE] obj$ctrs.TmatU.ret.anom.cal<-apply(TmatU.ret.anom.cal,2,mean) obj$stds.TmatU.ret.anom.cal<-sqrt(apply(TmatU.ret.anom.cal,2,var)) obj$TmatU.ret.anom.cal<-t((t(TmatU.ret.anom.cal)-obj$ctrs.TmatU.ret.anom.cal)) obj$TmatV.truncated<-obj$TmatV[,RetPCs,drop=FALSE] return(obj)} mann.fit<-function(controldata,datamatrices,svdproduct,logfile){ #controldata is the object produced by read.control #datamatrices is the object produced by format.datamatrices #svdproduct is the object produced by get.svd #logfile can be specified as a /path/filename string obj<-list() TmatU.ret.anom.fitted<-qr.solve(svdproduct$TmatU.ret.anom.cal,datamatrices$Ymat.calib.anomaly) TmatU.ret.anom.fitted.transpose<-t(TmatU.ret.anom.fitted) TmatU.ret.anom.cal.fittedtranspose<-qr.solve(TmatU.ret.anom.fitted.transpose, t(datamatrices$Ymat.calib.anomaly)) TmatU.ret.anom.precal.fittedtranspose<-qr.solve(TmatU.ret.anom.fitted.transpose, t(datamatrices$Ymat.anomaly[-match(controldata$years.calibration, datamatrices$YyearsA),])) TmatU.ret.anom.cal.fitted<-t(TmatU.ret.anom.cal.fittedtranspose) TmatU.ret.anom.precal.fitted<-t(TmatU.ret.anom.precal.fittedtranspose) Error<-svdproduct$TmatU.ret.anom.cal-TmatU.ret.anom.cal.fitted RMSEError.TmatU.fitted<-sqrt(mean((Error)^2)) #added: C. Ammann April 2005: Scale reconstructed PCs if RPC turned on (Default: TRUE) if(controldata$rpc){ stds.TmatU.ret.anom.cal.fitted<-sqrt(apply(TmatU.ret.anom.cal.fitted,2,var)) TmatU.fitted.norm.factor.cal<-svdproduct$stds.TmatU.ret.anom.cal/stds.TmatU.ret.anom.cal.fitted TmatU.ret.anom.cal.fitted<-t(t(TmatU.ret.anom.cal.fitted)*TmatU.fitted.norm.factor.cal) test.TmatU.ret.anom.cal.fitted.normalization<-(sqrt(apply(TmatU.ret.anom.cal.fitted,2,var)))/ svdproduct$stds.TmatU.ret.anom.cal TmatU.ret.anom.precal.fitted<-t(t(TmatU.ret.anom.precal.fitted)*TmatU.fitted.norm.factor.cal) } Tmat.anomaly.cal.fitted<-t(t(TmatU.ret.anom.cal.fitted)+svdproduct$ctrs.TmatU.ret.anom.cal)%*% svdproduct$TmatD.truncated%*%t(svdproduct$TmatV.truncated) Tmat.anomaly.precal.fitted<-t(t(TmatU.ret.anom.precal.fitted)+svdproduct$ctrs.TmatU.ret.anom.cal)%*% svdproduct$TmatD.truncated%*%t(svdproduct$TmatV.truncated) mean.Tmat.anomaly<-mean(datamatrices$Tmat.anomaly[match(controldata$years.calibration, controldata$years.instrumental),]) write(paste("CHECK: Instrumental Calibration Period mean : ",mean.Tmat.anomaly),file=logfile,append=TRUE) write(paste("CHECK: Reconstructed Calibration Period mean : ",mean(Tmat.anomaly.cal.fitted)),file=logfile,append=TRUE) #####Calculation of NH mean values Tmat.anomaly.cal.fitted.cos.stdn.corr<-t(t(Tmat.anomaly.cal.fitted)/ datamatrices$cos.lat.instrumental* datamatrices$grid.cell.stds* datamatrices$stds.Tmat.anomaly) Tmat.anomaly.precal.fitted.cos.stdn.corr<-t(t(Tmat.anomaly.precal.fitted)/ datamatrices$cos.lat.instrumental* datamatrices$grid.cell.stds* datamatrices$stds.Tmat.anomaly) Tmat.anomaly.cal.fitted.NH.mean.corrected<-Tmat.anomaly.cal.fitted.cos.stdn.corr[,1:707] %*% (datamatrices$cos.lat.instrumental[1:707]/ sum(datamatrices$cos.lat.instrumental[1:707])) Tmat.anomaly.precal.fitted.NH.mean.corrected<-Tmat.anomaly.precal.fitted.cos.stdn.corr[,1:707] %*% (datamatrices$cos.lat.instrumental[1:707]/ sum(datamatrices$cos.lat.instrumental[1:707])) #SM added obj$Tmat.anomaly.cal.fitted=Tmat.anomaly.cal.fitted #####Calculation of Beta Resolved Variance value on Grid Scaled_unweighted AND on NHem mean Scaled Tmat.anomaly.instrumental.cos.stdn.corr<-t(t(datamatrices$Tmat.anomaly)/ datamatrices$cos.lat.instrumental* datamatrices$grid.cell.stds* datamatrices$stds.Tmat.anomaly) Tmat.anomaly.cal.cos.stdn.corr<-Tmat.anomaly.instrumental.cos.stdn.corr[match(controldata$years.calibration, controldata$years.instrumental),] Beta.resolved.variance.grid.corrected<-(1-((sum((Tmat.anomaly.cal.fitted.cos.stdn.corr- Tmat.anomaly.cal.cos.stdn.corr)^2))/ (sum((Tmat.anomaly.cal.cos.stdn.corr)^2)))) write(paste("NHem_Grid_level Calibration-RE :",Beta.resolved.variance.grid.corrected),file=logfile,append=TRUE) Tmat.anomaly.cal.NH.mean.corrected<-apply(t(t(Tmat.anomaly.cal.cos.stdn.corr[,1:707])* datamatrices$cos.lat.instrumental[1:707]/ sum(datamatrices$cos.lat.instrumental[1:707])),1,sum) Beta.resolved.variance.NH.mean.corrected<-(1-((sum((Tmat.anomaly.cal.fitted.NH.mean.corrected- Tmat.anomaly.cal.NH.mean.corrected)^2))/ (sum((Tmat.anomaly.cal.NH.mean.corrected)^2)))) write(paste("NHem_Mean Calibration-RE :",Beta.resolved.variance.NH.mean.corrected),file=logfile,append=TRUE) #####Calculation of Beta Resolved Variance value on NH mean Scaled for Verification Period years.noncal.fitted<-datamatrices$YyearsA[-match(controldata$years.calibration,datamatrices$YyearsA)] Tmat.anomaly.verif.NH.mean.corrected<-datamatrices$Tmat.anomaly.verif[,1:172] %*% (datamatrices$cos.lat.verification[1:172]/ sum(datamatrices$cos.lat.verification[1:172])) Tmat.anomaly.verif.fitted.NH.corrected<-Tmat.anomaly.precal.fitted.cos.stdn.corr[match(controldata$years.verification,years.noncal.fitted),c(4:5,16:17,20:23,38:39,43:50,72:73,77:83,86:89,92:93,125:133,138:139,142:143,148:149,152:153,190:200,203:204,211:212,242:245,253:265,298:299,304:309,315:326,328:331,354:357,360:365,371:381,395,406:409,414:415,419,423:430,447:448,465:466,470,474:479,494,516:517,527:528,548,576:577,619:620,632:633,660:661,676:677,702:703)] Tmat.anomaly.verif.fitted.NH.mean.corrected<-Tmat.anomaly.verif.fitted.NH.corrected %*% (datamatrices$cos.lat.verification[1:172]/ sum(datamatrices$cos.lat.verification[1:172])) Verif.overall.mean.Tmat.anomaly.verif.NH.mean.corrected<-mean(Tmat.anomaly.verif.NH.mean.corrected) write(paste("Verification_Mean_Reconstructed : ",mean(Tmat.anomaly.verif.fitted.NH.mean.corrected)), file=logfile,append=TRUE) Verif.fitted.offset <- mean(Tmat.anomaly.verif.fitted.NH.mean.corrected)- Verif.overall.mean.Tmat.anomaly.verif.NH.mean.corrected write(paste("Verification_Mean_Offset from Instrumental : ",Verif.fitted.offset),file=logfile,append=TRUE) Beta.resolved.variance.NH.mean.corrected.verif<-(1-((sum((Tmat.anomaly.verif.fitted.NH.mean.corrected- Tmat.anomaly.verif.NH.mean.corrected)^2))/ (sum((Tmat.anomaly.verif.NH.mean.corrected)^2)))) write(paste("NHem_Mean Verification-RE :",Beta.resolved.variance.NH.mean.corrected.verif),file=logfile,append=TRUE) obj$Tmat.anomaly.cal.fitted.NH.mean.corrected<-Tmat.anomaly.cal.fitted.NH.mean.corrected obj$Tmat.anomaly.precal.fitted.NH.mean.corrected<-Tmat.anomaly.precal.fitted.NH.mean.corrected #SM inserts: obj$beta=c(Beta.resolved.variance.grid.corrected,Beta.resolved.variance.NH.mean.corrected,Beta.resolved.variance.NH.mean.corrected.verif) obj$estimator= ts( c(Tmat.anomaly.verif.fitted.NH.mean.corrected,Tmat.anomaly.cal.fitted.NH.mean.corrected),start=1854) #sparse in ver step; dense in cal step obj$observed= ts( c(Tmat.anomaly.verif.NH.mean.corrected,Tmat.anomaly.cal.NH.mean.corrected),start=1854) #sparse in ver step; dense in cal step obj$verif_stats=verification.stats( obj$estimator,obj$observed) return(obj)} writeout.res<-function(controldata,datamatrices,svdproduct,proxyfit,logfile){ #controldata is the object produced by read.control #datamatrices is the object produced by format.datamatrices #svdproduct is the object produced by get.svd #proxyfit is the object produced by mann.fit #logfile can be specified as a /path/filename string j<-datamatrices$label.proxy write(t(svdproduct$TmatV[,1:5]),ncol=5,file=paste(controldata$caselabel,"/EOFs_calib",j,".txt",sep="")) write(t(svdproduct$TmatU[,1:5]),ncol=5,file=paste(controldata$caselabel,"/PCs_calib",j,".txt",sep="")) write(svdproduct$TmatD,ncol=1,file=paste(controldata$caselabel,"/SVs_calib",j,".txt",sep="")) Tmat.NH.fitted.all<-numeric(length(datamatrices$YyearsA)) Tmat.NH.fitted.all[match(controldata$years.calibration, datamatrices$YyearsA)]<-as.vector(proxyfit$Tmat.anomaly.cal.fitted.NH.mean.corrected) Tmat.NH.fitted.all[-match(controldata$years.calibration, datamatrices$YyearsA)]<-as.vector(proxyfit$Tmat.anomaly.precal.fitted.NH.mean.corrected) write(Tmat.NH.fitted.all, ncol=1, file = paste(controldata$caselabel,"/", paste("NH.Tmat.anomaly.fitted_corrected.",j,".dat",sep=""),sep="")) write("",file=logfile,append=TRUE) write("",file=logfile,append=TRUE) return(NULL) } ##SM: THIS IS INSERTED verification.stats<-function(estimator,observed=sparse,calibration=c(1902,1980),verification=c(1854,1901)) { 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="pairwise.complete.obs") test<-cov(combine[index.cal,],use="pairwise.complete.obs") 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 }