# R code accompanying lecture 5 for 36-350, Fall 2012
  # See lecture for context


  # Example of predicting from a model

# Predict response values from a power-law scaling model
# Presumes: power-law model objects are lists containing a scaling factor
  # called "y0" and a scaling exponent called "a"
# Inputs: fitted power-law model (object), vector of values at which to make
  # predictions at (newdata)
  # Argument names follow conventions for other predict() functions
# Outputs: vector of predicted response values
predict.plm <- function(object, newdata) {
  # Check that object has the right components
  stopifnot("a" %in% names(object), "y0" %in% names(object))
  a <- object$a
  y0 <- object$y0
  # Sanity check the inputs
  stopifnot(is.numeric(a),length(a)==1)
  stopifnot(is.numeric(y0),length(y0)==1)
  stopifnot(is.numeric(newdata))
  return(y0*newdata^a)  # Actual calculation and return
}

	# Examples of plotting

# Plot fitted curve from power law model over specified range
# Inputs: list containing parameters (plm), start and end of range (from, to)
# Outputs: TRUE, silently, if successful
# Side-effect: Makes the plot
  # Take sanity-checking of parameters as read
  y0 <- plm$y0 # Extract parameters
  a <- plm$a
  f <- function(x) { return(y0*x^a) }
  curve(f(x),from=from,to=to)
  # Return with no visible value on the terminal
  invisible(TRUE)
}


	# With other arguments for curve()
          # illustrates the use of "..." to pass along extra arguments

# Plot fitted curve from power law model over specified range
# Inputs: list containing parameters (plm), start and end of range (from, to),
  # other graphical arguments for curve (...)
# Outputs: TRUE, silently, if successful
# Side-effect: Makes the plot
plot.plm.2 <- function(plm,...) {
  y0 <- plm$y0
  a <- plm$a
  f <- function(x) { return(y0*x^a) }
  # from and to are possible arguments to curve()
  curve(f(x),...)
  invisible(TRUE)
}


	# Using predict.plm()

# Plot fitted curve from power law model over specified range
# Inputs: list containing parameters (plm), start and end of range (from, to),
  # number of points for plotting (n), other graphical arguments for plot (...)
# Outputs: TRUE, silently, if successful
# Side-effect: Makes the plot
plot.plm.3 <- function(plm,from,to,n=101,...) {
  x <- seq(from=from,to=to,length.out=n)
  y <- predict.plm(object=plm,newdata=x)
  plot(x,y,...)
  invisible(TRUE)
}





	# Examples of recursion

# Calculate a factorial
# Input: a number (n)
# Output: the factorial of n
# Presumes: n is a single positive integer
my.factorial <- function(n) {
  if (n == 1) {  # Base case
    return(1)
  } else {  # Recursion
    return(n*my.factorial(n-1))
  }
}

# Calculate a Fibonacci number
# Input: a number (n)
# Output: the nth Fibonacci number
# Presumes: n is a single positive integer
fib <- function(n) {
  if ( (n==1) || (n==0) ) {  # Base cases
   return(1)
  } else {  # Recursion
   return (fib(n-1) + fib(n-2))
  }
}
   # EXERCISE: Try tracing through fib(5)



	# Our resource allocation example, made more modular and recursive

# Adjust an output plan until it's under budget and close to full utilization
# Inputs: vector of initial planned output levels (output), matrix giving
  # resources needed per unit of output (factory), vector of resources
  # available (available), vector of acceptable idle resource amounts
  # (slack), factor by which to adjust plans (tweak)
# Output: List containing acceptable output levels (output) and resources
  # required (needed)
planner <- function(output,factory,available,slack,tweak=0.1) {
  needed <- plan.needs(output,factory)
  if (all(needed <= available) && all(available-needed <= slack)) {
    return(list(output=output,needed=needed))
  }
  else {
    output <- adjust.plan(output,needed,available,tweak)
    return(planner(output,factory,available,slack))
  }
}


# Calculate resources needed by a given output plan
# Inputs: vector of output levels (output), matrix giving resources needed per
  # unit of output (factory)
# Output: column matrix giving resources needed
plan.needs <- function(output,factory) { factory %*% output }


# Make multiplicative adjustments to an unsatisfactory output plan
# Inputs: vector of planned output levels (output), resources needed by plan
  # (needed), vector of resources available (available), factor by which to
  # adjust plans
# Output: vector of new output levels
# Presumes: resources needed do not exactly match resources available
adjust.plan <- function(output,needed,available,tweak) {
  # If we use too much, reduce all planned output levels by the same factor
  if (all(needed >= available)) { return(output*(1-tweak)) }
  # If we use too little, increase all planned output levels
  if (all(needed < available)) { return((1+tweak)) }
  # If we use too much of some things and too little of others, randomly
    # change all output levels
    # See help(runif)
  return(output*runif(n=length(output),min=1-tweak,max=1+tweak))
  # N.B., should never call this when needed == available
}