#SOME UTILITIES 
##JONES GRIDCELL LOOKUP
	plotf=function(x,row=1,col=1,K=4,L=2,ylab="",scale=TRUE,xlim0=c(300,2000),tag=NULL) {
	  if(scale) ylim0=c(-2,2) else ylim0=range(x,na.rm=T)
 	  par0= c(3,0,0,0)
	  if(col==1) par0[2]=3 else par0[2]=0
	  if(row==K) par0[1]=3 else par0[1]= 0
	  if(row==1) par0[3]=2 else par0[3]=0
	  if(col==L) par0[4]=1 else par0[4]=0
	  par(mar=par0)
	  plot(x,type="l",col="grey80",xlim=xlim0,ylim=ylim0,axes=FALSE,ylab="",xlab="")
	  lines(f(x),col=2,lwd=2)
	  grid()
	  if(!is.null(tag))  text(xlim0[1],.9*ylim0[2],pos=4,font=2,tag)
	 abline(h=0,lty=3)
	  if(row==K) axis(side=1) else axis(side=1,labels=FALSE)
	  if (col==1) axis(side=2) else axis(side=2,labels=FALSE)
	  text(300,1.7,ylab,pos=4,font=2) 
	 box()
 } 	  


 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
	}

  ts.annavg=function(X,M=6) { 
        if(class(X)=="ts") {start0<-floor(tsp(X)[1]);X<-ts.array(X)} else 
        if( ((class(X)=="matrix")|(class(X)=="data.frame") ) &(ncol(X)==12)) start0<-as.numeric(row.names(X)[1])  else {X
                start0<-X[1,1]; X<-X[,2:13]}
        count<-apply(!is.na(X),1,sum);temp<-(count<M)
        ts.annavg<-apply(X,1,mean,na.rm=TRUE)
        ts.annavg[temp]<-NA
        ts.annavg<-ts(as.matrix(ts.annavg),start=start0)
        ts.annavg
	}
   
	
     fill.array<-function(X,method="column") {
	if (method=="row.name") {xyear<-as.numeric(row.names(X));X<-X[,c(13,1:12)] }
	if(class(X)=="matrix") {xyear<-X[,1];X<-X[,2:13]} 
	if (class(X)=="data.frame") {xyear<-X$year;X<-as.matrix(X[,2:13])}
	years<-min(xyear,na.rm=TRUE):max(xyear,na.rm=TRUE)
	fill.array<-array(NA,dim=c(length(years),12))
	temp<-!is.na( match(years,xyear) )
	fill.array[temp,]<-X
	row.names(fill.array)<-years
	dimnames(fill.array)[[2]]<-month0
	fill.array
	}

	ts.array<-function(x) {
	year1<- floor(tsp(x)[1]):floor(tsp(x)[2])
	month1<- round(12*(tsp(x)[1:2]%%1),0)
	x<-c(rep(NA,month1[1]),x ,rep(NA, 11-month1[2]))
	ts.array<-t(array(x,dim=c(12,length(x)/12)) )
	row.names(ts.array)<-year1
	ts.array
	}


#NORMAIZE
	anom_old=function(x,reference=1961:1990,M=15) {
		tsp0=floor(tsp(x))
 	   	month= 1+round(12*  c(time(x))%%1)
		temp=!is.na(match(reference,tsp0[1]:tsp0[2]))
		if(sum(temp)<M) anom=NA else {
		  K=floor(tsp(x)[1])
		  if( is.null(dim(x))) {
		   x=c( rep(NA,month[1]-1),x ,rep(NA,12- month[length(month)]))
		   test=t(array(x,dim=c(12,length(x)/12)) )
		   m0=apply(test[(reference[temp])-K+1,],2,mean,na.rm=TRUE)
		   test=scale(test,center=m0,scale=FALSE)
		   anom=round(ts(c(t(test)),start=c(K,1),freq=12),2)
		   anom} else {anom=x;
		   for(i in 1:ncol(x)) {test=t(array(x[,i],dim=c(12,nrow(x)/12)) )
		   m0=apply(test[(reference[temp])-K+1,],2,mean,na.rm=TRUE)
		   test=round(scale(test,center=m0,scale=FALSE),2)
		   anom[,i]=c(t(test))}}}
		   anom
		}
			

trim=function(x) window( x, start= min(time(x)[!is.na(x)]), end= max(time(x)[!is.na(x)]) )

anom=function(x,start=-999999,end=999999,Round=2) { #assumes time series
	  work= ts( rep(1:12,ceiling(tsp(x)[2]-tsp(x)[1]+1)), start= floor(tsp(x)[1]),freq=12) 
	  work=window(work,start=tsp(x)[1],end= tsp(x)[2] )
	  anom=month=factor(work)
	  temp= time(x)>=start & time(x)<= (end+1)
	  levels(anom)=unlist(tapply(x[temp],month[temp],mean,na.rm=TRUE))
	  anom=  round(x- as.numeric(as.character(anom)),Round) 
	return(anom)
}	


##annavg
	#CALCULATE ANNUAL AVERAGE FROM MONTHLY TIME SERIES
	#even number of months
   annavg<-function(x) {
	n<-length(x)
	m<-n-12*floor(n/12)
	if(m>0) x<-c(x, rep(NA,12-m))
	years<-length(x)/12
	x<-array(x,dim=c(12,years))
	annavg<-apply(x,2,mean,na.rm=TRUE)
	return(annavg)
	}

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

