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

  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]))
                                                 

#####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
    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)
  }