#' @title Find the intercept points
#' @name intercept_points
#' @usage intercept_points(xs, hs, h_primes)
#'
#' @param xs sorted x values
#' @param hs sorted hx values
#' @param h_primes sorted hx/dx values
#'
#' @return Returns the intercepts of the envelope function
intercept_points <- function(xs, hs, h_primes){
  intercepts = hs - h_primes*xs
  return(intercepts)
}




#' @title Find the intersection points
#' @name intersection_points
#' @usage intersection_points(xs, hs, h_primes, bounds)
#'
#' @param xs sorted x values
#' @param hs sorted hx values
#' @param h_primes sorted hx/dx values
#' @param bounds vector of length 2 of the upper and lower bound
#'
#' @return Returns the intersection value on the x axis of the envelope lines
intersection_points <- function(xs, hs, h_primes, bounds) {

  L <- length(h_primes)
  b <- intercept_points(xs, hs, h_primes)
  # y-intercepts b of envelope function line segments,
  # for j = 1,..,k-1, the tangents at x_j and x_j+1 intersect at z
  z <- (b[2:L] - b[1:L-1]) / (h_primes[1:L-1] - h_primes[2:L])
  z <- c(bounds[1], z, bounds[2])

  return(z)
}

#' @title Find u_k
#' @name u_k
#' @usage u_k(x, xs, hs, h_primes, bounds)
#'
#' @param x values that were sampled
#' @param xs sorted x values
#' @param hs sorted hx values
#' @param h_primes sorted hx/dx values
#' @param bounds vector of length 2 of the upper and lower bound
#'
#' @return Returns the likelihood of u_k, the upper envelope
u_k <- function(x, xs, hs, h_primes, bounds){
  zs = intersection_points(xs, hs, h_primes, bounds)
  j = which.max(zs > x)-1

  return(hs[j] + (x-xs[j])*h_primes[j])
}






#' @title Find l_k
#' @name l_k
#' @usage l_k(x, xs, hs)
#'
#' @param x values that were sampled
#' @param xs sorted x values
#' @param hs sorted hx values
#'
#' @return Returns the likelihood of l_k, the squeezing function
l_k <- function(x, xs, hs){
  if(x < xs[1] || x > xs[length(xs)]){

    return(0)
  }else{
    j = which.max(xs > x)-1
    return( ((xs[j+1]-x)*hs[j] + (x-xs[j])*hs[j+1])/(xs[j+1] - xs[j]) )

  }

}


#' @title Find h
#' @name h
#' @usage h(x)
#'
#' @param x values that were sampled
#'
#' @return returns the log distribution
h <- function(x) {
  #return (log(f(x, ...)))
  return(log(f(x, ...)))
}



#' @title Find h_prime
#' @name h_prime
#' @usage h_prime(x)
#'
#' @param x values that were sampled
#'
#' @return returns derivative of the log of the distribution
h_prime <- function(x) {
  return (numDeriv::grad(h, x))
}




