readmann =function(k, url= "d:/climate/data/mann.2008/proxy/A_allproxyoriginal.zip",id=idorig) {
  handle <- unz(url,paste(idorig$id[k],"ppd",sep=".") ,"r") #formerly V4.M02
  Data <- scan(handle);length(Data)
  Data=t(array(Data,dim=c(3,length(Data)/3)))
  Data=data.frame(Data[4:nrow(Data),])
  names(Data)=c("year","proxy","count")
  close(handle)
  readmann=Data
readmann}

#strip
	#FUNCTION TO LISTFILES AND REMOVE SUFFIXES
   strip<-function(h,n) {
	h<-unlist(strsplit(h, "\\."))
	h<-array(h,dim=c(n,length(h)/n))
	h<-t(h)
	strip<-h[,1]
	strip
	}


filter.combine.pad=function(x,a,M=NA,method="norm") {
	temp<-!is.na(x)  
   	w<-x[temp]
	N<-length(w)
	if(is.na(M)) M<-trunc (length(a)/2)+1
	if(method=="norm")  w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M)  )
	#if(method=="reflect") w<-c (  w[M:1],w,w[M:(N-M+1):N])
	if(method=="mannomatic") {w<-c( c (w[1]-(w[(M+1):2]-w[1])),w, c(w[N]- (w[(N-1):(N-M)]-w[N]) ) ) }
	  if(method=="mann2")  w<-c (rep(w[1],M),w,rep(w[N],M)  )
	   y<-filter(w,a,method="convolution",sides=2)
	y<-y[(M+1):(length(y)-M)]
	z<- rep(NA,length(x))
	z[temp]<-y
	if(class(x)=="ts") z<-ts(z,start=tsp(x)[1],freq=tsp(x)[3])
	   filter.combine.pad<-cbind(x,z)
	   #filter.combine<-y[,2]
	filter.combine.pad
	   }

##gaussian.filter.weights
   gaussian.filter.weights<-function(N,year) {
		#N number of points; year - is number of years
	sigma<-year/6
	i<-( (-(N-1)/2):((N-1)/2)   )  /sigma
	w<- ((2*pi)^-0.5) * exp(-0.5* i^2)
	gaussian.filter.weights<-w/sum(w)
	gaussian.filter.weights
	}

	#GAUSSIAN FILTER
	#see http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf
	#The appropriate Gaussian distribution therefore has standard deviation
	# s = ? / 6,
	# where ? is the desired wavelength at which the amplitude of frequency response is 0.5. The
	# filter weights are obtained by sampling the pdf of the standard normal distribution at
	# values  i/s  where i={-(N-1)/s, ..., 0, .. (N-1)/s }
	# MEKo has values to verify for case of N=21, year=10 #this is verified

##truncated.gauss.weights
   truncated.gauss.weights<-function(year) {
	a<-gaussian.filter.weights(2*year+1,year)
	temp<-(a>0.05*max(a))
	a<-a[temp]
	a<-a/sum(a)
	truncated.gauss.weights<-a
	truncated.gauss.weights
	}


##JONES
		#JONES LOOKUP
   jones<-function(lat,long) {
	i<-18-floor(lat/5);
	j<- 36+ceiling(long/5);
	jones<-72*(i-1) + j
	jones
	}

##jones.inv
   jones.inv<-function(k) {
	j<- k%%72
	temp<-k-j
	i<-(temp/72)+1
	long<-(j-36)*5-2.5
	lat<-(18-i)*5+2.5
	jones.inv<-c(lat,long)
	jones.inv
	}


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


  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
 	}




#this extracts an underlying decadal series from a linear smoothed time series

     extractf=function(x) {
 	r1=range(time(x))
	year=c(time(x))
 	y=diff(x)
 	z=round(diff(y),2)
 	temp=!(z==0)&!is.na(z)
 	time0=sort(c(r1,time(z)[temp]-1))
 	rx=round(x,2)
	temp=!is.na(match(year, time0));sum(temp)
	test=data.frame(time(x)[temp],rx[temp]/units)
	names(test)=c("year","proxy")
		#sum(!is.na(test$proxy[test$year>=1850])) #9
	extractf=test# list(data=test,x=x,y=y,z=z,r1=r1)
	extractf
	}

	##k=202;units=1;x=units*g(mannorig[[details$old[k] ]])


  #collects according to X$year cumulative to end	
    irrmatch=function(A,y) {
	gridded=data.frame(c(time(y)),y);names(gridded)=c("year","target")
  	h= function( x) if(x[1]<=max(A$year)) min(A$year[ A$year>=x[1] ] ) else NA
	gridded$d= apply(gridded,1,h)
	working=tapply(gridded$target,gridded$d,mean,na.rm=T)
	test=match(A$year,as.numeric(names(working)));temp=!is.na(test)
	A$grid=NA;A$grid[temp]=working[test[temp]]
	irrmatch=A
	irrmatch}

  


  #rescales
	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. FUNCTION TO RETURN 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
	}


