#1. filtfiltx (Roman Murelka)
   #this is an emulation of Matlab filtfilt (as R filtfilt executes a little differently. See http://www.climateaudit.org/?p=4038#comment-304170

   filtfiltx = function( b, a, x = NULL ) {
	if (is.null(x)) { x = a
	a = b$a
	b = b$b }
	len = length(x)
	nfilt = max(length(b),length(a))
	nfact = 3*(nfilt-1) # padding size
	if (len <= nfact) stop("Data must have length more than 3 times filter order.") # Not enough data
	b = c(b, rep(0,nfilt))[1:nfilt] # zero-pad shorter of b and a
	a = c(a, rep(0,nfilt))[1:nfilt] # ditto
	y = c(2*x[1]-rev(x[2:(nfact+1)]), x, 2*x[len]-rev(x[(len-nfact):(len-1)]))
	y = rev(y[1] + filter(b,a,y-y[1]))
	y = rev(y[1] + filter(b,a,y-y[1]))
	y[(nfact+1):(len+nfact)] } 


#2. mannsmooth
   #this function emulates M08 smoothing lowpassmin. It has been benchmarked against M08 examples.

    mannsmooth=function(x,bwf=bf,M=5,method="filtfiltx") { #filtfiltx
    	mannsmooth = x
 	temp=!is.na(x);N=sum(temp)
 	x=x[temp]; m0=mean(x)
 	if (method=="center") x=x-m0
 	mannextend=array(rep(x,9),dim=c(N,9))
 	dummy1=cbind(rep(mean(x),M),x[M:1],2*x[1]-x[M:1]); dummy1=cbind(dummy1,dummy1,dummy1)
	  #these are 3 padding cases at start: mean; y-reflection; xy-reflection (endpoint is center)
  	dummy2=array(c(rep(mean(x), 3*M),rep(x[(N-1):(N-M)],3),rep(2*x[N]-x[(N-1):(N-M)],3)),dim=c(M,9)); 
	   #three padding cases at end
   	mannextend=rbind(dummy1,mannextend,dummy2);dim(mannextend) #2023 9 
	mannextend=mannextend[,c(1,4,7,2,5,8,3,6,9)]  #to match iconb, icone order
  	working=apply(mannextend,2, function(x) filtfiltx(bwf,x) )  #this uses Roman M's implementation
 	if(method=="filtfilt") working=apply(mannextend,2, function(x) filtfilt(bwf,x) )  
 	test=working[(M+1):(N+M),]
   	residssq=apply(test-x,2,function(x) sum(x^2))  #var in original but doesn't matter
   	index=max((1:9)[residssq==min(residssq)]);index #the last value is chosen
              #tie breaker in Mann is the last series 
  	mannsmooth[temp]=test[,index]
  	if(method=="center") mannsmooth=mannsmooth+m0 
  	mannsmooth
 	}


#3. rescale
   #simple function to rescale one series against another based on period "index" 

	rescale=function (x,target, index=1850:1995) (x- mean(x[index-tsp(x)[1]+1],na.rm=T))*
           sd(target[index-tsp(target)[1]+1])/sd(x[index-tsp(x)[1]+1],na.rm=T)+mean(target[index-tsp(target)[1]+1])


#4. subcount
  #this is a utility function that is used to concisely perform some counting operations   

    subcount=function(x,y,f= function(x) length(x)) {
	  y=factor(y);levels(y)=tapply(x,y,f)
	  subcount=as.numeric(as.character(y))
	  subcount}

    #for example, suppose that you can have multiple proxies in the same gridcell and you want to attach a weight to
    #	