#' @title Find abcissa
#' @name set_abcissa
#' @usage set_abcissa(f, h, h_prime, bounds, ...)
#'
#' @param f values that were sampled
#' @param h sorted x values
#' @param h_prime sorted hx/dx values
#' @param bounds vector of length 2 of the upper and lower bound
#' @param ... further arguments passed to or from other methods
#'
#' @return generates the sorted x, sorted hx, and sorted hx/dx values
set_abcissa <- function(f, h, h_prime, bounds, ...) {

  # variables
  dist_from_bound = 0.1
  time_limit <- 10

  # the optim function errored or took hours to run if the bounds were too big
  # so setting a time limit of 10 seconds if optim cant find the max
  setTimeLimit(cpu = time_limit, elapsed = time_limit, transient = TRUE)
  on.exit({
    setTimeLimit(cpu = Inf, elapsed = Inf, transient = FALSE)
  })

  # setting initial x value
  if(bounds[1] == -Inf & bounds[2] == Inf){
    x_vals = -1
  }
  else if (bounds[1] == -Inf){
    x_vals = bounds[2]-dist_from_bound
  }
  else if (bounds[2] == Inf){
    x_vals = bounds[1]+dist_from_bound
  }
  else{
    x_vals <- tryCatch({
      optimize(f, interval=bounds, maximum=TRUE, ...)$maximum

    }, error = function(e) {
      print("Time limit reached, probably for optimize to search in Inf bounds")
      print(e$message)
      return((bounds[2]-bounds[1])/2)
    })
  }


  # create initial h and hprime values
  h_values <- h(x_vals)
  h_prime_vals <- h_prime(x_vals)

  # errored when 0 slope
  if(h_prime_vals[1] == 0){
    h_prime_vals[1] = 1e-7
  }

  # - value = downward slope
  # + value = upward slope
  # 0 = und slope
  slope_dir <- sign(h_prime_vals[1])
  enough_values <- FALSE
  conditions_unsatisfied <- TRUE

  # while loop to keep searching for points until we have enough + and - slopes
  while(conditions_unsatisfied) {

    L <- length(h_prime_vals)

    # if the slope is - then we stop looking for points on the right side
    # instead look for points on the left side where the slope is +
    # to make the points we add poinst to the end if the starting slope is -
    # and add points to the front if starting slope is +
    insert_idx <- ifelse(slope_dir < 0, 1, length(x_vals) + 1)

    # add a slope_dir value to the previous x value
    x_new <- x_vals[ifelse(slope_dir < 0, 1, L)] + slope_dir

    # if x_new > upper or less than lower bounds
    # bound break loop
    if (x_new > bounds[2] || x_new < bounds[1]) {
      break
    }


    # new h and h prime to add
    h_new <- h(x_new)
    h_prime_new <- h_prime(x_new)

    if (h_prime_new == 0 || is.null(h_prime_new)) {
      slope_dir <- slope_dir * 2
      next
    }


    # appending values
    x_vals <- append(x_vals, x_new, insert_idx)
    h_values <- append(h_values, h_new, insert_idx)
    h_prime_vals <- append(h_prime_vals, h_prime_new, insert_idx)

    if (bounds[2] - bounds[1] < 25) {
      # when bounds are really small, it increments too fast
      # so increment by + 2 instead of * 2
      slope_dir <- 2
    } else {
      slope_dir <- slope_dir * 2
    }

    opposite_vals_added <- ifelse(slope_dir > 0, sum(h_prime_vals < 0), sum(h_prime_vals > 0))
    enough_values <- opposite_vals_added >= 3
    conditions_unsatisfied <- !enough_values
  }

  return(list(xs=x_vals,
              hs=h_values,
              h_primes=h_prime_vals))
}



#' @title Find sample x's
#' @name sample_step
#' @usage sample_step(xs, hs, h_primes, bounds)
#'
#' @param xs sorted x values
#' @param hs sorted hx values
#' @param h_primes sorted hx/dx values
#' @param bounds vector of length 2 of the upper and lower bound
#'
#' @return samples x values of the distribution
sample_step <- function(xs, hs, h_primes, bounds){

  probs = c()

  zs = intersection_points(xs, hs, h_primes, bounds)
  ics = intercept_points(xs, hs, h_primes)

  for(j in 1:length(zs)-1){
    #prob = (exp(h_primes[j]*zs[j+1])-exp(h_primes[j]*zs[j]))*exp(ics[j])
    #prob = (h_primes[j]*(.5*zs[j+1]^2)+ics[j]*zs[j+1]) - (h_primes[j]*(.5*zs[j]^2)+ics[j]*zs[j])
    # if(h_primes[j] == 0){
    #   prob = (zs[j+1]-zs[j])*exp(ics[j])
    # }else{
    #   prob = (exp(h_primes[j]*zs[j+1]+ics[j])-exp(h_primes[j]*zs[j]+ics[j]))/h_primes[j]
    # }
    prob = (exp(h_primes[j]*zs[j+1]+ics[j])-exp(h_primes[j]*zs[j]+ics[j]))/h_primes[j]
    #prob = (exp(h_primes[j] * x - h_primes[j] * zs[j+1]) * h_primes[j]) / (1 - exp(h_primes[j] * zs[j] - h_primes[j] * zs[j+1]))
    probs = append(probs, prob)
  }

  probs = probs/sum(probs)


  inds = 1:length(probs)
  #chunk index
  section_ind = sample(inds, size = 1, prob = probs)
  #sample index from probs to determine section

  # Sample u = Uniform(0, 1)
  u = runif(1)

  # Invert i^th piecewise exponential CDF to get a sample from that interval
  if(h_primes[section_ind] == 0){
    return(u * (zs[section_ind+1] - zs[section_ind]) + zs[section_ind])
  }
  else{
    x = log(u * exp(h_primes[section_ind] * zs[section_ind+1]) + (1 - u) * exp(h_primes[section_ind] * zs[section_ind]))
    x = x / h_primes[section_ind]

    return(x)
  }
}