mannstats= 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)
	combine=data.frame(combine)
	mean.obs.cal<-mean ( combine$observed[index.cal],na.rm=TRUE  ) #xmean
       	mean.obs.ver<-mean(  combine$observed[index.ver],na.rm=TRUE  ) # mx
	resid.ver= sum((combine$estimator[index.ver]-combine$observed[index.ver])^2)
	ssx.ver = sum( combine$observed[index.ver]^2) ############## mean is NOT subtracted for some reason
	betamulti=1-resid.ver/ssx.ver
	rcal=cor(combine$observed[index.cal],combine$estimator[index.cal],use="pairwise.complete.obs")
	rver=cor(combine$observed[index.ver],combine$estimator[index.ver],use="pairwise.complete.obs")
	R2.cal=rcal^2; R2.ver=rver^2
	RE.ver=1-resid.ver/sum( (combine$estimator[index.ver]-mean.obs.cal)^2)
	CE= 1-resid.ver/sum( (combine$estimator[index.ver]-mean.obs.ver)^2)
	mannstats<-data.frame(RE.ver=RE.ver,R2.cal=R2.cal,R2.ver=R2.ver,CE=CE)#
	mannstats
	}
	
	#mannstats(x,target)


  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}

   manniancps=function(k,criterion=passing$whole,outerlist=outerlist.mbh,smoothmethod="mann",lat.adjustment= -1,info=instrumental.info,verbose="default") {
	# outerlist=outerlist.mbh;smoothmethod="mann";outermethod="mann";info=instrumental.info	

	tempi=!is.na(proxy[period_MBH[k]-tsp(proxy)[1]+1,])&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) info$mannid[match(x,info$jones)] )
		#this is subset of applicable proxies
	idsort=sort(unique(c(unlist(outerlist1000))));idsort #21  #this is MAnn id (not jones id)

     #make sub-data.frame showing of 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)
	Data$regrid=info$regrid[match(Data$jones,info$jones)]	
	Data$regrid.count= subcount(Data$mannid,Data$regrid,f=function(x) length(unique(x)))
	Data$sd=info$sd[Data$mannid]
	Data$lat=info$lat[Data$mannid]#-lat.adjustment	
	Data$mu= cos( Data$lat*pi/180)	
	Data$weight= Data$sd* Data$gridcount^-1  *Data$mu * Data$regrid.count^-1
	regid=sort(unique(Data$regrid));K=length(regid);K;regid
 		
    #make working proxy subset: two methods "sensible" and "mann"
     temp=(time(proxy) >=period_MBH[k])&(time(proxy)<=1995)
     if(smoothmethod=="mann"){
	working=ts(proxy[temp,tempi],start=period_MBH[k]);#dim(working) #997 25
	temp=(time(working) >=1850)&(time(working)<=1995)
       #smooth working subset
 	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 ) )

	temp=(time(proxys)>=1850)&(time(proxys)<=1995) ;sum(temp)
	   #his is period in gridboxcps
	mean.proxy=apply(proxys[temp,],2,mean,na.rm=T);names(mean.proxy)=dimnames(proxys)[[2]]; mean.proxy["385"] # 7.585251 
	sd.proxy=apply(proxys[temp,],2,sd,na.rm=T);names(sd.proxy)=dimnames(proxys)[[2]]; sd.proxy["385"]
	working = ts( scale(proxys,center=mean.proxy,scale=sd.proxy),start=tsp(proxys)[1])
 	##apply(working[temp,],2,mean) #vary
	      #this makes average of available proxies for utilized gridcells (NOT regridded)
	}
    if(smoothmethod=="sensible") {		
	temp=(time(proxy) >=period_MBH[k])&(time(proxy)<=1995)
  	working=ts(proxyscaled[temp,tempi],start=period_MBH[k]);#dim(working) #997 25
	temp=(time(working) >=1850)&(time(working)<=1995)
	}

   #make average of scaled proxies in applicable gridcells
	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])
	temp=(time(grids) >=1850)&(time(grids)<=1995)

  #rescale gridcell proxy average to gridcell instrumental target	
	mean.obs= info$mean[idsort];sd.obs=info$sd[idsort]; c(mean.obs[9],sd.obs[9]) # 0.03850556 0.20116227
	mean.est=apply(grids[(1850:1995)-tsp(grids)[[1]]+1,],2,mean);mean.est #all zero
	sd.est=apply(grids[(1850:1995)-tsp(grids)[[1]]+1,],2,sd);sd.est  #some are 1; some are less than 1
	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)
	temp=(time(grids) >=1850)&(time(grids)<=1995)

  #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
   
  #make NH and SH composites
	index=sort(unique( Data$regrid));index  # 38 77 92
		#regrid[index,]
	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
  
	manniancps= ts( array(NA,dim=c(nrow(rescaled),2)),start=period_MBH[k])
	dimnames(manniancps)[[2]]=c("NH","SH")
	manniancps[,"NH"]=apply(rescaled[,nhtemp],1,function(x) weighted.mean(x,regrid.weights[nhtemp],na.rm=T) )
	manniancps[,"SH"]=apply(rescaled[,!nhtemp],1,function(x) weighted.mean(x,regrid.weights[!nhtemp],na.rm=T) )

   ##RECON FOR ALL VERSIONS
	#estimate$recon=ts(array(NA,dim=c(nrow(estimate$raw),ncol(instr.smooth)) ),start=tsp(estimate$raw)[1]) #1000 1995
	#for(i in 1:8) estimate$recon[,i]=rescale(x=estimate$raw[,1+floor((i-1)/4)],target=instr.smooth[,i],index=1850:1995)
	#manniancps=estimate
	if(verbose=="verbose") manniancps=list(idsort=idsort,Data=Data,grids=grids,regts=regts,rescaled=rescaled,recon=manniancps,working=working)
	manniancps
   } #end -k loop


 #takes one list and derive cps re-scaled versions 
    cpsf=function(raw,calibration) {
	 cpsf=ts(array( NA,dim=c(nrow(raw),ncol(instr.smooth)) ),start=tsp(raw)[1]) #1000 1995
	 dimnames(cpsf)[[2]]= c("NH_CRU","NH_iCRU", "NH_HAD",  "NH_iHAD", "SH_CRU", "SH_iCRU", "SH_HAD",  "SH_iHAD")
	 for(j in 1:8) cpsf[,j]=rescale(x=raw[,1+floor((j-1)/4)],target=instr.smooth[,j],index=calibration) #earlym
	cpsf}