# Rather trivial example

# Make circumference-calculators for non-standard values of pi
# Inputs: effective value of pi (ratio.to.diameter)
# Outputs: Function which calculates circumference from diameter
make.noneuclidean <- function(ratio.to.diameter=pi) {
  circumference <- function(d) { return(ratio.to.diameter*d) }
  return(circumference)
}

# Will produce a ``not found'' error
circumference(10)
kings.i <- make.noneuclidean(3)
kings.i(10)
formals(kings.i)
body(kings.i)
environment(kings.i)
circumference(10) # Still not found

	# Less trivial example

# Make a linear predictor based on sample data
# Inputs: vectors for independent and dependent variable (x and y)
# Output: function which calculates expected response at new values of x
make.linear.predictor <- function(x,y) {
  linear.fit <- lm(y~x)
  predictor <- function(x) {
   return(predict(object=linear.fit,newdata=data.frame(x=x)))
  }
  return(predictor)
}

independent.variable <- runif(10)
dependent.variable <- 7 + 3*independent.variable   # No noise
my.predictor <- make.linear.predictor(x=independent.variable,
  y=dependent.variable)
all.equal(my.predictor(independent.variable),dependent.variable)
   # Wouldn't work with noise in dependent variable
rm(independent.variable,dependent.variable)
independent.variable  # Gives not-found error
dependent.variable    # Ditto
my.predictor(5)       # Should be 5*3+7=22
my.predictor(runif(10))



	# Even less trivial example: build the gradient operator, instead
	# of just a function that evaluates the gradient at a point

# Gradient operator
# Inputs: function (f), optional extras (...)
# Outputs: the function which calculates the gradient of f
# Presumes: f takes a vector argument and is numeric-valued
  # our gradient() function from last lecture exists
nabla <- function(f,...) {
  g <- function(x,...) { gradient(f=f,x=x,...) }
  return(g)
}

# Presume that our mse.for.optimization()  function from last time is around
mse.gradient <- nabla(mse.for.optimization)
# Next two should match
mse.gradient(c(6611,0.15),deriv.steps=c(1,1e-6))
gradient(mse.for.optimization,c(6611,0.15),c(1,1e-6))
# Check that ellipses are working properly by changing a default value
gradient(mse.for.optimization,c(6611,0.15),c(1,1e-6),Y=2*gmp$pcgmp)
mse.gradient(c(6611,0.15),deriv.steps=c(1,1e-6),Y=2*gmp$pcgmp)

	# As mentioned last time, the first-difference method for approximating
	# partial derivatives is not very good, the numDeriv package has
	# better approximation algorithms in its grad() function

# Gradient operator
# Inputs: function (f), optional extras (...)
# Outputs: the function which calculates the gradient of f
# Presumes: f takes a vector argument and is numeric-valued
  # the numDeriv package is on the system
del <- function(f,...) {
  # If the numDeriv library isn't loaded, load it or die trying
  require(numDeriv)  
  # Invoke its grad() function in our function definition
  g <- function(x,...) { grad(func=f,x=x, ...)}
  return(g)
}



	# Long example: write surface() as an analog to curve(), using
	# the standard graphics function contour() to do actual plots

	# First attempt
# Plot contours of a two-argument function
# Inputs: function (f), limits for x and y (from.x,to.x,from.y,to.y), number
  # of values to space between limits (n.x,n.y), optional extras for contour()
  # via ...
# Outputs: Invisibly, a list with the x and y coordinates, and the matrix of
  # z coordinates
