#  File containing my commonly used R functions.
#
#  source("/home/waichler/R/functions1.R")

# Linear extrapolation.  Get y as a function of some x, given x1, x2, y1, y2.
# R's approx() does interpolation, but I couldn't find a built-in function to
# do the general case of finding a point along a line.
EXTRAPOLATE <- function (x, x1, x2, y1, y2) {
  m <- (y2-y1)/(x2-x1)  # slope
  b <- y2 - m*x2        # intercept
  y <- m*x + b
  return(y)
}

# This function is intended to help graphs more readable by replacing
# consecutive zero values in the vector with NA.  The beginning and end
# of a nonzero sequence are maintained at zero, however.  This allows one
# to show, for example, an ephemeral hydrograph without all the zero-valued
# line or symbols showing between flow events.
ELIMINATE.SEQUENTIAL.ZEROS <- function (vec) {
  zero.ind <- which(vec == 0.0)
  vec[zero.ind] <- NA
  for(i in 2:(length(vec)-1)) {
    if(is.na(vec[(i-1)]) == FALSE && vec[(i-1)] != 0.0 && is.na(vec[i]) == TRUE) {
      vec[i] = 0.0
    }
    if (is.na(vec[i]) == TRUE && is.na(vec[(i+1)]) == FALSE) {
      vec[i] = 0.0
    }
  }
  return(vec)
}  # end ELIMINATE.SEQUENTIAL.ZEROS()

# This function is intended to help graphs more readable by replacing
# all zero values in the vector with NA.  This is a suitable approach for
# precipitation or other variables represented as points or
# height bars (vertical lines).
ELIMINATE.ALL.ZEROS <- function (vec) {
  zero.ind <- which(vec == 0.0)
  vec[zero.ind] <- NA
  return(vec)
}  # end ELIMINATE.ALL.ZEROS()

# This function figures out the best range and tick mark placing for
# the y1 and y2 axes.  The y1 settings are determined first, then the
# y2 tickmarks and labels are determined; y2 tick marks are level with those
# for y1.  If you are going to plot multiple variables on the same axis,
# concatenate those vectors as y1 and/or y2 before calling this function.
# If you want to fix a given limit, use y*min and/or y*max; otherwise, set
# these to NA.  num.int = number of desired intervals.  There is some very
# simple logic to help pretty() deal with situations where the intervals are
# not whole numbers, but no doubt this won't always work right.
SET.YAXES <- function (y1, y1min, y1max, y2, y2min, y2max, num.int) {
  # y1 axis first
  if (is.na(y1min) == FALSE) {
    miny1 <- y1min
  } else {
    miny1 <- min(y1)
  }
  if (is.na(y1max) == FALSE) {
    maxy1 <- y1max
  } else {
    maxy1 <- max(y1)
  }
  
  if( (maxy1-miny1) <= 2) {
    cat(file="", "Warning, SET.YAXES() may not work correctly for small range in y1.  Check results.\n")
    miny1 <- 1000 * miny1
    maxy1 <- 1000 * maxy1
    y1ticks <- pretty(miny1:maxy1, n=num.int)
    y1ticks <- y1ticks/1000
    miny1 <- min(y1ticks)
    maxy1 <- max(y1ticks)
  } else {
    y1ticks <- pretty(miny1:maxy1, n=num.int)
    miny1 <- min(y1ticks)
    maxy1 <- max(y1ticks)
  }
  
  # y2 axis
  if (is.na(y2min) == FALSE) {
    miny2 <- y2min
  } else {
    miny2 <- min(y2)
  }
  if (is.na(y2max) == FALSE) {
    maxy2 <- y2max
  } else {
    maxy2 <- max(y2)
  }

  if( (maxy2-miny2) <= 2) {
    cat(file="", "Warning, SET.YAXES() may not work correctly for small range in y2.  Check results.\n")
    miny2 <- 1000 * miny2
    maxy2 <- 1000 * maxy2
    y2ticks <- pretty(miny2:maxy2, n=num.int)
    y2ticks <- y2ticks/1000
    miny2 <- min(y2ticks)
    maxy2 <- max(y2ticks)
  } else {
    y2ticks <- pretty(miny2:maxy2, n=num.int)
    miny2 <- min(y2ticks)
    maxy2 <- max(y2ticks)
  }
  
  scaledy2 <- (y2 - miny2) / (maxy2 - miny2) 
  scaledy2 <- (scaledy2 * (maxy1 - miny1)) + miny1
  
  scaledy2ticks <- (y2ticks - miny2) / (maxy2 - miny2) 
  scaledy2ticks <- (scaledy2ticks * (maxy1 - miny1)) + miny1
  
  return(miny1, maxy1, y1ticks, scaledy2, miny2, maxy2, y2ticks, scaledy2ticks)
}  # end SET.YAXES()

