# gradient descent: find local minimum of function of one variable from initial guess
# numerically estimate the derivative
@everywhere function gradient_descent(
  f::Function,
  x0::Float64,
  delta::Float64, # shift is delta*df
  dx::Float64, # finite difference for numerical derivative evaluation
  maxiter::Int64 # interrupt and fail after maxiter steps
)
  counter=0

  # init
  x=x0

  while counter<maxiter
    # value at x and around
    val=f(x)
    valm=f(x-dx)
    valp=f(x+dx)
    # quit if minimum
    if(val<valm && val<valp)
      return(x,val)
    end

    # derivative
    df=(valp-val)/dx
    # step
    x=x-delta*df
    counter+=1
  end

  # fail
  return(Inf,Inf)
end

# Newton algorithm to compute extrema
# numerically estimate the derivatives
@everywhere function newton_maximum(
  f::Function,
  x0::Float64,
  dx::Float64, # finite difference for numerical derivative evaluation
  maxiter::Int64,
  tolerance::Float64,
  maxstep::Float64 # maximal size of step
)
  counter=0

  # init
  x=x0

  while counter<maxiter
    # value at x and around
    val=f(x)
    valm=f(x-dx)
    valp=f(x+dx)

    # derivative
    dfp=(valp-val)/dx
    dfm=(val-valm)/dx
    # second derivative
    ddf=(dfp-dfm)/dx

    #@printf(stderr,"% .15e % .15e % .15e % .15e\n",x,val,dfp,ddf)

    if abs(dfp/ddf)<tolerance
      # check it is a local maximum
      if ddf<0
	return(x,val)
      else
	return(Inf,Inf)
      end
    end

    # step
    step=dfp/abs(ddf)

    # step too large
    if abs(step)>maxstep
      step=maxstep*sign(step)
    end

    x=x+step

    # fail if off to infinity
    if x==Inf || x==-Inf
      return(x,val)
    end
    counter+=1
  end

  # fail
  return(Inf,Inf)
end