##VARIATION WITH "CORRECT' DATA

library(fields)
library(signal)
library("R.matlab");library("Rcompression")
library(xlsReadWrite)
 	source("d:/climate/scripts/mann.2008/utilities.v6.txt")
  	period_MBH=seq(1800,200,-100);M=length(period_MBH);M
   	use0="pairwise.complete.obs"	

####################
##real_data
############################

#real_mxd
	download.file("http://www.climateaudit.info/data/rutherford05/info.rutherford.tab","temp",mode="wb");load("temp")
	download.file("http://www.climateaudit.info/data/rutherford05/mxd_gridded.tab","temp",mode="wb");load("temp")
	schwein=(1:1209)[details$code==7500]; length(schwein)
	Basics$proxy[,schwein]= NA
	index=match( details$jones[schwein],info.rutherford$jones)
	Basics$proxy[(1400:1994)+1,schwein]= mxd_gridded[,index]	

#actual tornetrask
	download.file("http://www.climateaudit.info/data/briffa/briffa.raw.tab","temp",mode="wb");load("temp")
	torn=briffa.raw[,"TORN"]#
	torn=window(torn,end=max(time(torn)[!is.na(torn)]))
	tsp(torn)
	grep("torn",details$id) #1070
	Basics$proxy[,1070]=NA
	Basics$proxy[(1:1993)+1,1070]=torn

#emulation
	outerlist_method="mann";smoothmethod="mann";lat.adjustment= -1;screenmethod="mann";center_smooth="mann"
	notilj =  is.na(match( 1:1209,1061:1064) )
	K=17
	name0=c("latem","earlym","whole")
		   #used to get NOtilj_nobristle
	criterion=notilj&noluter##;sum(criterion) #1134
	   #used to get NOtilj_nobristle
	emulation=rep(list(NA),K); names(emulation)=period_MBH[1:17]
	for (k in 1:17) {
	emulation[[k]]=rep(list(NA),3); names(emulation[[k]])=name0
	for (j in 1:3) {
	    emulation[[k]][[j]]=manniancps(k,case=name0[j],criterion,verbose="verbose")  #verbose gives interim steps for benchmark
	   }#j
	}
	length(emulation) #17

#cps recons and stats for 4 cases 
	name0=c("latem","earlym","whole")
	target_name=c("CRU","iCRU","HAD","iHAD")
	hemi_name=c("NH","SH")
	name1=c(outer(target_name,hemi_name,function(x,y) paste(x,y,sep="_") ) )
	Target= data.frame(id= name1,target=rep(target_name,2), hemi=as.character(gl(2,4,labels=hemi_name)) )
	cps=rep(list(NA),4); names(cps)=target_name
	for (n in 1:4) {
	  cps[[n]]=rep(list(NA),17)
	  for(k in 1:17) {
	    cps[[n]][[k]]= rep(list(NA),3); names(cps[[n]][[k]])=name0
	      for (j in 1:3) {
	        x=Basics$cal[[ name0[j] ]]
	        cps[[n]][[k]][[ name0[j] ]]= cpsf( emulation[[k]][[ name0[j] ]]$recon,target=Target$target[n], index= x[1]:x[2])
	  }
	} 
	}

	Stat=rep(list(NA),8); names(Stat)=Target$id
	name1=c("century", c(t( outer(c("earlym","latem","avg","adj"),c("RE","CE","r2"), function(x,y) paste(x,y,sep="_") ) )) )[1:12]
	for (j in 1:8) {
	 Stat[[j]]=statf(emulationx=cps[[Target$target[j] ]], target_id=Target$target[j], hemi=Target$hemi[j])
  	}
	
	Real_tornetrask=list(emulation=emulation, cps=cps,stat=Stat)
	save(Real_tornetrask,file="d:/climate/scripts/mann.2008/notes/Real_tornetrask.tab")