##split
   split<-function(h,char) {
	y<-unlist(strsplit(h,paste("\\",char,sep="")))
	y<-array(y,dim=c(2,length(y)/2))
	split<-t(y)
	split<-split[,1]
	split
	}


##smooth and pad: requires no gaps 
filter.combine.pad=function(x,a=truncated.gauss.weights(21),M=NA,method="norm") {
   temp<-!is.na(x) 
   w<-x[temp]
   N<-length(w)
   if(is.na(M)) M<-trunc (length(a)/2)
   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)  )
   if(method=="linear"){
	fm1=lm (w[1:M]~I(1:M)); 
	fm2= lm (w[N+1-(1:M)]~I(N+1-(1:M))); 
   	w<-c (fm1$coef[1]+fm1$coef[2]*( -(M-2):0) ,w,fm2$coef[1]+fm2$coef[2]*( N+M-((M-1):1)) )
     }
    	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
	}

 maketable=function(X,align=-999,widths=-999,out= "d:/temp/temp.txt"){
	N=nrow(X);M=ncol(X); 
	if (align[1]== -999) right=rep(TRUE,M) else right=align
	code=rep("",M);code[right]=" <p align=&quotright&quot>"
	if(widths[1]== -999) pct= floor(100/M) else pct=widths

	write.table("<table border=&quot1&quot style=&quotborder-collapse: collapse&quot bordercolor=&quot#111111&quot width=&quot80%&quot id=&quotAutoNumber1&quot>",file=out,col.names=FALSE,row.names=FALSE,quote=FALSE)
   	for (j in 1:N){
	 write.table("<tr>",append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE)
	  for(i in 1:M){
	    write.table( paste("<td width=&quot",pct[i],"%&quot>",code[i],X[j,i],"</td>",sep=""),append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE)
   	    }
	 write.table("</tr>",append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE)
	 }
	write.table("</table>", append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE)
     }
	
#TO EXTEND SERIES BY PERSISTENCE
extend.persist<-function(tree) {
	extend.persist<-tree
	for (j in 1:ncol(tree) )  {
		test<-is.na(tree[,j])
		end1<-max ( c(1:nrow(tree)) [!test])
		test2<-( c(1:nrow(tree))>end1) & test
		extend.persist[test2,j]<-tree[end1,j]
		}
	extend.persist
	}	

make.ts=function(A,h=identity,interp =TRUE) {
  names(A)=c("year","proxy")
  A=A[!is.na(A$year),]	
  if(! interp) {
    time0=min(A$year):max(A$year)
    test=match(time0, A$year)
    temp=!is.na(test) 
    x=rep(NA,length(time0))
    x[temp]= A$proxy[test[temp]]
    return(x) 
    }
  if(interp) {	
    h=approxfun(A$year,A$proxy) 
    start0=ceiling(min(A$year));end0=floor(max(A$year))
    x=ts( h(start0:end0),start= start0)
    return(x)
   }
}


make.tsm=function(A) {
  	names(A)=c("year","month","proxy")
  	A=A[!is.na(A$year),]	
	target= sort(outer( min(A$year):max(A$year),1:12, function(x,y) 100*x+y) )
	A$time=100*A$year+A$month
	target=data.frame(time=target)
	target$x=NA
	target$x=A$proxy[match(target$time,A$time)]
	y=ts(target$x,start=min(A$year),freq=12)
	return(y)
}
	



#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
  use0="pairwise.complete.obs"

  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=use0)
	
	test<-cov(combine[index.cal,],use=use0)
	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
	}



trend =function(x,method="AR1" ) { 

  if (method=="short") return(round(lm(x~ I((1:length(x))/120))$coef[2],3))

   if (method=="ols") {
	fm= lm(x~ I((1:length(x))/120))
  	work= round( c(unlist(summary(fm)$coef[2,1:3]),fm$df)	[c(1,4,2,3)],4)
	names(work)=c("trend","df","se","t")
	return(work)
	}
   if (method=="AR1") {
	fm= lm(x~ I((1:length(x))/120))
	se.ols=summary(fm)$coef[2,2]
  	N=length(x)
	(trend.obs= fm$coef[2])# 1.116659 
  	fm.arima=arima(fm$residuals,order=c(1,0,0));fm.arima
	r=fm.arima$coef[1]
        neff= N * (1-r)/(1+r) ;neff # 102.2358 
  	se.obs= sqrt((N-2)/(neff-2))* se.ols; se.obs # 0.071758 
   	trend= round(c(trend=unlist(trend.obs), neff=neff, se=se.obs,t=trend.obs/se.obs),4)
	names(trend)=c("trend","neff","se","t")
	trend
}}


read.fwf2=function(working,widths,colClasses="character") {
	widths1=c(0,cumsum(widths));K=length(widths1)
	widths1=cbind(widths1[1:(K-1)]+1,widths1[2:K])
	K=K-1
	N=length(working)
	X= array(NA,dim=c(N,K)); X=data.frame(X)
	for(i in 1:K) X[,i]=substr(working,widths1[i,1],widths1[i,2])
	temp=colClasses=="numeric"
	for (i in (1:K)[temp]) X[,i]=as.numeric(X[,i])
	X[,temp][X[,temp]==na_char]=NA
	for (i in (1:K)[!temp]) X[,i]=gsub(" +$","",X[,i])
	return(X)
	}



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

     extractf=function(x,units=1,tol=2) {
 	r1=range(time(x))
	year=c(time(x))
 	y=diff(x)
 	z=round(diff(y),tol)
 	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}