# Presumes: f takes a vector (of length 2!) and returns a single number
surface.0 <- function(f,from.x=0,to.x=1,from.y=0,to.y=1,n.x=101,
  n.y=101,...) {
  # Make sequences with the right limits and sizes
  x.seq <- seq(from=from.x,to=to.x,length.out=n.x)
  y.seq <- seq(from=from.y,to=to.y,length.out=n.y)
  # Make an (n.x*n.y) by 2 data frame with all combinations of x and y values
  plot.grid <- expand.grid(x=x.seq,y=y.seq)
  # Apply f to reach row of the plotting-grid dataframe
  z.values <- apply(plot.grid,1,f)
  # Reshape into an n.x by n.y matrix
  z.matrix <- matrix(z.values,nrow=n.x)
  # Make the contour plot
  contour(x=x.seq,y=y.seq,z=z.matrix,...)
  # Invisibly, return the values we calculated
  invisible(list(x=x.seq,y=y.seq,z=z.matrix))
}

# First graphic
	# Note use of anonymous function
surface.0(function(p){return(sum(p^3))},from.x=-1,from.y=-1)


	# Second attempt at surface()
	# Use substitute() and eval() to hold off on evaluating the
	# expression until the right values are assembled
	# idea taken from source to curve()
# Plot contours of a two-argument function
# Inputs: R expression (expr), limits for x and y (from.x,to.x,from.y,to.y),
  # number of values to space between limits (n.x,n.y), optional extras for
  # contour() via ...
# Outputs: Invisibly, a list with the x and y coordinates, and the matrix of
  # z coordinates
# Presumes: expr is a valid R expression involving x and y, any other names in
  # it can be found in the environment the function is called from
surface.1 <- function(expr,from.x=0,to.x=1,from.y=0,to.y=1,n.x=101,
  n.y=101,...) {
  # Make sequences with the right limits and sizes
  x.seq <- seq(from=from.x,to=to.x,length.out=n.x)
  y.seq <- seq(from=from.y,to=to.y,length.out=n.y)
  # Make an (n.x*n.y) by 2 data frame with all combinations of x and y values
  plot.grid <- expand.grid(x=x.seq,y=y.seq)
  # Recover the unevaluated, but evaluable, R expression from the first
    # argument
  unevaluated.expression <- substitute(expr)
  # ... and evaluate it in the environment set up by plot.grid
  z.values <- eval(unevaluated.expression,envir=plot.grid)
  # Reshape the vector of values into a matrix
  z.matrix <- matrix(z.values,nrow=n.x)
  # Make the contour plot
  contour(x=x.seq,y=y.seq,z=z.matrix,...)
   # Return everything we calculated
  invisible(list(x=x.seq,y=y.seq,z=z.matrix))
}

# Second figure
surface.1(abs(x^3)+abs(y^3),from.x=-1,from.y=-1)


	# Third attempt, using function creation and outer
# Plot contours of a two-argument function
# Inputs: R expression (expr), limits for x and y (from.x,to.x,from.y,to.y),
  # number of values to space between limits (n.x,n.y), optional extras for
  # contour() via ...
# Outputs: Invisibly, a list with the x and y coordinates, and the matrix of
  # z coordinates
# Presumes: expr is a valid R expression involving x and y, any other names in
  # it can be found in the environment the function is called from
surface.2 <- function(expr,from.x=0,to.x=1,from.y=0,to.y=1,n.x=101,
  n.y=101,...) {
  # Make appropriate coordinate sequences
  x.seq <- seq(from=from.x,to=to.x,length.out=n.x)
  y.seq <- seq(from=from.y,to=to.y,length.out=n.y)
  # Recover the actual R expression from the first argument
  unevaluated.expression <- substitute(expr)
  # Define a new function that evaluates that expression
  z <- function(x,y) {
    return(eval(unevaluated.expression,envir=list(x=x,y=y)))
  }
  # Run that function on all combinations of coordinates
  z.values <- outer(X=x.seq,Y=y.seq,FUN=z)
  # Reshape the vector of z.values into a matrix
  z.matrix <- matrix(z.values,nrow=n.x)
  # Make the contour plot
  contour(x=x.seq,y=y.seq,z=z.matrix,...)
  # Return all the numbers we calculated
    # How would you also return the function z as well?
  invisible(list(x=x.seq,y=y.seq,z=z.matrix))
}

# Third plot
surface.2(x^4-y^4,from.x=-1,from.y=-