D <- function(f,delta=.000001){ function(x){ (f(x+delta) - f(x-delta))/(2*delta)} }


# spring length is 1

spring = function(x,k=2,g=9.8){-g*x + (k/2)*(x-1)^2}

# two springs in series
springs = function(x, k=2, g=9.8){ sum(cumsum(-g*x) + (k/2)*(x-1)^2)}

evalAtGrid = function(x,y,fun){ 
   res = matrix(0,length(x), length(y));
   for (k in 1:length(x)){
      for (j in 1:length(y)) {
         res[k,j] = fun(c(x[k],y[j]));
      }
   }
   res
}
 
grad = function(fun){ function(x) {
  res = matrix(0,length(x),1);
  for (k in 1:length(x)) {
        offset = 0*x;
        offset[k] =  .000001;
	res[k] = (fun(x+offset) - fun(x-offset))/(2*.000001);
  }
  res
}
}

plotarrow = function(gradfun,x) {
  foo = gradfun(x);
  arrows(x[1], x[2], x[1]+foo[1], x[2]+ foo[2])
}

census = data.frame(list(year=seq(1790,1990,by=10),
pop=c( 3.929, 5.308, 7.240, 9.638, 12.866, 17.069,
 23.192, 31.443, 38.558, 50.156, 62.948, 75.996,
 91.972, 105.711, 122.775, 131.669, 150.697, 
179.323, 203.185, 226.546, 248.710)))
  
censusfit = function(A,Y,m){ sum((census$pop - A/(1+exp(-(m*(census$year-Y)))))^2)}
censuserror = function(x){ censusfit(x[1], x[2], x[3])}

# returns a function that will take a vector c(A,Y,m) and return the
# sum of squared residuals.
logisticerror = function(t,pop){ function(x){
	sum((pop - x[1]/(1+exp(-x[3]*(t-x[2]))))^2)}}

# for fitting a simple line y = mx + b; b is the first parameter, m is the second.
linearErrorFun = function(x,y){ function(p){
	sum( (y - p[1] - p[2]*x)^2 ) }}


# for fitting a function of the form z = f(x,y) up to second order polynomial.
# This could be an in-class demo.
quaderror = function(x,y,z) { function(p){
	sum((z - p[1] -p[2]*x - p[3]*y - p[4]*x*y - p[5]*x*x - p[6]*y*y)^2)}}