# Function to convert input in y2 scale to y1 scale.  Useful if you want to
# add a horizontal line to plot, etc.
CONVERT.Y2 <- function (miny1, maxy1, miny2, maxy2, y2.coord) {
  y1.coord <- (y2.coord - miny2) / (maxy2 - miny2) 
  y1.coord <- (y1.coord * (maxy1 - miny1)) + miny1
  return(y1.coord)
}

# water year
WY <- function (month, calyear) {return (ifelse (month >= 10, calyear+1, calyear)) }

# Function to aggregate irregular timespans in the vector "vec".
# ind are the breakpoints for summing.
AGGREGATE.INFLOW <- function (vec, ind) {
  agg.vec <- NULL
  for(i in 2:length(ind)) {
    if(i < length(ind)) {
      agg.vec <- c(agg.vec, sum(vec[ ind[(i-1)]:(ind[i]-1) ]))
    } else {
      agg.vec <- c(agg.vec, sum(vec[ ind[(i-1)]:ind[i] ]))
    }
  }
  return(agg.vec)
}

# Function to compute annual average of a vector, weighted by the timespan
# in days that corresponds to each value.
ANNUAL.AVERAGE <- function (timespan, values) {
  total.days <- sum(timespan)
  weights <- timespan/total.days
  weighted.average <- sum(weights * values)
  return(weighted.average)
}

# Function to compute annual average of a vector, weighted by the timespan
# in days that corresponds to each value.
ANNUAL.AVERAGE <- function (timespan, values) {
  total.days <- sum(timespan)
  weights <- timespan/total.days
  weighted.average <- sum(weights * values)
  return(weighted.average)
}

# Partition water infiltrating into ground into quickflow, shallow gw,
# and deep gw.  The fraction for each input infil.water value is returned.
# infil.water = water entering ground (vector)
# lower = fraction entering designated pool at the lowest value of infil.water
# middle = fraction entering designated pool at the median value of infil.water
# upper = fraction entering designated pool at the highest value of infil.water
FLOW.FRACTION <- function (infil.water, lower, middle, upper) {
  lower <- lower/100
  middle <- middle/100
  upper <- upper/100
  min.infil.water <- min(infil.water)
  max.infil.water <- max(infil.water)
  med.infil.water <- median(infil.water)
  fractions <- ifelse(infil.water <= med.infil.water,
       (lower + (middle - lower) * (infil.water - min.infil.water)/(med.infil.water - min.infil.water)),
       (middle + (upper - middle) * (infil.water - med.infil.water)/(max.infil.water - med.infil.water))
                     )
  return(fractions)
}


# Get the midpoints of a set of classes, defined by minimum, maximum (extremes),
# and number of classes.
GET.CLASS.VALUES <- function (num.classes, min.value, max.value) {
  midpoints <- NULL
  increment <- (max.value - min.value) / num.classes
  value <- min.value
  while (value < max.value) {
    midpoints <- c(midpoints, mean(c(value, (value + increment))))
    value <- value + increment
  }
  return(midpoints)  
}  

# Get the set of regularly spaced values for a specified r=range and n=number of points.
# If n=1, then minimum is returned.  If n=2, then min and max are returned.
# If n>2, then min, max, and (n-2) points that are spaced (max-min)/(n-1) apart are returned.
GET.REG.SPACED.VALUES <- function (lower, upper, num.points) {
  if(num.points == 1) {
    return(lower)
  } else if (num.points == 2) {
    return(c(lower, upper))
  } else if (num.points > 2) {
    vec <- lower
    value <- lower
    increment <- (upper - lower) / (num.points - 1)
    while (value + increment < upper) {
      value <- value + increment
      vec <- c(vec, value)
#      value <- value + increment
    }
    vec <- c(vec, upper)
    return(vec)
  } else {
    cat("ERROR--function GET.REG.SPACED.VALUES() doesn't work with input given.\n")
  }
}  