#4.VERIFICATION PERIOD STATISTICS
	#calibration and verificaiton are both vectors of length 2 with start and end period
	#estimator and observed are time series which may have different start and end periods

  verification.stats<-function(estimator,observed=sparse,calibration=c(1850,1949),verification=c(1950,1995))  {
        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)
	mean.obs.cal<-mean ( combine[index.cal,"observed"],na.rm=TRUE  ) #xmean
       	mean.obs.ver<-mean(  combine[index.ver,"observed"],na.rm=TRUE  ) # mx
	cov.ver<-cov ( combine[index.ver,],use=use0) #was test.ver
	cov.cal<-cov(combine[index.cal,],use=use0) #was test
	R2.cal<-(cov.cal[1,2]*cov.cal[2,1])/(cov.cal[1,1]*cov.cal[2,2])
	R2.ver<-(cov.ver[1,2]*cov.ver[2,1])/(cov.ver[1,1]*cov.ver[2,2])
	RE.cal<-1-sum( (combine[index.cal,"estimator"]- combine[index.cal,"observed"]) ^2,na.rm=TRUE)/sum( (mean.obs.cal- 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( (mean.obs.cal- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	CE<- 1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mean.obs.ver- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	verification.stats<-data.frame(RE.cal=RE.cal,RE.ver=RE.ver,R2.cal=R2.cal,R2.ver=R2.ver,CE=CE)#
	verification.stats
	}

#5. MANNIANCPS
   #base case returns NH and SH reconstructions, which are rescaled against instrumental targets in function cpsf
   # verbose: default - results in simplied return; "verbose" - gives intermediates for diagnostics   

   # case: three options- whole, earlym, latem
   # outerlist: two options- outerlist.mbh; outerlist.sensible
   # criterion:	this permits calculations using dendro, noluter, etc.
   # smoothmethod: two options- mann (piecemeal smoothing); sensible (smooth entire series)
   # lat.adjustment: default is quirky Mann displacement in which gridcells are centered 1 deg S of "typical" 2.5 deg centers;
   # M08 input data Basics is required; also passing; cal; ver; period_MBH
 #  k=9;case="latem";criterion=rep(TRUE,1209);outerlist=outerlist.mbh;smoothmethod="mann";lat.adjustment= -1;verbose="default";Info=Basics

   manniancps=function(k,case="whole",criterion=rep(TRUE,1209),outerlist=outerlist.mbh,smoothmethod="mann",lat.adjustment= -1,verbose="default",Info=Basics) {
	proxy=Info$proxy; info=Info$instrumental.info;passing=Info$passing;period_MBH=Info$period_MBH
	calibration=Info$cal[[paste(case)]][1]:Info$cal[[paste(case)]][2]
 	verification=Info$ver[[paste(case)]][1]:Info$ver[[paste(case)]][2]
	annual=Info$criteria$annual
	alias=c("r1850_1995","r1896_1995","r1850_1949");names(alias)=c("whole","earlym","latem")
 
	tempi=!is.na(proxy[period_MBH[k]-tsp(proxy)[1]+1,])&passing[[paste(case)]]&criterion;sum(tempi) ##25 
	    #this is screened so that the proxies are present in the period and meet late, early, whoe criteria 
	outerlist1000=  sapply(outerlist[tempi], function(x) instrumental.info$mannid[match(x,instrumental.info$jones)] )
		#this is subset of applicable proxies
	idsort=sort(unique(c(unlist(outerlist1000))));#idsort  #this is MAnn id (not jones id)

    #make data.frame with relevant gridcell information
	Data=NULL
	for(i in 1:length(outerlist1000)) Data=rbind(Data,data.frame(rep(names(outerlist1000)[[i]],length(outerlist1000[[i]])),outerlist1000[[i]]))
	names(Data)=c("proxy","mannid")
	Data$jones=info$jones[match(Data$mannid,info$mannid)]
	Data$gridcount=subcount(Data$proxy,Data$jones) 
		#this places on each line the count of proxies in that gridcell
	Data$regrid=info$regrid[match(Data$jones,info$jones)]	
	        #this looks up the regrid number for each gridcell  
	Data$lat=info$lat[Data$mannid]#-lat.adjustment
	Data$regrid.lat=info$regrid.lat[Data$mannid]#-lat.adjustment
		#this uses the regrid.lat for later weighting	
	Data$mu= cos( Data$lat*pi/180)	
	Data$flip=1
	regid=sort(unique(Data$regrid));K=length(regid); #K;regid #	15
		#the set of unique regrid cell; and count
 
    #smooth and subset proxy to proxys: two methods "sensible" and "mann"
     if(smoothmethod=="mann"){
 	orientation_lf= !annual & !sign(Info$details$rtable.r1850_1995lf)==sign(Info$details[,paste("rtable.",alias[paste(case)],"lf",sep="")]) ; #sum(orientation_lf&tempb) #16 gross; 6 net
	orientation_hf= annual & !sign(Info$details$rtable.r1850_1995)==sign(Info$details[,paste("rtable.",alias[paste(case)],sep="")]);#c(sum(orientation_hf),sum(orientation_hf&tempb)) 
	proxy[,(orientation_hf|orientation_lf)&tempb]= - proxy[,(orientation_hf|orientation_lf)&tempb]
	Data$flip[!is.na(match(Data$proxy, (1:1209)[(orientation_hf|orientation_lf)&tempb]))]=- Data$flip[!is.na(match(Data$proxy, (1:1209)[(orientation_hf|orientation_lf)&tempb]))]
	  #weird orientation patch in case orientations are flipped

    	tempp=(time(proxy) >=period_MBH[k])&(time(proxy)<=1995)
	working=ts(proxy[tempp,tempi],start=period_MBH[k]);dim(working) #996 25
 	sd.proxy.long=apply(proxy[tempp,tempi],2,sd,na.rm=T);names(sd.proxy.long)=dimnames(proxy)[[2]][tempi];
	temp=!is.na(match(time(working),calibration));#sum(temp)
    	proxys=working;	 
	cutfreq=.05;ipts=10  #ipts set as 10 in Mann lowpass
	bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad
	proxys[,!annual[tempi]]=apply(working[,!annual[tempi]],2,function(x) mannsmooth(x,M=npad,bwf=bf ))
	cutfreq=.1;ipts=10  #ipts set as 10 in Mann lowpass
	bf=butter(ipts,2*cutfreq,"low");npad=1/(2*cutfreq);npad
	proxys[,annual[tempi]]=apply(working[,annual[tempi]],2,function(x) mannsmooth(x,M=5,bwf=bf ) )
           #this finishes piecemeal M08 smooth on truncated data
	}

    if(smoothmethod=="sensible") {		
	tempp=(time(proxy) >=period_MBH[k])&(time(proxy)<=1995)
	proxys=ts(Basics$proxys[temp,tempi],start=period_MBH[k]);#dim(working) #997 25
	sd.proxy.long=apply(proxy[tempp,tempi],2,sd,na.rm=T);names(sd.proxy.long)=dimnames(proxy)[[2]][tempi];
   	   #this is alternative using once-for-all smooth
	   #doesn't bother with Mannian flipping
	}

   #scale the smoothed proxies in short segment calibration period
	temp=!is.na(match(time(proxys),calibration)) ;#sum(temp)   #his is calibration period 
	mean.proxy=apply(proxys[temp,],2,mean,na.rm=T);names(mean.proxy)=dimnames(proxys)[[2]]; 
	sd.proxy=apply(proxys[temp,],2,sd,na.rm=T);names(sd.proxy)=dimnames(proxys)[[2]];
	working = ts( scale(proxys,center=mean.proxy,scale=sd.proxy),start=tsp(proxys)[1])
	Data$sd.proxy.long=sd.proxy.long[paste(Data$proxy)]
	Data$sd.proxy.cal=sd.proxy[paste(Data$proxy)]
           #this makes average of available proxies for utilized gridcells (this is NOT the regridded version)

    #make average of available proxies for utilized gridcells (NOT regridded)
	grids=array(NA,dim=c(nrow(working),length(idsort) ) ) #length idsort
	dimnames(grids)[[2]]=idsort
	average= function(X) if( is.null(ncol(X) )) X else apply(X,1,mean,na.rm=T)
	for (i in 1:length(idsort) ) {
	  tempr= !is.na(unlist(lapply(outerlist1000,function (A) match(idsort[i],A)) ) ) # ;sum(tempr) 
		#this picks out the relevant series in gridcell
	  grids[,i]=average( working[,tempr])
	 }
	grids=ts(grids,start=period_MBH[k])	    
	Data$gridswt=Data$flip*(Data$sd.proxy.long/Data$sd.proxy.cal)*(1/Data$gridcount)

    #rescale gridcell proxy average to gridcell instrumental target	
	mean.obs= info[idsort,paste("mean",case,sep="_")]; #this locates right match in instrumental.info
	sd.obs=info[idsort,paste("sd",case,sep="_")]; names(sd.obs)=idsort
	mean.est=apply(grids[calibration-tsp(grids)[[1]]+1,],2,mean);range(mean.est) #-2.758761e-15  2.145853e-15
	sd.est=apply(grids[calibration-tsp(grids)[[1]]+1,],2,sd);range(sd.est)  # 0.5655191 1.0000000

	regts=array(NA,dim=c(nrow(grids),ncol(grids)) ); dimnames(regts)[[2]]=dimnames(grids)[[2]]
	for( i in 1:ncol(grids)) {
         	regts[,i]=(grids[,i]-mean.est[i])*sd.obs[i]/sd.est[i]+mean.obs[i]
	 }
	regts=ts(regts,start=tsp(grids)[1]);dim(regts)
	Data$regtswt=sd.obs[paste(Data$mannid)]/sd.est[paste(Data$mannid)]

    #regrid N of 30N into 15x15 gridcells
	weighted.average=function(X,w) if( is.null(ncol(X) )) X else apply(X,1,function(x) weighted.mean(x,w,na.rm=T) )
	rescaled=array(NA,dim=c(nrow(regts),length(regid))) #15
	dimnames(rescaled)[[2]]=regid;K=length(regid)
 	for(i in 1:K) {
	  tempr= !is.na(match(info$regrid[idsort],dimnames(rescaled)[[2]][i])) 
	  rescaled[,i]=weighted.average( regts[,tempr],info$mu[idsort[tempr]])
	 }
	rescaled=ts(rescaled,start=tsp(regts)[1],end=1995) #ends one year early

	Data$regrid.count= subcount(Data$mannid,Data$regrid,f=function(x) length(unique(x)))
		#this counts the number of gridcells in each regrid cell and associates count with each line
	Data$regrid.test=subcount(Data$mannid,Data$regrid,f=function(x) sum ( sapply(unique(x), function(y) cos(Data$lat[!is.na(match(Data$mannid,y))][1]*pi/180)  ) ) )
	   #this gives total of gridcell weights
	Data$regridwt=cos(Data$lat*pi/180)/Data$regrid.test
	   #this weights each gridcell in the regrid by cos(latitude)

    #make NH and SH composites (the final output of this function
	###index=sort(unique( Data$regrid));# index  # 38 77 92
	###regrid.weights=cos( sapply(index, function(x) unique(info$regrid.lat[info$regrid==x])  )*pi/180);names(regrid.weights)=index
	####nhtemp=	sapply(index, function(x) unique(info$regrid.lat[info$regrid==x]) >=0 ) ;  nhtemp

	nhtemp= Regrid$lat[regid]>0  
	regrid.weights=cos(Regrid$lat[regid]*pi/180);names(regrid.weights)=regid

	composite= ts( array(NA,dim=c(nrow(rescaled),2)),start=period_MBH[k])
	dimnames(composite)[[2]]=c("NH","SH")
	composite[,"NH"]=apply(rescaled[,nhtemp],1,function(x) weighted.mean(x,regrid.weights[nhtemp],na.rm=T) )
	composite[,"SH"]=apply(rescaled[,!nhtemp],1,function(x) weighted.mean(x,regrid.weights[!nhtemp],na.rm=T) )
	Data$weights=Data$gridswt*Data$regtswt*Data$regridwt*cos(Data$regrid.lat*pi/180)

	if(verbose=="verbose") manniancps=list(idsort=idsort,Data=Data,grids=grids,regts=regts,rescaled=rescaled,recon=composite,working=working)
	if(verbose=="default") manniancps=composite

	manniancps
    } #end manniancps

	
#6. cpsf
   #combines manniancps composite (two hemispheres) and instrumental target name to yield NH, SH and GL CPS reconstructions
   #raw is NH, SH composite	 

    cpsf=function(raw,target="iHAD",index=1850:1995){
   	 cpsf=ts(array(NA,dim=c(nrow(raw),3)),start=tsp(raw)[1]);dimnames(cpsf)[[2]]=c("NH","SH","GL")
	 cpsf[,1] = rescale(raw[,"NH"],instr.smooth[,paste(target,"NH",sep="_")],index)
	 cpsf[,2] = rescale(raw[,"SH"],instr.smooth[,paste(target,"SH",sep="_")],index)
	 cpsf[,3]=apply(cpsf[,1:2],1,mean)
	 cpsf}



##	raw.mbh=manniancps(k=11,criterion=passing$whole, outerlist=outerlist.mbh,lat.adjustment= -1,
#	   smoothmethod="mann",verbose="default") ;tsp(raw.mbh)
##	cps.mbh=cpsf(raw.mbh)