vers <- 'VERSION: CADSWES Study Manager, Stochastic Hydrology Simulator R Scripts, Andrew Verdin, 5-15-2012'
library(zoo)
library(grid)
library(reshape)
library(ggplot2)
library(msm)
library(e1071)


is.wholenumber <- function(x, tol = .Machine$double.eps^0.5){
		abs(x - round(x)) < tol
	}

v_cc_knn <- function(xpath,sim.len=NULL,nsims=NULL,filepath=NULL,output.rdata=NULL,mean.change=NULL,returnpath=NULL,filename=NULL) {  
# 4/27/2012 4:45PM
## mean.change should be a decimal between 0 and 1.  
############# the decimal corresponding to the % mean change 0.2 = 20% mean reduction
## x should be a time series of streamflows 
if(is.character(xpath)) {
	x <- scan(xpath)
}else{
	x <- xpath
}


if(is.null(filepath)) filepath = paste(getwd(),'/',sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(is.null(returnpath)) returnpath = paste(filepath,"/successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: K-Nearest Neighbor",file=returnpath,append=TRUE)
write(paste("PARAMETER: historic streamflows: ",xpath,sep=""),file=returnpath,append=TRUE)
	plot.dir <- paste(filepath,'/plots',sep="")
	output.dir <- paste(filepath,'/output',sep="")
if(is.null(sim.len)) sim.len = length(x)
write(paste("PARAMETER: trace length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate trace length value",file=returnpath,append=TRUE)
	stop("Inappropriate trace length value",call.=FALSE)
}
if(is.null(nsims)) nsims = 1200
write(paste("PARAMETER: trace count: ",nsims,sep=""),file=returnpath,append=TRUE)
if(nsims < 1 | nsims > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate trace count value",file=returnpath,append=TRUE)
	stop("Inappropriate trace count value",call.=FALSE)
}
if(is.null(output.rdata)) {
	if(is.null(mean.change)) {
		output.rdata = '/knn-sim.RData'
	}else{
		output.rdata = '/knn-cc-sim.RData'
	}
}
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: output file: ",output.rdata,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: mean change: ",mean.change,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)
  
	sims.knn <- matrix(NA,nrow=nsims,ncol=sim.len)

	message('Simulating...')
	pb <- txtProgressBar(1,nsims,style=3)
 
if(is.null(mean.change)){ 
for(i in 1:nsims){
    setTxtProgressBar(pb,i)
	sims.knn[i,] <- knn.sim(x,sim.len)
  }
}else{
for(i in 1:nsims){
	red <- seq(1,1+mean.change,length.out=sim.len)
    setTxtProgressBar(pb,i)
	sims.knn[i,] <- knn.sim(x,sim.len)*red
  }
}

  
	close(pb)
	save(sims.knn,file=paste(filepath,output.rdata,sep=""))
	write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)

traces <- rep("trace",nsims)
counts <- 1:nsims
traceFolders <- paste(traces,counts,sep="")

for (av in 1:length(traceFolders)){
dir.create(paste(filepath,traceFolders[av],sep="/"))
write.table(sims.knn[av,], file=paste(filepath,traceFolders[av],"/",filename,sep=""), row.names = F, col.names = F)
write(paste("OUTPUT: traces file: ",filepath,traceFolders[av],"/",filename,sep=""),file=returnpath,append=TRUE)
}



	
if(!file.exists(plot.dir)) dir.create(plot.dir)

##################################
# PDF of traces AND observed
##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/allsim-pdf-knn.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'allsim-pdf-knn.pdf'), height=4, width = 6)
	d <- density(sims.knn[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims.knn)){
		d <- density(sims.knn[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/allsim-pdf-knn-cc.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'allsim-pdf-knn-cc.pdf'), height=4, width = 6)
	d <- density(sims.knn[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims.knn)){
		d <- density(sims.knn[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}

	##################################
	# Basic Statistics KNN
	##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/basic-stats-knn.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'basic-stats-knn.pdf'), height=4, width = 6)
	stats.knn <- sim.stats.basic(sims.knn,x)
	boxplot.stats.basic(stats.knn,x)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/basic-stats-knn-cc.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'basic-stats-knn-cc.pdf'), height=4, width = 6)
	stats.knn <- sim.stats.basic(sims.knn,x)
	boxplot.stats.basic(stats.knn,x)
	dev.off()
}
	##################################
	# Drought Statistics KNN
	##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/drought-stats-knn.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'drought-stats-knn.pdf'), height=4, width = 6)
	stats.knn.dr <- sim.stats.drought(sims.knn)
	boxplot.stats.drought(stats.knn.dr,x)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/drought-stats-knn-cc.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,'drought-stats-knn-cc.pdf'), height=4, width = 6)
	stats.knn.dr <- sim.stats.drought(sims.knn)
	boxplot.stats.drought(stats.knn.dr,x)
	dev.off()
}
write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}

v_cc_hg <- function(xpath,paleopath,m=NULL,sim.len=NULL,nsims=NULL,thresh=NULL,threshclass=NULL,sampleclass='conditional',filepath=NULL,output.rdata=NULL,mean.change=NULL,returnpath=NULL,filename=NULL) {   
# 4/11/2012 2:25PM
## mean.change should be a decimal between 0 and 1.  
############# the decimal corresponding to the % mean change 0.2 = 20% mean reduction
## x should be a time series of streamflows 

if(is.character(xpath)) {
	x <- scan(xpath)
}else{
	x <- xpath
}

if(is.character(paleopath)) {
	paleo <- scan(paleopath)
}else{
	paleo <- paleopath
}

if(is.null(filepath)) filepath = paste(getwd(),'/',sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(is.null(returnpath)) returnpath = paste(filepath,"/successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: Homogeneous Markov Chain",file=returnpath,append=TRUE)
write(paste("PARAMETER: historic streamflows: ",xpath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: paleo streamflows: ",paleopath,sep=""),file=returnpath,append=TRUE)
	plot.dir <- paste(filepath,'/plots',sep="")
	output.dir <- paste(filepath,'/output',sep="")
if(is.null(m)) m = 3
tpmtest <- tpm(paleo,m)
write(paste("PARAMETER: states count: ",m,sep=""),file=returnpath,append=TRUE)
if(m < 2 | m > 3) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Specify 2 or 3 states only",file=returnpath,append=TRUE)
	stop("Specify 2 or 3 states only",call.=FALSE)
}
if(is.null(sim.len)) sim.len = length(x)
write(paste("PARAMETER: trace length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 100000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate sim.len value",file=returnpath,append=TRUE)
	stop("Inappropriate sim.len value",call.=FALSE)
}
if(is.null(nsims)) nsims = 1200
write(paste("PARAMETER: trace count: ",nsims,sep=""),file=returnpath,append=TRUE)
if(nsims < 1 | nsims > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate nsims value",file=returnpath,append=TRUE)
	stop("Inappropriate nsims value",call.=FALSE)
}
if(is.null(output.rdata)) {
	if(is.null(mean.change)){
		output.rdata = paste('hg-sim-',m,'-state.RData',sep="")
	}else{
		output.rdata = paste('hg-cc-sim-',m,'-state.RData',sep="")
	}
}

write(paste("PARAMETER: thresholds class: ",threshclass,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: thresholds: ",thresh,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: output file: ",output.rdata,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: mean change: ",mean.change,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)  

	
		# these are all the values in a particular state 
	pools.to <- vector('list',m)
	if(is.null(thresh)) {
		for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m)==(i-1))
		tpmat <- tpm(paleo,m,limType='prob')
	}
	if(!is.null(thresh)) {
		if(length(thresh)!=(m-1)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Number of states does not match number of threshold values",file=returnpath,append=TRUE)
			write("i.e. 3 states -> 2 threshold values",file=returnpath,append=TRUE)
			stop("Number of states does not match number of threshold values (i.e. 3 states -> 2 threshold values")
		}
		if(min(thresh)<0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(threshclass=="raw") {
		if(min(thresh) < min(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) > max(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='raw',rawvals=thresh)==(i-1))
			qnts <- integer(length(thresh))
			for (i in 1:length(thresh)) qnts[i] <- length(x[x<thresh[i]])/length(x)
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=qnts)
		}
		if(threshclass=="quantile") {
		if(min(thresh) <= 0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) >= 1) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='quantile',quants=thresh)==(i-1))
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=thresh)
		}
		}
		
	#if(min(tpmat)<0.01) {
	#	write("ADVISE ADJUSTING THRESHOLD VALUES.  
	
	sims <- sims.yr <- matrix(NA,nrow=nsims,ncol=sim.len)

	message('Simulating...')
	pb <- txtProgressBar(1,nsims,style=3)
  
  if(is.null(mean.change)){

for(i in 1:nsims){
    setTxtProgressBar(pb,i)
    ind <- simulate.mc.paleo(x, tpm(paleo,m), pools.to, sim.len, weight=sampleclass, index = TRUE)
    sims[i,] <- x[ind]
    sims.yr[i,] <- time(x)[ind]
  }
  }else{
for(i in 1:nsims){
	red <- seq(1,1+mean.change,length.out=sim.len)
    setTxtProgressBar(pb,i)
    ind <- simulate.mc.paleo(x, tpm(paleo,m), pools.to, sim.len, weight=sampleclass, index = TRUE)
    sims[i,] <- x[ind]*red
    sims.yr[i,] <- time(x)[ind]
}
}
  



	close(pb)
	save(sims,sims.yr,file=paste(filepath,output.rdata,sep=""))

	write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)

traces <- rep("trace",nsims)
counts <- 1:nsims
traceFolders <- paste(traces,counts)

for (av in 1:length(traceFolders)){
dir.create(paste(filepath,traceFolders[av],sep="/"))
write.table(sims[av,],file=paste(filepath,"/",traceFolders[av],"/",filename,sep=''), row.names = F, col.names = F)
write(paste("OUTPUT: traces file: ",filepath,"/",traceFolders[av],"/",filename,sep=""),file=returnpath,append=TRUE)
write.table(sims.yr[av,],file=paste(filepath,"/",traceFolders[av],"/",filename,".Years",sep=''), row.names = F, col.names = F)
write(paste("OUTPUT: years file: ",filepath,"/",traceFolders[av],"/",filename,".Years",sep=""),file=returnpath,append=TRUE)
}	