# Return a vector of indices designating the min, max, and regularly-spaced intermediate points
# between the minimum and maximum of the supplied vector.
# If the index whose value is closest to the target is already selected, then
# the index of the next closest value is selected.
GET.INTERVAL.IND <- function (num.intermediate.points, vec) {
  ind <- which.min(vec)
  max.value <- max(vec)
  min.value <- min(vec)
  increment <- (max.value - min.value) / (num.intermediate.points + 1)
  value <- min.value
  for (i in 1:num.intermediate.points) {
    diff <- abs( vec - (value + increment) )
    order.diff <- order(diff)
    for (j in order.diff) {
      if( !is.element(j, c(ind, which.max(vec))) ) {
        ind <- c(ind, j)
        break
      }
    }
    value <- value + increment
  }
  ind <- c(ind, which.max(vec))
  return(ind)  
}  


# Function to get moving window of a vector, such as moving average or moving sum.
# Type of aggregation can be "moving.window" or "discrete".
# Input:  vector, number of timesteps framing a window (should be an ODD number if
# window is centered type), aggregating function, and type ("backward" or "centered".
# Output: vector of same length.
WINDOW.AGGREGATE <- function (vec, type.aggregation, numw, type, aggfun) {
  numsteps <- length(vec)

  if(type.aggregation == "moving.window") {
    if(type == "centered") {  # centered window
      numw2 <- floor(numw/2) # number of timesteps on either side of center step
      # index noting end of partial windows (near start of vector)
      ind.end.partial <- numw2
      # index noting beginning of partial windows (near end of vector)
      ind.beg.partial <- numsteps - floor(numw/2)
      boxcar <- NULL

      # First and last numw/2 timesteps are partial windows
      for (i in 1:ind.end.partial) { # first steps
        cmd <- paste(sep="", aggfun, "(vec[1:(i+numw2)])")
        boxcar[i] <- eval(parse(text=cmd))
      }
      for (i in ind.beg.partial:numsteps) {  # last steps
        cmd <- paste(sep="", aggfun, "(vec[(i-numw2):numsteps])")
        boxcar[i] <- eval(parse(text=cmd))
      }

      # All other timesteps
      for(i in (ind.end.partial+1):(ind.beg.partial-1)) {
        cmd <- paste(sep="", aggfun, "(vec[(i-numw2):(i+numw2)])")
        boxcar[i] <- eval(parse(text=cmd))      
      }
      # end centered window
    } else if (type == "backward") {
      # backward-looking window
      # index noting end of partial windows (near start of vector)
      ind.end.partial <- numw - 1
      boxcar <- NULL

      # First (numw-1) timesteps are partial windows
      for (i in 1:ind.end.partial) { # first steps
        cmd <- paste(sep="", aggfun, "(vec[1:i])")
        boxcar[i] <- eval(parse(text=cmd))
      }

      # All other timesteps
      for(i in (ind.end.partial+1):numsteps) {
        cmd <- paste(sep="", aggfun, "(vec[(i-numw+1):i])")
        boxcar[i] <- eval(parse(text=cmd))      
      }
      # end backward-looking window
    }  
    return(boxcar)
    # end moving window type of aggregation  
  } else if (type.aggregation == "discrete") {  
    # discrete aggregation--no timestep used more than once
    boxcar <- aggregate.ts(vec, ndeltat=numw, FUN=aggfun)
    return(boxcar)
  }
}


# Function to get moving average of vector1, using a window based on vector 2.
# Input:  vector1 (the target data to average); vector 2 (the units with which you
# make the window, typically Julian days; and half-width, the
# number of units on either side of the framing a window.
# Output: vector of same length.
MOVING.AVERAGE <- function (vec1, vec2, half.win.width) {
  numpoints <- length(vec1)

  boxcar <- rep(NA, numpoints)

  for(i in 1:numpoints) {
    ind <- which(vec2 >= (vec2[i] - half.win.width) & vec2 <= (vec2[i] + half.win.width))
    boxcar[i] <- mean(vec1[ind])
  }
  return(boxcar)
}


# Function to check for clustering of values within some width.
# Input:  vector (the units with which you make the window, typically Julian days;
# and half-width, the number of units on either side to frame the window.
# The input vector must be chronologically ordered!
# Output: vector of same length.
CHECK.TIME.CLUSTERING <- function (vec, half.win.width) {
  numpoints <- length(vec)

  flag <- rep("", numpoints)

  for(i in 1:numpoints) {
    ind <- which(vec >= (vec[i] - half.win.width) & vec <= (vec[i] + half.win.width))
    if(length(ind) > 1 ) flag[i] <- "TC"
  }
  return(flag)
}