if(!file.exists(plot.dir)) dir.create(plot.dir)
##################################
# PDF of traces AND observed
##################################
if(is.null(mean.change)){
	write(paste("PLOTS: ",plot.dir,"/allsim-pdf-hg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/allsim-pdf-hg-',m,'-state.pdf',sep="")), height=4, width = 6)
	d <- density(sims[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims)){
		d <- density(sims[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}else{
	write(paste("PLOTS: ",plot.dir,"/allsim-pdf-hg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/allsim-pdf-hg-cc-',m,'-state.pdf',sep="")), height=4, width = 6)
	d <- density(sims[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims)){
		d <- density(sims[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}

##################################
# Basic Statistics
##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/basic-stats-hg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/basic-stats-hg-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats <- sim.stats.basic(sims,x)
	boxplot.stats.basic(stats,x)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/basic-stats-hg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/basic-stats-hg-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats <- sim.stats.basic(sims,x)
	boxplot.stats.basic(stats,x)
	dev.off()
}
##################################
# Drought Statistics
##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/drought-stats-hg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/drought-stats-hg-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.dr <- sim.stats.drought(sims)
	boxplot.stats.drought(stats.dr,x)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/drought-stats-hg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/drought-stats-hg-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.dr <- sim.stats.drought(sims)
	boxplot.stats.drought(stats.dr,x)
	dev.off()
}
##################################
# Simulation Transition Prob boxplots
##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/tp-hg-sim-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/tp-hg-sim-',m,'-state.pdf',sep='')), height=4, width = 6)
	sim.tpms <- lapply(as.data.frame(t(sims)),tpm,m)
	sim.trans <- vector('list',m^2)
	for(i in 1:m^2)
		sim.trans[[i]] <- as.vector(sapply(sim.tpms,'[',i))
  
	boxplot(sim.trans,cex=.5)
	points(as.vector(tpm(paleo,m)),col=2,pch=2)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/tp-hg-sim-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/tp-hg-sim-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	sim.tpms <- lapply(as.data.frame(t(sims)),tpm,m)
	sim.trans <- vector('list',m^2)
	for(i in 1:m^2)
		sim.trans[[i]] <- as.vector(sapply(sim.tpms,'[',i))
  
	boxplot(sim.trans,cex=.5)
	points(as.vector(tpm(paleo,m)),col=2,pch=2)
	dev.off()
}
##################################
# Paleo Transition Prob boxplots
##################################
if(is.null(mean.change)){
write(paste("PLOTS: ",plot.dir,"/tp-hg-paleo-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/tp-hg-paleo-',m,'-state.pdf',sep='')), height=4, width = 6)
	paleo.tpms <- vector('list',length(paleo)-sim.len)
	for(i in 1:(length(paleo)-sim.len))
		paleo.tpms[[i]] <- tpm(paleo[i:(i+sim.len-1)],m)
  
	paleo.trans <- vector('list',m^2)
	for(i in 1:m^2)
		paleo.trans[[i]] <- as.vector(sapply(paleo.tpms,'[',i))

	boxplot(paleo.trans,cex=.5)
	points(as.vector(tpm(paleo,m)),col=2,pch=2)
	dev.off()
}else{
write(paste("PLOTS: ",plot.dir,"/tp-hg-paleo-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/tp-hg-paleo-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	paleo.tpms <- vector('list',length(paleo)-sim.len)
	for(i in 1:(length(paleo)-sim.len))
		paleo.tpms[[i]] <- tpm(paleo[i:(i+sim.len-1)],m)
  
	paleo.trans <- vector('list',m^2)
	for(i in 1:m^2)
		paleo.trans[[i]] <- as.vector(sapply(paleo.tpms,'[',i))

	boxplot(paleo.trans,cex=.5)
	points(as.vector(tpm(paleo,m)),col=2,pch=2)
	dev.off()
}
##################################
# Paleo Transition Prob time series
##################################
write(paste("PLOTS: ",plot.dir,"/tp-ts-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/tp-ts-',m,'-state.pdf',sep='')))
	grid.newpage() 
	pushViewport(viewport(layout = grid.layout(m, 1))) 
	vplayout <- function(x, y) 
      viewport(layout.pos.row = x, layout.pos.col = y)

	for(i in 1:m){
		paleo.trans.df <- as.data.frame(paleo.trans[1:m+((i-1)*m)])
		rownames(paleo.trans.df) <- -1:(end(ts(paleo,-1))[1] - sim.len)
		names(paleo.trans.df) <- paste(i,1:m,sep='')
		paleo.df <- melt(as.matrix(paleo.trans.df))
		paleo.df$X2 <- factor(paleo.df$X2)
		names(paleo.df) <- c('Year', 'Transition', 'Transition Probability')
		p <- qplot(Year,`Transition Probability`,data=paleo.df,color=Transition,geom='line')
		p <- p + theme_bw() +scale_color_manual(values= gray.colors(9))
		print(p, vp = vplayout(i, 1 ))
	}
	dev.off()

write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}

v_cc_nhg <- function(xpath,paleopath,m=NULL,sim.len=NULL,nsims=NULL,thresh=NULL,threshclass=NULL,sampleclass='conditional',filepath=NULL,output.rdata=NULL,ccPath=NULL,high.var=FALSE,low.var=FALSE,high.mean=FALSE,low.mean=FALSE,mean.change=NULL,returnpath=NULL,filename=NULL) {   
# 5/7/2012 8:45AM
## mean.change should be a decimal between 0 and 1.  
############# the decimal corresponding to the % mean change 0.2 = 20% mean reduction
## x should be a time series of streamflows 
## read the paleo data (tree rings!)
if(is.character(xpath)) {
	x <- scan(xpath)
}else{
	x <- xpath
}
if(is.character(paleopath)){
	paleo <- scan(paleopath)
}else{
	paleo <- paleopath
}

if(is.null(filepath)) filepath = paste(getwd(),'/',sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(is.null(returnpath)) returnpath = paste(filepath,"/successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: Non-Homogeneous Markov Chain",file=returnpath,append=TRUE)
write(paste("PARAMETER: historic streamflows: ",xpath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: paleo streamflows: ",paleopath,sep=""),file=returnpath,append=TRUE)
	plot.dir <- paste(filepath,'/plots',sep="")
	output.dir <- paste(filepath,'/output',sep="")
if(is.null(m) && is.null(thresh)) m = 2
if(is.null(m) && !is.null(thresh)) m = length(thresh)+1
write(paste("PARAMETER: states count: ",m,sep=""),file=returnpath,append=TRUE)
if(is.null(sim.len)) sim.len = length(x)
write(paste("PARAMETER: traces length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 100000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate sim.len value",file=returnpath,append=TRUE)
	stop("Inappropriate sim.len value",call.=FALSE)
}
if(is.null(nsims)) nsims = 1200
write(paste("PARAMETER: traces count: ",nsims,sep=""),file=returnpath,append=TRUE)
if(nsims < 1 | nsims > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate nsims value",file=returnpath,append=TRUE)
	stop("Inappropriate nsims value",call.=FALSE)
}
if(is.null(output.rdata)) {
	if(high.var==FALSE && is.null(mean.change)){
		output.rdata = paste('nhg-sim-',m,'-state.RData',sep="")
	}else{
		output.rdata = paste('nhg-cc-sim-',m,'-state.RData',sep="")
	}
}

write(paste("PARAMETER: thresholds class: ",threshclass,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: thresholds: ",thresh,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: thresholds class: ",threshclass,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: output file: ",output.rdata,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: paleo variance weighting: ",high.var,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: mean change: ",mean.change,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)  


	# these are all the values in a particular state 
	pools.to <- vector('list',m)
	if(is.null(thresh)) {
		for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m)==(i-1))
		tpmat <- tpm(paleo,m,limType='prob')
	}
	if(!is.null(thresh)) {
		if(length(thresh)!=(m-1)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Number of states does not match number of threshold values",file=returnpath,append=TRUE)
			write("i.e. 3 states -> 2 threshold values",file=returnpath,append=TRUE)
			stop("Number of states does not match number of threshold values (i.e. 3 states -> 2 threshold values")
		}
		if(min(thresh)<0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(threshclass=="raw") {
		if(min(thresh) < min(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) > max(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='raw',rawvals=thresh)==(i-1))
			qnts <- integer(length(thresh))
			for (i in 1:length(thresh)) qnts[i] <- length(x[x<thresh[i]])/length(x)
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=qnts)
		}
		if(threshclass=="quantile") {
		if(min(thresh) <= 0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) >= 1) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='quantile',quants=thresh)==(i-1))
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=thresh)
		}
		}
		
	if(min(tpmat)==0) { 
	write("Error: One or more 3-state transition probabilities equal to zero (0)",file=returnpath,append=TRUE)
	write("Continuing with 2-state...Threshold now set at median",file=returnpath,append=TRUE)	
	m=2
		# recalculate the values for two-states
	pools.to <- vector('list',m)
	if(is.null(thresh)) {
		for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m)==(i-1))
		tpmat <- tpm(paleo,m,limType='prob')
	}
	if(!is.null(thresh)) {
		if(length(thresh)!=(m-1)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Number of states does not match number of threshold values",file=returnpath,append=TRUE)
			write("i.e. 3 states -> 2 threshold values",file=returnpath,append=TRUE)
			stop("Number of states does not match number of threshold values (i.e. 3 states -> 2 threshold values")
		}
		if(min(thresh)<0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(threshclass=="raw") {
		if(min(thresh) < min(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) > max(x)) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='raw',rawvals=thresh)==(i-1))
			qnts <- integer(length(thresh))
			for (i in 1:length(thresh)) qnts[i] <- length(x[x<thresh[i]])/length(x)
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=qnts)
		}
		if(threshclass=="quantile") {
		if(min(thresh) <= 0) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate lower threshold value",file=returnpath,append=TRUE)
			stop("Inappropriate lower threshold value")
		}
		if(max(thresh) >= 1) {
			write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
			write("Error: Inappropriate upper threshold value")
			stop("Inappropriate upper threshold value")
		}
			for(i in 1:m) pools.to[[i]] <- which(ntile.ts(x,m,limit.type='quantile',quants=thresh)==(i-1))
			tpmat <- tpm(paleo,m,limType='quantile',quantiles=thresh)
		}
		}
	}
	
	for(pool in 1:m) {
		if(length(pools.to[[pool]]) < sqrt(x)) {
			write("WARNING: The number of observations is dangerously small in one or more states",file=returnpath,append=TRUE)
			write("WARNING: Advised altering thresholds or continue with 2-states",file=returnpath,append=TRUE)
		}
	}
	
	
	sims.nh <- sims.nh.yr <- matrix(NA,nrow=nsims,ncol=sim.len)

	message('Simulating...')
	pb <- txtProgressBar(1,nsims,style=3)
  
for(i in 1:nsims){
    setTxtProgressBar(pb,i)
	if(is.null(ccPath)){
		window <- sample(1:(length(paleo)- sim.len),1)
	}else{
		if(high.var==TRUE & low.var==FALSE & high.mean==FALSE & low.mean==FALSE){
			if(is.character(ccPath)){
				ccx <- scan(ccPath)
				}else{
				ccx <- ccPath
				}
		p_var <- rollapply(paleo,width=sim.len,FUN=var)
		p_var <- p_var[1:(length(p_var)-1)]
		N <- length(p_var)
		index <- wts <- 1:N
		index <- rev(index[order(p_var)])
		wts <- 1/wts
		wts <- wts/sum(wts)
		window <- sample(1:(length(paleo)- sim.len),1,prob=wts)
			window <- window:(window+sim.len-1)
			ind <- simulate.mc.paleo(ccx, tpm(paleo[window],m), pools.to, sim.len, weight=sampleclass, index = TRUE)
				if(is.null(mean.change)){
					sims.nh[i,] <- ccx[ind]
				}else{
					red <- seq(1,1+mean.change,length.out=sim.len)
					sims.nh[i,] <- ccx[ind]*red
				}
					sims.nh.yr[i,] <- time(ccx)[ind]
		}
		}
		if(high.var==FALSE & low.var==TRUE & high.mean==FALSE & low.mean==FALSE){
			if(is.character(ccPath)){
				ccx <- scan(ccPath)
				}else{
				ccx <- ccPath
				}
		p_var <- rollapply(paleo,width=sim.len,FUN=var)
		p_var <- p_var[1:(length(p_var)-1)]
		N <- length(p_var)
		index <- wts <- 1:N
		index <- index[order(p_var)]
		wts <- 1/wts
		wts <- wts/sum(wts)
		window <- sample(1:(length(paleo)- sim.len),1,prob=wts)
			window <- window:(window+sim.len-1)
			ind <- simulate.mc.paleo(ccx, tpm(paleo[window],m), pools.to, sim.len, weight=sampleclass, index = TRUE)
				if(is.null(mean.change)){
					sims.nh[i,] <- ccx[ind]
				}else{
					red <- seq(1,1+mean.change,length.out=sim.len)
					sims.nh[i,] <- ccx[ind]*red
				}
					sims.nh.yr[i,] <- time(ccx)[ind]
		}
		
		if(high.var==FALSE & low.var==FALSE & high.mean==TRUE & low.mean==FALSE){
			if(is.character(ccPath)){
				ccx <- scan(ccPath)
				}else{
				ccx <- ccPath
				}
		p_var <- rollapply(paleo,width=sim.len,FUN=mean)
		p_var <- p_var[1:(length(p_var)-1)]
		N <- length(p_var)
		index <- wts <- 1:N
		index <- rev(index[order(p_var)])
		wts <- 1/wts
		wts <- wts/sum(wts)
		window <- sample(1:(length(paleo)- sim.len),1,prob=wts)
			window <- window:(window+sim.len-1)
			ind <- simulate.mc.paleo(ccx, tpm(paleo[window],m), pools.to, sim.len, weight=sampleclass, index = TRUE)
				if(is.null(mean.change)){
					sims.nh[i,] <- ccx[ind]
				}else{
					red <- seq(1,1+mean.change,length.out=sim.len)
					sims.nh[i,] <- ccx[ind]*red
				}
					sims.nh.yr[i,] <- time(ccx)[ind]
		}
		
		if(high.var==FALSE & low.var==FALSE & high.mean==FALSE & low.mean==TRUE){
			if(is.character(ccPath)){
				ccx <- scan(ccPath)
				}else{
				ccx <- ccPath
				}
		p_var <- rollapply(paleo,width=sim.len,FUN=mean)
		p_var <- p_var[1:(length(p_var)-1)]
		N <- length(p_var)
		index <- wts <- 1:N
		index <- index[order(p_var)]
		wts <- 1/wts
		wts <- wts/sum(wts)
		window <- sample(1:(length(paleo)- sim.len),1,prob=wts)
		    window <- window:(window+sim.len-1)
			ind <- simulate.mc.paleo(ccx, tpm(paleo[window],m), pools.to, sim.len, weight=sampleclass, index = TRUE)
				if(is.null(mean.change)){
					sims.nh[i,] <- ccx[ind]
				}else{
					red <- seq(1,1+mean.change,length.out=sim.len)
					sims.nh[i,] <- ccx[ind]*red
				}
					sims.nh.yr[i,] <- time(ccx)[ind]
		}
	
    window <- window:(window+sim.len-1)
	ind <- simulate.mc.paleo(x, tpm(paleo[window],m), pools.to, sim.len, weight=sampleclass, index = TRUE)
if(is.null(mean.change)){
    sims.nh[i,] <- x[ind]
}else{
	red <- seq(1,1+mean.change,length.out=sim.len)
	sims.nh[i,] <- x[ind]*red
}
    sims.nh.yr[i,] <- time(x)[ind]
 }
  
	close(pb)
	save(sims.nh,sims.nh.yr,file=paste(filepath,output.rdata,sep=""))
	write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)

traces <- rep("trace",nsims)
counts <- 1:nsims
traceFolders <- paste(traces,counts)

for (av in 1:length(traceFolders)){
dir.create(paste(filepath,traceFolders[av],sep="/"))
write.table(sims.nh[av,], file=paste(filepath,"/",traceFolders[av],"/",filename,sep=''), row.names = F, col.names = F)
write.table(sims.nh.yr[av,], file=paste(filepath,"/",traceFolders[av],"/",filename,".Years",sep=''), row.names = F, col.names = F)
write(paste("OUTPUT: traces file: ",filepath,"/",traceFolders[av],"/",filename,sep=""),file=returnpath,append=TRUE)
write(paste("OUTPUT: years file: ",filepath,"/",traceFolders[av],"/",filename,".Years",sep=""),file=returnpath,append=TRUE)
}

if(!file.exists(plot.dir)) dir.create(plot.dir)

##################################
# PDF of traces AND observed
##################################
if(is.null(mean.change)){
	write(paste("PLOTS: ",plot.dir,"/allsim-pdf-nhg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/allsim-pdf-nhg-',m,'-state.pdf',sep="")), height=4, width = 6)
	d <- density(sims.nh[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims.nh)){
		d <- density(sims.nh[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}else{
	write(paste("PLOTS: ",plot.dir,"/allsim-pdf-nhg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/allsim-pdf-nhg-cc-',m,'-state.pdf',sep="")), height=4, width = 6)
	d <- density(sims.nh[1,])
	plot(d,col="grey",main="PDF: simulated vs. observed",xlab="",ylab="",ylim=c(0,max(d$y)*2))
	for(den in 2:nrow(sims.nh)){
		d <- density(sims.nh[den,])
		lines(d,col="grey")
	}
	d <- density(x)
	lines(d,col="red")
	dev.off()
}

##################################
# Basic Statistics NH
##################################
if(high.var==FALSE && is.null(mean.change)){
	write(paste("PLOTS: ",plot.dir,"/basic-stats-nhg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/basic-stats-nhg-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.nh <- sim.stats.basic(sims.nh,x)
	boxplot.stats.basic(stats.nh,x)
	dev.off()
}else{
	write(paste("PLOTS: ",plot.dir,"/basic-stats-nhg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/basic-stats-nhg-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.nh <- sim.stats.basic(sims.nh,x)
	boxplot.stats.basic(stats.nh,x)
	dev.off()
}
##################################
# Drought Statistics NH
##################################
if(high.var==FALSE && is.null(mean.change)){
	write(paste("PLOTS: ",plot.dir,"/drought-stats-nhg-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/drought-stats-nhg-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.nh.dr <- sim.stats.drought(sims.nh)
	boxplot.stats.drought(stats.nh.dr,x)
	dev.off()
}else{
	write(paste("PLOTS: ",plot.dir,"/drought-stats-nhg-cc-",m,"-state.pdf",sep=""),file=returnpath,append=TRUE)
	pdf(file.path(plot.dir,paste('/drought-stats-nhg-cc-',m,'-state.pdf',sep='')), height=4, width = 6)
	stats.nh.dr <- sim.stats.drought(sims.nh)
	boxplot.stats.drought(stats.nh.dr,x)
	dev.off()
}
write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}


annual.to.monthly <- function(mntpath, first.year, ann=NULL, sim=NULL, sim.len=NULL, filepath=NULL, returnpath=NULL, filename=NULL){
# 4/11/2012 2:25PM
## mnt should be vector with length = nyear*12 
## this code only applicable to disag one site from annual to monthly
if(is.character(mntpath)) {
	mnt <- scan(mntpath)
}else{
	mnt <- mntpath
}

if(is.null(filepath)) filepath = paste(getwd(),"/",sep="")
plot.dir <- paste(filepath,"/plots/",sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(!file.exists(plot.dir)) dir.create(plot.dir)
if(is.null(returnpath)) returnpath = paste(filepath,"/successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: Temporal Disaggregation",file=returnpath,append=TRUE)
write(paste("PARAMETER: monthly flow: ",mntpath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: first year: ",first.year,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: annual flow: ",ann,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: simulated flow for disag: ",sim,sep=""),file=returnpath,append=TRUE)


# sum monthly to annual (unless user has input)
if(is.null(ann)){
an <- rep(NA,(length(mnt)/12))
m=1;n=12
for (i in 1:length(an)) {
	an[i] <- sum(mnt[m:n])
	m=m+12;n=n+12	}
	ann <- an
}else{
if(is.character(ann)){
an <- scan(ann)
ann <- an
}else{
an <- ann
ann <- an
}
}

# initialize month matrix (nrow = nyear, ncol = 12 (months))
mntmat <- prbmat <- matrix(NA,nrow=length(ann),ncol=12)
m=1;n=12
for (i in 1:length(ann)){
mntmat[i,] <- mnt[m:n]
m=m+12;n=n+12
}

# initialize probability matrix (dimensions match month matrix)
prbmat <- mntmat/ann

# set number of years to simulate
if(is.null(sim.len)) sim.len = length(ann)
write(paste("PARAMETER: traces length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate value for Trace Length",file=returnpath,append=TRUE)
	stop("Inappropriate value for Trace Length",call.=FALSE)
}
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)

# simulate annual streamflow data (unless user has input)
if(is.null(sim)){
annsim <- ar(ann-mean(ann))
sim <- arima.sim(list(ar=annsim$ar),n=sim.len,sd=sqrt(annsim$var))
sim <- sim + mean(ann)
sim[sim<0]=0
}else{
if(is.character(sim)){
sims <- scan(sim)
sim <- sims
}else{
sims <- sim
sim <- sims
}
}

# initialize disagg matrix and corresponding year matrix
years <- matrix(NA,nrow=sim.len,ncol=1)
dismat <- matrix(NA,nrow=sim.len*12,ncol=1)
dismat2 <- matrix(NA,nrow=sim.len,ncol=12)

# KNN selection process for first year of simulated data
k <- sqrt(length(ann))
flow <- sim[1]
d <- abs(ann-flow)
delta <- cbind(seq(first.year,by=1,length.out=length(ann)),d)
delta_sort <- cbind(delta[,1][order(delta[,2])], sort(delta[,2]))
kmatrix = delta_sort[1:k,1:2] #selects the "k-nearest-neighbors" from Delta_sort 
weight = matrix(nrow=k, ncol=1) # defines matrix for weights
rnk=rank(kmatrix[,2]) #ranks distances for purpose of generating weights
	for(i in 1:k){
		weight[i,1] = 1/(rnk[i]) #fills weighting matrix
	}
z = sum(weight) # sums weights 
weights = weight/z	#divides weights by sum of weights so cumulative probability = 1
N=sample(kmatrix[,1], 1, replace = TRUE, prob=weights) #Selects a year to be "nearest neighbor"
pos=N-(first.year-1) # index for selected yr

dismat[1:12,1] <- flow*prbmat[pos,]
dismat2[1,] <- flow*prbmat[pos,]
years[1] <- N

m=13
n=24
# now disagg the rest of the years of simulated data
for (ij in 2:sim.len){

flow <- sim[ij]
d <- abs(ann-flow)
delta <- cbind(seq(first.year,by=1,length.out=length(ann)),d)
delta_sort <- cbind(delta[,1][order(delta[,2])], sort(delta[,2]))
kmatrix = delta_sort[1:k,1:2] #selects the "k-nearest-neighbors" from Delta_sort 
weight = matrix(nrow=k, ncol=1) # defines matrix for weights
rnk=rank(kmatrix[,2]) #ranks distances for purpose of generating weights
	for(i in 1:k){
		weight[i,1] = 1/(rnk[i]) #fills weighting matrix
	}
z = sum(weight) # sums weights 
weights = weight/z	#divides weights by sum of weights so cumulative probability = 1
N=sample(kmatrix[,1], 1, replace = TRUE, prob=weights) #Selects a year to be "nearest neighbor"
pos=N-(first.year-1) # index for selected yr


dismat[m:n,1] <- flow*prbmat[pos,]
dismat2[ij,] <- flow*prbmat[pos,]
m=m+12
n=n+12
years[ij] <- N

}

if(is.null(filename))	filename <- 'site1'

write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)
write(paste("OUTPUT: disaggregated flat file (rows = Jan,Feb,etc): ",filepath,"/",filename,sep=""),file=returnpath,append=TRUE)
write.table(dismat,file=paste(filepath,filename,sep="/"),row.names=FALSE,col.names=FALSE)

write(paste("PLOTS: ",plot.dir,filename,".pdf",sep=""),file=returnpath,append=TRUE)
pdf(file=paste(plot.dir,filename,".pdf",sep=""))
boxplot(dismat2,main="Temporal Disaggregation",names=month.abb)
lines(apply(mntmat,2,mean),col="red")
dev.off()

write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}

spatial.disag <- function(mntpath=NULL, annpath=NULL, indsite, first.year=NULL, nsite, sim=NULL, sim.len=NULL, filepath=NULL, returnpath=NULL, filename=NULL){
# 4/25/2012 1:52PM
## mnt should be matrix with nrow = nyears*12, ncol = nsites 
## this code only applicable to disag multiple sites
if(!is.null(mntpath)){
if(is.character(mntpath)) {
	mnt <- matrix(scan(mntpath),ncol=nsite,byrow=T)
}else{
	mnt <- mntpath
}

if(is.null(filepath)) filepath = paste(getwd(),"/",sep="")
plot.dir <- paste(filepath,"/plots/",sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(!file.exists(plot.dir)) dir.create(plot.dir)
if(is.null(returnpath)) returnpath = paste(filepath,"successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: Spatial Disaggregation",file=returnpath,append=TRUE)
write(paste("PARAMETER: monthly flow: ",mntpath,sep=""),file=returnpath,append=TRUE)
#write(paste("PARAMETER: first year: ",first.year,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: number of sites: ",nsite,sep=""),file=returnpath,append=TRUE)
if(nsite < 1 | nsite > 10000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate value for Site Count",file=returnpath,append=TRUE)
	stop("Inappropriate value for Site Count",call.=FALSE)
}
write(paste("PARAMETER: index site: ",indsite,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: simulated flow for disag: ",sim,sep=""),file=returnpath,append=TRUE)

# determine number of sites
######## NSITE IS NOW REQUIRED
#if(is.null(nsite)) nsite=ncol(mnt)

# read in INDEX GAGE
if(is.character(indsite)){
isite <- scan(indsite)
}else{
isite <- indsite
}

# simulate streamflow data (unless user has input)
if(is.null(sim)){
indsim <- ar(isite-mean(isite))
sim <- arima.sim(list(ar=indsim$ar),n=sim.len,sd=sqrt(indsim$var))
sim <- sim + mean(isite)
sim[sim<0]=0
}else{
if(is.character(sim)){
sims <- scan(sim)
sim <- sims
}else{
sims <- sim
sim <- sims
}
}

# calculate spatial pcts
sp_prbmat <- matrix(NA,nrow=nrow(mnt),ncol=nsite)
sp_prbmat <- mnt/isite

# set length to simulate
sim.len = length(sim)
write(paste("PARAMETER: traces length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 100000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate sim.len value",file=returnpath,append=TRUE)
	stop("Inappropriate sim.len value",call.=FALSE)
}
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)

nyear = nrow(mnt)/12

# initialize disagg matrix and corresponding year matrix
years <- matrix(NA,nrow=sim.len,ncol=1)
dismat <- matrix(NA,nrow=sim.len*12,ncol=1)
dismat2 <- matrix(NA,nrow=sim.len,ncol=12)
disarray <- array(NA,dim=c(sim.len,12,nsite))
mntarray <- array(NA,dim=c(nyear,12,nsite))

sp_dismat <- matrix(NA,nrow=sim.len,ncol=nsite)

k <- sqrt(length(isite))
for (jk in 1:sim.len){

iflow <- sim[jk]
#print(iflow)
d <- abs(isite-iflow)
delta <- cbind(seq(1,by=1,length.out=length(isite)),d)
delta_sort <- cbind(delta[,1][order(delta[,2])], sort(delta[,2]))
kmatrix = delta_sort[1:k,1:2] #selects the "k-nearest-neighbors" from d_sort 
weight = matrix(nrow=k, ncol=1) # defines matrix for weights
rnk=rank(kmatrix[,2]) #ranks distances for purpose of generating weights
	for(i in 1:k){
		weight[i,1] = 1/(rnk[i]) #fills weighting matrix
	}
z = sum(weight) # sums weights 
weights = weight/z	#divides weights by sum of weights so cumulative probability = 1
N=sample(kmatrix[,1], 1, replace = TRUE, prob=weights) #Selects a value to be "nearest neighbor"
pos=N  # index for selected neighbor
#print(jk)

#print(iflow*sp_prbmat[pos,])
#print(as.numeric(iflow)*as.numeric(sp_prbmat[pos,]))
sp_dismat[jk,] <- as.numeric(iflow)*as.numeric(sp_prbmat[pos,])
}



for(bc in 1:nsite){
	m=1
	n=12
	for(ab in 1:(sim.len/12)){
		disarray[ab,,bc] <- sp_dismat[m:n,bc]
		m=m+12
		n=n+12
	}
	m=1
	n=12
	for(abc in 1:nyear){
		mntarray[abc,,bc] <- mnt[m:n,bc]
		m=m+12
		n=n+12
	}
}

if(is.null(filename)) {
sites <- rep('site',nsite)
filename <- paste(sites,1:nsite,sep="")
}

write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)
write(paste("OUTPUT: spatial disag flat file: ",filepath,"/",nsite,"_site_spatial_disag.txt",sep=""),file=returnpath,append=TRUE)
write.table(sp_dismat,file=paste(filepath,"/",nsite,"_site_spatial_disag.txt",sep=""),row.names=FALSE,col.names=FALSE)
#return(sp_dismat)

for(jj in 1:nsite){
write(paste("OUTPUT: spatial disag file for ",filename[jj],": ",filepath,"/",filename[jj],sep=""),file=returnpath,append=TRUE)
write.table(sp_dismat[,jj],file=paste(filepath,"/",filename[jj],sep=""),row.names=FALSE,col.names=FALSE)
}

for(jj in 1:nsite){
write(paste("PLOTS: ",plot.dir,filename[jj],".pdf",sep=""),file=returnpath,append=TRUE)
pdf(file=paste(plot.dir,filename[jj],".pdf",sep=""))
#boxplot.stats(disarray[,,jj],apply(mntarray[,,jj],2,mean),lc=2)
boxplot(disarray[,,jj],names=month.abb)
lines(apply(mntarray[,,jj],2,mean),col="red")
title(main=paste("Spatial Disag Check - ",filename[jj],sep=""), outer=TRUE, line=-2)
dev.off()
}

write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}


if(!is.null(annpath)){
if(is.character(annpath)) {
	ann <- matrix(scan(annpath),ncol=nsite,byrow=T)
}else{
	ann <- annpath
}

if(is.null(filepath)) filepath = paste(getwd(),"/",sep="")
plot.dir <- paste(filepath,"/plots/",sep="")
if(!file.exists(filepath)) dir.create(filepath)
if(!file.exists(plot.dir)) dir.create(plot.dir)
if(is.null(returnpath)) returnpath = paste(filepath,"/successReturnPath.txt",sep="")
write(x=vers,file=returnpath)
write(x=paste("STARTED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
write("FUNCTION: Spatial Disaggregation",file=returnpath,append=TRUE)
write(paste("PARAMETER: annual flow: ",annpath,sep=""),file=returnpath,append=TRUE)
#write(paste("PARAMETER: first year: ",first.year,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: number of sites: ",nsite,sep=""),file=returnpath,append=TRUE)
if(nsite < 1 | nsite > 100000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate value for Site Count",file=returnpath,append=TRUE)
	stop("Inappropriate value for Site Count",call.=FALSE)
}
write(paste("PARAMETER: index site: ",indsite,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: simulated flow for disag: ",sim,sep=""),file=returnpath,append=TRUE)

# determine number of sites
######## NSITE IS NOW REQUIRED
#if(is.null(nsite)) nsite=ncol(mnt)

# read in INDEX GAGE
if(is.character(indsite)){
isite <- scan(indsite)
}else{
isite <- indsite
}

# simulate streamflow data (unless user has input)
if(is.null(sim)){
indsim <- ar(isite-mean(isite))
sim <- arima.sim(list(ar=indsim$ar),n=sim.len,sd=sqrt(indsim$var))
sim <- sim + mean(isite)
sim[sim<0]=0
}else{
if(is.character(sim)){
sims <- scan(sim)
sim <- sims
}else{
sims <- sim
sim <- sims
}
}

# calculate spatial pcts
sp_prbmat <- matrix(NA,nrow=nrow(ann),ncol=nsite)
sp_prbmat <- ann/isite

# set length to simulate
sim.len = length(sim)
write(paste("PARAMETER: traces length: ",sim.len,sep=""),file=returnpath,append=TRUE)
if(sim.len < 1 | sim.len > 100000) {
	write("FUNCTION STATUS: FAILURE",file=returnpath,append=TRUE)
	write("Error: Inappropriate sim.len value",file=returnpath,append=TRUE)
	stop("Inappropriate sim.len value",call.=FALSE)
}
write(paste("PARAMETER: output directory: ",filepath,sep=""),file=returnpath,append=TRUE)
write(paste("PARAMETER: return file: ",returnpath,sep=""),file=returnpath,append=TRUE)

nyear = nrow(ann)

# initialize disagg matrix and corresponding year matrix
years <- matrix(NA,nrow=sim.len,ncol=1)
dismat <- matrix(NA,nrow=sim.len*12,ncol=1)
dismat2 <- matrix(NA,nrow=sim.len,ncol=12)
disarray <- array(NA,dim=c(sim.len,12,nsite))
mntarray <- array(NA,dim=c(nyear,12,nsite))

sp_dismat <- matrix(NA,nrow=sim.len,ncol=nsite)

k <- sqrt(length(isite))
for (jk in 1:sim.len){

iflow <- sim[jk]
#print(iflow)
d <- abs(isite-iflow)
delta <- cbind(seq(1,by=1,length.out=length(isite)),d)
delta_sort <- cbind(delta[,1][order(delta[,2])], sort(delta[,2]))
kmatrix = delta_sort[1:k,1:2] #selects the "k-nearest-neighbors" from d_sort 
weight = matrix(nrow=k, ncol=1) # defines matrix for weights
rnk=rank(kmatrix[,2]) #ranks distances for purpose of generating weights
	for(i in 1:k){
		weight[i,1] = 1/(rnk[i]) #fills weighting matrix
	}
z = sum(weight) # sums weights 
weights = weight/z	#divides weights by sum of weights so cumulative probability = 1
N=sample(kmatrix[,1], 1, replace = TRUE, prob=weights) #Selects a value to be "nearest neighbor"
pos=N  # index for selected neighbor
#print(jk)

#print(iflow*sp_prbmat[pos,])
#print(as.numeric(iflow)*as.numeric(sp_prbmat[pos,]))
sp_dismat[jk,] <- as.numeric(iflow)*as.numeric(sp_prbmat[pos,])
}


if(is.null(filename)) {
sites <- rep('site',nsite)
filename <- paste(sites,1:nsite,sep="")
}

write("FUNCTION STATUS: SUCCESS",file=returnpath,append=TRUE)
write(paste("OUTPUT: spatial disag flat file: ",filepath,"/",nsite,"_site_spatial_disag.txt",sep=""),file=returnpath,append=TRUE)
write.table(sp_dismat,file=paste(filepath,"/",nsite,"_site_spatial_disag.txt",sep=""),row.names=FALSE,col.names=FALSE)
#return(sp_dismat)

for(jj in 1:nsite){
write(paste("OUTPUT: spatial disag file for ",filename[jj],": ",filepath,"/",filename[jj],sep=""),file=returnpath,append=TRUE)
write.table(sp_dismat[,jj],file=paste(filepath,"/",filename[jj],sep=""),row.names=FALSE,col.names=FALSE)
}

for(jj in 1:nsite){
write(paste("PLOTS: ",plot.dir,filename[jj],".pdf",sep=""),file=returnpath,append=TRUE)
pdf(file=paste(plot.dir,filename[jj],".pdf",sep=""))
#boxplot.stats(disarray[,,jj],apply(mntarray[,,jj],2,mean),lc=2)
boxplot(disarray[,,jj],names=month.abb)
lines(apply(mntarray[,,jj],2,mean),col="red")
title(main=paste("Spatial Disag Check - ",filename[jj],sep=""), outer=TRUE, line=-2)
dev.off()
}

write(x=paste("COMPLETED: ",as.character(Sys.time()),sep=""),file=returnpath,append=TRUE)
}
}

binary.ts <- 
function(x, limit = median(x), tie = 1 ){
    # returns a binary representation of x, 1 above limit, 0 below
    #
    # Ex. 
    #
    # binary.ts(1:10)
    ## [1] 0 0 0 0 0 1 1 1 1 1
    # binary.ts(1:10,8)
    ## [1] 0 0 0 0 0 0 0 1 1 1
    # binary.ts(1:10,8,tie=0)
    ## [1] 0 0 0 0 0 0 0 0 1 1
    
    b <- integer(length(x))
    
    filter <- if(tie == 1) x >= limit else x > limit
    
    #only need to set the 1's because b is already 0's
    b[filter] <- 1L
    
    if(class(x) == 'ts') 
        return(ts(b,start=start(x),end=end(x))) 
    else 
        return(b)
}

ntile.ts <- 
function(x, n, limit.type = 'prob', quants = NULL, rawvals = NULL, tie = 1, altobs = NULL ){
    # returns an integer vector corresponding to n states broken by equal 
    # probability or equal distance
    #
    limit <- 
    if(limit.type == 'prob') 
      quantile(x,seq(0,1,1/n))
    else if(limit.type == 'equal')
      seq(min(x),max(x),by=diff(range(x))/n)
	else if(limit.type == 'quantile')
	  quantile(x,c(0,quants,1))
	else if(limit.type == 'raw')
	  c(min(x),rawvals,max(x))
	  

    if(!is.null(altobs)) limit <- quantile(altobs,seq(0,1,1/n))
    
    b <- integer(length(x))
    
    for(i in 1:(n+1)){
        filter <- 
        if(tie == 1) 
            x >= limit[i] & x <= limit[i+1]
        else 
            x > limit[i] & x <= limit[i+1]
    
        #only need to set the 1's because b is already 0's
        b[filter] <- as.integer(i-1)
    }
    
    if(class(x) == 'ts') 
        return(ts(b,start=start(x),end=end(x))) 
    else 
        return(b)
}

binarys <- function(i,align=F){
    # returns the binary representation of integer vector i (coerced) 
    # as a character vector
    #
    # Ex.
    #
    # binarys(1:10)
    ## [1] "1"    "10"   "11"   "100"  "101"  "110"  "111"  "1000" "1001" "1010"
    # binarys(1:10,align=T)
    ## [1] "0001" "0010" "0011" "0100" "0101" "0110" "0111" "1000" "1001" "1010"
    #
    
    bb <- function(i) if (i) paste(bb(i %/% 2), i %% 2, sep="") else ""
        
        #number of binary digits of largest number
    width <- ceiling(log(max(i)+1,base=2))
    
        #get the character vector
    i <- sapply(as.integer(i),bb)
    
    if(align)
        return(sprintf(paste('%0',width,'d',sep=''),as.integer(i)))
    else
        return(i)
}

rows.equal <- function(x, y){
    
    if(is.null(dim(x))) x <- rbind(x)
    which( as.logical( 
            apply(x, 1, all.equal, y)))
}

#' Returns a first order transition probability matrix 
#'  for the given data and number of states
#' 
#' @param x a vector 
#' @param ns the number of states for the tpm. If NULL, it is assumed that the
#'           series is already converted to states.
tpm <- function(x,ns=NULL,limType='equal',quantiles=NULL,rawvals=NULL){
    
    require(msm)
    
    if(is.null(ns)){
      ns <- max(x)
      states <- x
      if(length(unique(states)) > 20) stop('Too many states, specify a smaller number.')
    }else{
      states <- ntile.ts(x,ns,limit.type=limType,quants=quantiles,rawvals=rawvals)
    }
    
    st <- statetable.msm(state,data=list(state=states))
    st/apply(st,1,sum)
    
}

trprob <- function(lag,lead,ns = max(lead)+1){
    
    require(gtools)
    require(Hmisc)
    
    q <- ncol(lag)
    r <- ncol(lead)
    
    states.lag <- gtools::permutations(ns,q,0:(ns-1),repeats.allowed=T)
    states.lead <- gtools::permutations(ns,r,0:(ns-1),repeats.allowed=T)
    nfrom <- nrow(states.lag)
    nto <- nrow(states.lead)
    
    nc <- ns^(q+r)
    
    tpm <- list(
        from = matrix(NA,nrow=nc,ncol=q),
        to = matrix(NA,nrow=nc,ncol=r), 
        p = matrix(NA,nrow=nfrom,ncol=nto),
        cor = matrix(NA,nrow=nfrom,ncol=nto),
        pts = matrix(NA,nrow=nfrom,ncol=nto),
        sig = matrix(NA,nrow=nfrom,ncol=nto),
        pval = matrix(NA,nrow=nfrom,ncol=nto),
        gcv = matrix(NA,nrow=nfrom,ncol=nto),
        pools = list()
        )
        
    pb <- txtProgressBar(1,nrow(states.lag),style=3)
    
    n <- 0
    tpm$pools$to <- tpm$pools$from <- list()
    rownames(tpm$p) <- character(nfrom)
    colnames(tpm$p) <- character(nto)
    
    for(i in 1:nfrom){
        setTxtProgressBar(pb,i)
        
            #all rows starting in a particular state
        pool <- rows.equal(lag,states.lag[i,])
        tpm$pools$from[[i]] <- pool
        tpm$pools$to[[i]] <- list()
    
        attr(tpm$pools$from[[i]],'state.from') <- states.lag[i,]
        
        for(j in 1:nto){
          
            n <- n + 1
              # record the current to and from states
            tpm$from[n,] <- states.lag[i,]
            tpm$to[n,] <- states.lead[j,]
            
            #matches <- rows.equal(cbind(lead[pool,]), states.lead[j,])
            this.state <- rbind(states.lead[j,])
            
              # turn vectors (single state to) into matricies, 
              # rbind handles the case of one entry only in the pool, 
              # but multiple to states 
            this.pool <- if(r==1) cbind(lead[pool,]) else rbind(lead[pool,])
            
              # If we get this condition, something is wrong.
            if(ncol(this.state) != ncol(this.pool))browser()
            
              # find the matching rows
            matches <- 
            if(nrow(this.pool) == 0){
              # then we're empty
              integer(0)
            }else{
              find.matches(this.state, this.pool, maxmatch=nrow(lag))$matches
            }
            
              # find.matches returns a single 0 if there are no matches
          if(length(matches) == 1) if( matches == 0 ) matches <- integer(0)
          
            # record the transition probability
            possible <- nrow(this.pool)
            tpm$p[i,j] <- length( matches ) / possible
              # set the transition state names, this sets them way too much
              # but it doesnt hurt anything
            rownames(tpm$p)[i] <- paste(tpm$from[n,],collapse='')
            colnames(tpm$p)[j] <- paste(tpm$to[n,],collapse='')
                
                #save the indexes of the states from and to
            tpm$pools$to[[i]][[j]] <- as.vector(matches)
          attr(tpm$pools$to[[i]],'state.from') <- states.lag[i,]
          attr(tpm$pools$to[[i]][[j]],'state.to') <- states.lead[j,]
           #browser()
        }
    }
    close(pb)
    return(tpm)
    
}

binary.combos <- function(n){
    as.matrix( expand.grid( rep( list(c(T,F)), n ) ) )
}

sim.pqr <- function(tr.paleo, conditional.pool = FALSE){
    
    if(conditional.pool){
            #simulate transition
        rand <- runif(1)
        which.to <- rank(c(rand,cumsum(this.p)))[1]
        state.to <- tr.paleo$to[which.to,]


            #conditional pool, quantiles given state.to
        pool.from <- tr.paleo$pools$from[[ which.from[which.to] ]]
        pool.to <- tr.paleo$pools$to[[ which.from[which.to] ]]
        pool.from <- matrix(pool.from,,q)
        pool.to <- matrix(pool.to,,r)
    }
    
}

mylag <- 
function(x,lag,docor=FALSE){

    if(lag>length(x)) 
        warning("Lag is larger than input vector length, returning NA's") 

    if(lag<0)  lagn = c(rep(NA,abs(lag)),x[-(length(x):(length(x)+lag+1))])
    if(lag==0) lagn = x    
    if(lag>0)  lagn = c(x[-(1:lag)],rep(NA,lag))

    remove = !is.na(lagn) & !is.na(x)
    if(docor){
        return(cor(x[remove],lagn[remove]))
    }else{
        return(lagn)
    }
}

lagc <- 
function(x,lag=1){

    if(lag>length(x)) 
        warning("Lag is larger than input vector length, returning NA's") 

    if(lag<0)  lagn = c(rep(NA,abs(lag)),x[-(length(x):(length(x)+lag+1))])
    if(lag==0) lagn = x    
    if(lag>0)  lagn = c(x[-(1:lag)],rep(NA,lag))

    remove = !is.na(lagn) & !is.na(x)
    return(cor(x[remove],lagn[remove]))
}


oneoverk <- function(k){
    
    W <- numeric(k)
    for(j in 1:k)
        W[j] <- 1/j
    W <- W/sum(W)
    W
    
}

############################################################
############################################################
# sigcor:
#       Returns the minimum significant correlation at the 
#       given alpha value
sigcor <- function(n,alpha=.05){
 
    k <- qnorm(alpha)^2/(n-2)
    sqrt(k/(k+1))   
    
}

#WAVELET  1D Wavelet transform with optional singificance testing
#
#   [WAVE,PERIOD,SCALE,COI] = wavelet(Y,DT,PAD,DJ,S0,J1,MOTHER,PARAM)
#
#   Computes the wavelet transform of the vector Y (length N),
#   with sampling rate DT.
#
#   By default, the Morlet wavelet (k0=6) is used.
#   The wavelet basis is normalized to have total energy=1 at all scales.
#
#
# INPUTS:
#
#    Y = the time series of length N.
#    DT = amount of time between each Y value, i.e. the sampling time.
#
# OUTPUTS:
#
#    WAVE is the WAVELET transform of Y. This is a complex array
#    of dimensions (N,J1+1). FLOAT(WAVE) gives the WAVELET amplitude,
#    ATAN(IMAGINARY(WAVE),FLOAT(WAVE) gives the WAVELET phase.
#    The WAVELET power spectrum is ABS(WAVE)^2.
#    Its units are sigma^2 (the time series variance).
#
#
# OPTIONAL INPUTS:
# 
# *** Note *** setting any of the following to -1 will cause the default
#               value to be used.
#
#    PAD = if set to 1 (default is 0), pad time series with enough zeroes to get
#         N up to the next higher power of 2. This prevents wraparound
#         from the end of the time series to the beginning, and also
#         speeds up the FFT's used to do the wavelet transform.
#         This will not eliminate all edge effects (see COI below).
#
#    DJ = the spacing between discrete scales. Default is 0.25.
#         A smaller # will give better scale resolution, but be slower to plot.
#
#    S0 = the smallest scale of the wavelet.  Default is 2*DT.
#
#    J1 = the # of scales minus one. Scales range from S0 up to S0*2^(J1*DJ),
#        to give a total of (J1+1) scales. Default is J1 = (LOG2(N DT/S0))/DJ.
#
#    MOTHER = the mother wavelet function.
#             The choices are 'MORLET', 'PAUL', or 'DOG'
#
#    PARAM = the mother wavelet parameter.
#            For 'MORLET' this is k0 (wavenumber), default is 6.
#            For 'PAUL' this is m (order), default is 4.
#            For 'DOG' this is m (m-th derivative), default is 2.
#
#
# OPTIONAL OUTPUTS:
#
#    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
#           to the SCALEs.
#
#    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J1
#            where J1+1 is the total # of scales.
#
#    COI = if specified, then return the Cone-of-Influence, which is a vector
#        of N points that contains the maximum period of useful information
#        at that particular time.
#        Periods greater than this are subject to edge effects.
#        This can be used to plot COI lines on a contour plot by doing:
#
#              contour(time,log(period),log(power))
#              plot(time,log(coi),'k')
#
#----------------------------------------------------------------------------
#   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
#
#   This software may be used, copied, or redistributed as long as it is not
#   sold and this copyright notice is reproduced on each copy made. This
#   routine is provided as is without any express or implied warranties
#   whatsoever.
#
# Notice: Please acknowledge the use of the above software in any publications:
#    ``Wavelet software was provided by C. Torrence and G. Compo,
#      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
#
# Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
#            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
#
# Please send a copy of such publications to either C. Torrence or G. Compo:
#  Dr. Christopher Torrence               Dr. Gilbert P. Compo
#  Research Systems, Inc.                 Climate Diagnostics Center
#  4990 Pearl East Circle                 325 Broadway R/CDC1
#  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
#  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
#----------------------------------------------------------------------------
wavelet = function(Y, DT = 1, pad = 0, dj = 0.25, param = 6) {
    #Y is time series to be analyzed
    #DT is timestep for annual data, 1
    #pad data ? 0=F, 1=T
    #dj= spacing between discrete scales (.25)
    #param = wavenumber (6)
    s0 = 2 * DT
    #s0=DT
    n1 = length(Y)
    J1 = floor((log2(n1 * DT/s0))/dj)
    #....construct time series to analyze, pad if necessary
    x = Y - mean(Y)
    if (pad == 1) {
        base2 = trunc(log(n1)/log(2) + 0.4999)  # power of 2 nearest to N
        x = c(x, rep(0, 2^(base2 + 1) - n1))
    }
    n = length(x)
    #....construct wavenumber array used in transform [Eqn(5)]
    k = (1:trunc(n/2))
    k = k * ((2 * pi)/(n * DT))
    k = c(0, k, -rev(k[1:floor((n - 1)/2)]))
    #....compute FFT of the (padded) time series
    f = fft(x)  # [Eqn(3)]
    #....construct SCALE array & empty PERIOD & WAVE arrays
    scale = s0 * 2^((0:J1) * dj)
    period = scale
    wave = matrix(data = 0, ncol = n, nrow = J1 + 1)  # define the wavelet array
    wave = as.complex(wave)  # make it complex
    wave = matrix(data = wave, ncol = n, nrow = J1 + 1)
    # loop through all scales and compute transform
    for (a1 in 1:(J1 + 1)) {
        scl = scale[a1]
        #source('wave_bases_r.txt')
        nn = length(k)
        k0 = param
        expnt = -(scl * k - k0)^2/(2 * (k > 0))
        norm = sqrt(scl * k[2]) * (pi^(-0.25)) * sqrt(nn)  # total energy=N   [Eqn(7)]
        daughter = norm * exp(expnt)
        daughter = daughter * (k > 0)  # Heaviside step function
        fourier_factor = (4 * pi)/(k0 + sqrt(2 + k0^2))  # Scale-->Fourier [Sec.3h]
        coi = fourier_factor/sqrt(2)  # Cone-of-influence [Sec.3g]
        dofmin = 2  # Degrees of freedom
        out <- list(daughter = daughter, fourier_factor = fourier_factor, coi = coi, 
            dofmin = dofmin)
        daughter = out$daughter
        fourier_factor = out$fourier_factor
        coi = out$coi
        dofmin = out$dofmin
        wave[a1, ] = fft((f * daughter), inverse = TRUE)/(length(f * daughter))  # wavelet transform[Eqn(4)]
    }
    period = fourier_factor * scale
    coi = coi * c(1:(floor(n1 + 1)/2), rev(1:floor(n1/2))) * DT
    wave = wave[, 1:n1]  # get rid of padding before returning
    power = abs(wave)^2
    list(wave = wave, period = period, scale = scale, power = power, coi = coi)
    # end of code
} 

rl.count <- function(x,ns){
    
    rl <- rle(ntile.ts(x,ns))$lengths
    df <- as.data.frame(table(rl))
    count <- df$Freq
    lens <- as.numeric(levels(df$rl))
    
    return(list(rl=rl,len=lens,count=count))
}

rl.count.lower <- function(x,ns){
    
    ql <- quantile(x,1/ns)
    x[x > ql] <- NA
    x[!is.na(x)] <- 1
    rl <- rle(as.vector(x))
    df <- as.data.frame(table(rl$lengths[!is.na(rl$values)]))
    count <- df$Freq
    lens <- as.numeric(levels(df[,1]))
    
    return(list(rl=rl$lengths,len=lens,count=count))
}

rl.count.upper <- function(x,ns){
    
    qu <- quantile(x,1-1/ns)
    x[x < qu] <- NA
    x[!is.na(x)] <- 1
    rl <- rle(as.vector(x))
    df <- as.data.frame(table(rl$lengths[!is.na(rl$values)]))
    count <- df$Freq
    lens <- as.numeric(levels(df[,1]))
    
    return(list(rl=rl$lengths,len=lens,count=count))
}


boxplot.stats <- function(sims,x,lc=1){
    
    
    layout(rbind(1:5))
    boxplot(apply(sims,lc,mean))
    points(mean(x),col='Blue',pch=2)
    mtext("Mean",1,1)
    boxplot(apply(sims,lc,sd))
    points(sd(x),col='Blue',pch=2)
    mtext('SD',1,1)
    boxplot(apply(sims,lc,min))
    points(min(x),col='Blue',pch=2)
    mtext('Min',1,1)
    boxplot(apply(sims,lc,max))
    points(max(x),col='Blue',pch=2)
    mtext('Max',1,1)
    boxplot(apply(sims,lc,mylag,1,docor=T))
    points(mylag(x,1,docor=T),col='Blue',pch=2)
    mtext('Lag 1 AC',1,1)
    abline(h=)
    
}

boxplot.stats.basic <- function(stats1, x){
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:6))
    
    boxplot(stats1$mean)
    points(mean(x),col='Blue')
    mtext("Mean",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats1$sd)
    points(sd(x),col='Blue')
    mtext('SD',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats1$skew)
    points(skewness(x),col='Blue')
    mtext('Skewness',3,1)
    
    boxplot(stats1$min)
    points(min(x),col='Blue')
    mtext('Min',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats1$max)
    points(max(x),col='Blue')
    mtext('Max',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats1$lag1)
    points(mylag(x,1,docor=T),col='Blue')
    mtext('Lag 1 AC',3,1)
    
}

boxplot.stats.basic2 <- function(stats1, stats2, x, names){
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:6))
    
    boxplot(cbind(stats1$mean,stats2$mean),names=names)
    abline(h=mean(x),col='Blue')
    mtext("Mean",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$sd,stats2$sd),names=names)
    abline(h=sd(x),col='Blue')
    mtext('SD',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$skew,stats2$skew),names=names)
    abline(h=skewness(x),col='Blue')
    mtext('Skewness',3,1)
    
    boxplot(cbind(stats1$min,stats2$min),names=names)
    abline(h=min(x),col='Blue')
    mtext('Min',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$max,stats2$max),names=names)
    abline(h=max(x),col='Blue')
    mtext('Max',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$lag1,stats2$lag1),names=names)
    abline(h=mylag(x,1,docor=T),col='Blue')
    mtext('Lag 1 AC',3,1)
    
}

boxplot.stats.basic3 <- function(stats1, stats2, stats3, x, names ){
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:3,4:6))
    
    boxplot(cbind(stats1$mean,stats2$mean,stats3$mean),names=names)
    abline(h=mean(x),col='Blue')
    mtext("Mean",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$sd,stats2$sd,stats3$sd),names=names)
    abline(h=sd(x),col='Blue')
    mtext('SD',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$skew,stats2$skew,stats3$skew),names=names)
    abline(h=skewness(x),col='Blue')
    mtext('Skewness',3,1)
    
    boxplot(cbind(stats1$min,stats2$min,stats3$min),names=names)
    abline(h=min(x),col='Blue')
    mtext('Min',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$max,stats2$max,stats3$max),names=names)
    abline(h=max(x),col='Blue')
    mtext('Max',3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$lag1,stats2$lag1,stats3$lag1),names=names)
    abline(h=mylag(x,1,docor=T),col='Blue')
    mtext('Lag 1 AC',3,1)
    
}

boxplot.stats.drought <- function(stats,x){
    
    obs.stats <- stats.drought(x,median(x))
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:4))
    
    boxplot(stats$dri)
    points(obs.stats$avgdri,col='Blue',pch=2)
    mtext("Ave DRI",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats$max.runlen)
    points(obs.stats$maxrunlen,col='Blue',pch=2)
    mtext("Max Run Length",1,1)
    
    boxplot(stats$max.def)
    points(obs.stats$maxdeficit,col='Blue',pch=2)
    mtext("Max Deficit",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats$ave.def)
    points(obs.stats$avgdeficit,col='Blue',pch=2)
    mtext("Ave Deficit",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
}

boxplot.stats.drought3 <- function(stats1, stats2, stats3, x, names){
    
    obs.stats <- stats.drought(x,median(x))
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:4))
    
    boxplot(cbind(stats1$dri,stats2$dri,stats3$dri),names=names)
    abline(h=obs.stats$avgdri,col='Blue')
    mtext("Ave DRI",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$max.runlen,stats2$max.runlen,stats3$max.runlen),names=names)
    abline(h=obs.stats$maxrunlen,col='Blue')
    mtext("Max Run Length",3,1)
    
    boxplot(cbind(stats1$max.def,stats2$max.def,stats3$max.def),names=names)
    abline(h=obs.stats$maxdeficit,col='Blue')
    mtext("Max Deficit",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(cbind(stats1$ave.def,stats2$ave.def,stats3$ave.def),names=names)
    abline(h=obs.stats$avgdeficit,col='Blue')
    mtext("Ave Deficit",3,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
}

boxplot.stats.drought2 <- function(stats1,stats2,x){
    
    obs.stats <- stats.drought(x,median(x))
    
    par(mar=c(5,3.2,3,.5)+.1)
    layout(rbind(1:4))
    
    boxplot(stats$dri)
    boxplot(stats$dri)
    points(obs.stats$avgdri,col='Blue',pch=2)
    mtext("Ave DRI",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats$max.runlen)
    points(obs.stats$maxrunlen,col='Blue',pch=2)
    mtext("Max Run Length",1,1)
    
    boxplot(stats$max.def)
    points(obs.stats$maxdeficit,col='Blue',pch=2)
    mtext("Max Deficit",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
    boxplot(stats$ave.def)
    points(obs.stats$avgdeficit,col='Blue',pch=2)
    mtext("Ave Deficit",1,1)
    mtext("Flow Volume (KAF)",2,2,cex=.7)
    
}


matexp <- function(x,m){
    mat <- x
    for(i in 1:m)
        mat <- mat %*% x
    mat
}

count.run.lengths <- function(sims, count.max=25, ns=3){
    
    nsims <- nrow(sims)
    sim.len <- ncol(sims)
    runs.mat <- matrix(NA,ncol=sim.len,nrow=nsims)
    for(i in 1:nsims){
        runs <- rl.count(sims[i,],ns)
        runs.mat[i,1:length(runs$rl)] <- runs$rl
    }
    runs.mat   
}

sim.stats.basic <- function(sims, x, count.max=25, ns=3, wavelet.res = 115){
    
    require(e1071)
    
    stats <- list()
    
    stats$mean <- apply(sims,1,mean)
    stats$sd <- apply(sims,1,sd)
    stats$skew <- apply(sims,1,skewness)
    stats$min <- apply(sims,1,min)
    stats$max <- apply(sims,1,max)
    stats$lag1 <- apply(sims,1,mylag,1,docor=T)
    stats$lag2 <- apply(sims,1,mylag,2,docor=T)
    
    stats
}

sim.stats <- function(sims, x, count.max=25, ns=3, wavelet.res = 115){
    
    require(sm)
    require(e1071)
    
    wavelet.res <- length(wavelet(sims[1,], pad=1, dj=.05)$scale)
    nsims <- nrow(sims)
    sim.len <- ncol(sims)
    stats <- list()
    stats$runs <- stats$runs.upper <- stats$runs.lower <- 
      matrix(NA,ncol=sim.len,nrow=nsims)
    stats$pdf.x <- seq(from=0,to=max(x)*1.25,len=50)
    stats$pdf <- matrix(NA,ncol=length(stats$pdf.x),nrow=nsims)
    stats$spec <- stats$freq <- matrix(NA,ncol=length(spectrum(sims[1,],plot=F)$spec),nrow=nsims)
    stats$dens <- stats$dens.lower <- stats$dens.upper <- 
      matrix(NA,ncol=512,nrow=nsims)
    stats$dens.x <- stats$dens.x.lower <- stats$dens.x.upper <- numeric(512)
    stats$counts <- matrix(0,ncol=count.max,nrow=nsims)
    stats$pow <- matrix(0,ncol=wavelet.res,nrow=nsims)

    pb <- txtProgressBar(1,nsims,style=3)
    for(i in 1:nsims){
        setTxtProgressBar(pb,i)
        runs <- rl.count(sims[i,],ns)
        runs.upper <- rl.count.upper(sims[i,],ns)
        runs.lower <- rl.count.lower(sims[i,],ns)
        #browser()
        stats$runs[i,1:length(runs$rl)] <- runs$rl
        stats$runs.upper[i,1:length(runs.upper$rl)] <- runs.upper$rl
        stats$runs.lower[i,1:length(runs.lower$rl)] <- runs.lower$rl
        d <- density(stats$runs[i,],bw=1,from=1,to=count.max,na.rm=T)
        stats$dens[i,] <- d$y
        d <- density(stats$runs.lower[i,],bw=1,from=1,to=count.max,na.rm=T)
        stats$dens.lower[i,] <- d$y
        d <- density(stats$runs.upper[i,],bw=1,from=1,to=count.max,na.rm=T)
        stats$dens.upper[i,] <- d$y
        d.sm <- sm.density(sims[i,], display='none', eval.points = stats$pdf.x)
        stats$pdf[i,] <- d.sm$estimate
        stats$counts[i,runs$len] <- runs$count
        stats$pow[i,] <- apply(wavelet(sims[i,], pad=1, dj=.05)$power,1,sum)
        #browser()
        spec <- spectrum(sims[i,],plot=F)
        stats$spec[i,] <- spec$spec
        stats$freq[i,] <- spec$freq
    }
    close(pb)
    
    stats$dens.x <- d$x
    
    stats$mean <- apply(sims,1,mean)
    stats$sd <- apply(sims,1,sd)
    stats$skew <- apply(sims,1,skewness)
    stats$min <- apply(sims,1,min)
    stats$max <- apply(sims,1,max)
    stats$lag1 <- apply(sims,1,mylag,1,docor=T)
    stats$lag2 <- apply(sims,1,mylag,2,docor=T)
    
    stats
}

knn.sim <- function(x,sim.len=length(x),nn=floor(sqrt(length(x)))){
    
    W <- cumsum(oneoverk(nn))
    sim <- numeric(sim.len)
    ll <- laglead(x,1,1)
        # random srarting year
    this.neighbor <- sample(1:length(x),1)
    sim[1] <- x[this.neighbor]
    for(i in 2:sim.len){
            # neighbors to the current value
            # dont use the last value because there is no following value
        this.x <- ll$lag[!(ll$lag == sim[i-1])]
        xsample <- c(sim[i-1],this.x)
        neighbors <- order(as.matrix(dist(xsample))[,1])[-1][1:nn]
        
        r1 <- runif(1)
        this.neighbor <- neighbors[rank(c(r1,W))[1]]
        
            #choose following value to neighbor selected
        val <- xsample[this.neighbor]
        if(is.na(val))browser()
        sim[i] <- ll$lead[ll$lag == val]
    }
    sim
}

laglead <- 
function(x,q,r){

    ## laglead:
    ##   This function gives all the chunks of a vector of length q 
    ##   followed by chunks of length r. It gives back a list of two 
    ##   matricies, lag and lead, essentially the lag q matrix 
    ##   and the subsequent r values. 
    ##
    ##   Inputs:
    ##      x - A vetor 
    ##      q - lead length
    ##      r - lag length 
    ##
    ## Ex.
    ## 
    ##    laglead(1:10,3,2)
    ##   $x
    ##        [,1] [,2] [,3]
    ##   [1,]    1    2    3
    ##   [2,]    2    3    4
    ##   [3,]    3    4    5
    ##   [4,]    4    5    6
    ##   [5,]    5    6    7
    ##   [6,]    6    7    8
    ##   
    ##   $y
    ##        [,1] [,2]
    ##   [1,]    4    5
    ##   [2,]    5    6
    ##   [3,]    6    7
    ##   [4,]    7    8
    ##   [5,]    8    9
    ##   [6,]    9   10
    
    
    n <- length(x)

  lag <- matrix(NA, n-q-r+1, q)
    lead <- matrix(NA, n-q-r+1, r)

    for(i in 1:(n-q-r+1)){
        lag[i,] <- x[i:(i+q-1)]
        lead[i,] <- x[(i+q):(i+q+r-1)]
    }
    return(list(lag=lag,lead=lead))
}

statdist <- function(tpm){
    
    m <- nrow(tpm)
    ones <- rbind(rep(1,m))
    I <- diag(rep(1,m))
    U <- matrix(rep(1,m^2),m)
    as.vector(ones %*% solve(I - tpm + U))
    
}

simulate.mc.paleo <- function(x, tpm, pools.to, sim.len, index = TRUE, 
    weight = 'oneoverk', burn = 20, stationary = TRUE){
  
  delta <- 
    if(stationary)
      statdist(tpm)
    else
      delta <- c(1,rep(0,m-1))
  
  m <- nrow(tpm)
  
  cumsums <- tpm
  x.pools <- pools.to
  npool <- integer(m)
  for(i in 1:m){
    cumsums[i,] <- cumsum(tpm[i,])
    x.pools[[i]] <- x[pools.to[[i]]]
    npool[i] <- length(x.pools[[i]])
  }
  
  sim <- ind <- numeric(sim.len)
  
    # initial state
  from <- rank(c(runif(1),cumsum(delta)))[1]
  x.from <- sample(x[ntile.ts(x,m) == (from-1)],1)
  
  for(i in (burn+1):(sim.len+burn)){
  
    to <- rank(c(runif(1),cumsums[from,]))[1]
    totrans <- pools.to[[to]]-1
	totrans <- totrans[totrans > 0 & totrans < (length(x))]
	#print(totrans)
    np <- npool[to]
	nptrans <- length(x[totrans])
	
    x.to.ind <-
      if(weight == 'equal')
        knn(x.from,x.pools[[to]],k=np,W=cumsum(rep(1/np,np)), ret = 'ind')
      else if(weight == 'conditional')
		knn(x.from,x[totrans],k=nptrans,W=cumsum(rep(1/nptrans,nptrans)), ret = 'ind')
	  else
        knn(x.from,x.pools[[to]],k=np, ret = 'ind')
    
	  
	  
	  if(weight == 'conditional') {
		if(x.to.ind > (length(x)-1)) {
			x.to.ind <- knn(x.from,x.pools[[to]],k=np,W=cumsum(rep(1/np,np)),ret='ind')
			x.to <- x.pools[[to]][x.to.ind]
			}
		x.to <- x[x.to.ind+1]
	  }else{
		x.to <- x.pools[[to]][x.to.ind]
		}
    #print(x.to.ind)
	#print(x.to)
	ind[i - burn] <- which(x == x.to)
	sim[i - burn] <- x.to
    
    from <- to
    x.from <- x.to
    
  }
  if(index) return(ind) else return(sim)
  
}

knn <- function(val, x, n = 1, k = NULL, W = NULL, ret='val'){
  
  if(is.null(k)) k <- floor(sqrt(length(x)))
  if(is.null(W)) W <- cumsum(oneoverk(k))
  
  d <- sqrt((x-val)^2)
  neighbors <- order(d)[1:k]
  
  vals <- numeric(n)
  for(i in 1:n){
    rand <- runif(1)
    this.neighbor <- rank(c(rand,W))[1]
    vals[i] <- if(ret == 'val') x[neighbors][this.neighbor] else which(x == x[neighbors][this.neighbor])
  }
  vals
  
}


sim.stats.drought <- function(sims){
    
    nsim <- nrow(sims)
    #tpm <- lapply(as.data.frame(t(sims)),tpm,3)
    
    dri <- max.runlen <- max.def <- ave.def <- numeric(nsim)
    for (j in 1:nsim){
      z <- stats.drought(sims[j,],median(sims[j,]))
     dri[j] <- z$avgdri
     max.runlen[j] <- z$maxrunlen
     max.def[j] <- z$maxdeficit
     ave.def[j] <- z$avgdeficit
    }
    
    #list(tpm=tpm, 
    list(dri = dri,  max.runlen = max.runlen, max.def = max.def, 
        ave.def = ave.def)
}


stats.drought <- function(datavec, threshold) {
    #threshold -- threshold for defining drought
    #calculate spell length
    #y <- datavec - mean(datavec)
    y <- datavec - threshold
    n <- length(datavec)
    z <- rep(0, n)  #0 implies a surplus
    ysaved <- rep(0, n)  #save y values to calculate maximum deficit
    dryr = 0  #counter for number of dry/below average years
    for (i in 1:n) {
        if (y[i] < 0) {
            z[i] = i
            dryr = dryr + 1
            ysaved[dryr] = abs(y[i])
        }
    }
    runlen = 0
    store = rep(NA, n)
    z1 <- z2 <- rep(NA, n)
    k = 0
    for (i in 1:n) {
        if (z[i] == i) {
            runlen = runlen + 1
            store[i] = runlen
        }
        else {
            runlen = 0
        }
        if ((runlen == 0) && (i != 1)) {
            #if ((runlen == 0 ) && (z[i-1] != 0)) {
            if (z[i - 1] != 0) {
                k = k + 1
                #print(cbind((i-1),store[i-1]))
                z1[k] <- (i - 1)
                z2[k] <- store[i - 1]
            }
        }
    }  #i
    #zz1 <- subset(z1,z1) #end index of drought
    #zz2 <- subset(z2,z2) #spell length
    zz1 = na.omit(z1)
    zz2 = na.omit(z2)
    nspell = k  #number of drought spells
    #print(cbind(zz1,zz2))
    #calculate average drought intensity
    dintensity <- rep(NA, nspell)
    for (i in 1:nspell) {
        sumdef = 0
        endidx = zz1[i]
        begidx = endidx - zz2[i] + 1
        for (j in begidx:endidx) {
            sumdef = sumdef + abs(y[j])
        }
        dintensity[i] <- sumdef/zz2[i]
    }
    avgdri <- mean(dintensity)  #average drought run intensity
    #avgdri <- median(dintensity) #median drought run intensity
    maxrunlen <- max(zz2)  #maximum drought run length
    maxdeficit <- max(ysaved)  #maximum deficit
    avgdeficit <- sum(ysaved)/dryr  #Ave Deficit
    outlist = list(dintensity = dintensity, avgdri = avgdri, dryr = dryr, maxrunlen = maxrunlen, 
        maxdeficit = maxdeficit, avgdeficit = avgdeficit, sumy = sum(ysaved), leny = length(ysaved))
    #outlist=list('dryr'=dryr,'maxrunlen'=maxrunlen)
    return(outlist)
}