oct files for BFGS, numeric derivatives

// ========================= FiniteDifference ==========================
//  finite differences for numeric differentiation
//  the formulae here are based on those in Ox 3.20 (Doornik)
//  which references Rice
#include <oct.h>
DEFUN_DLD(FiniteDifference, args, ,"FiniteDifference, C++ version\n\
differences for NumGradient and NumHessian")
{
    double x = args(0).double_value();
    int order = args(1).int_value();
	int test;
	double eps, SQRT_EPS, DIFF_EPS, DIFF_EPS1, DIFF_EPS2, diff, d;

	eps = 2.2204e-16; // eps is machine precision - this is for i686
	SQRT_EPS = sqrt(eps); 
	DIFF_EPS = exp(log(eps)/2);
	DIFF_EPS1 = exp(log(eps)/3);
	DIFF_EPS2 = exp(log(eps)/4);
	if (order == 0) diff = DIFF_EPS;
	else if (order == 1) diff = DIFF_EPS1;
	else diff = DIFF_EPS2;

	test = (fabs(x) + SQRT_EPS) * SQRT_EPS > diff;
	
	if (test)
	{
		d = (fabs(x) + SQRT_EPS) * SQRT_EPS;
	}
	else
	{
		d = diff;
	}		

	return octave_value(d);
}

// ============================ BFGSMin.m =====================================
// 
//  Minimize a function of the form 
// 
//    [obj_value, ...] = f(arg, otherargs)
// 
//    or
// 
//    [obj_value, ...] = f(arg)
// 
//    with respect to arg. The first return of f() must be the value of f().
//    f() may have other returns, but they are not used here.
// 
//  Calling syntax:
// 
//    [theta, obj_value, iters, convergence] = BFGSMin(f, arg, otherargs, control)
//  or
//    [theta, obj_value, iters, convergence] = BFGSMin(f, arg, {1}, control)
// 
//    NOTE: the last argument (control) is optional
// 
// 
//  Input arguments:
// 
//    REQUIRED
//    f - (string) the function to minimize
//    arg - the arg we min with respect to, a vector
//    otherargs - a cell array of other arguments of func, or a placeholder, e.g., {1}
// 
//    OPTIONAL
//    control - 3x1 vector
// 	   * 1st elem. controls maximum iterations
// 		   scalar > 0 - max. iters
// 		   or -1 for infinity (default)
// 	   * 2nd elem. is verbosity control 
// 		   0 = no results printed (default)
// 		   1 = intermediate results
// 		   2 = only last iteration results
// 	   * 3rd elem. specifies convergence criterion
// 		   1 = function, gradient, and parameter change
// 			   are all tested (default)
// 		   0 = only function conv tested
// 
//  Outputs:
//    theta - the minimizer
//    obj_value - the minimized funtion value
//    iters - number of iterations used
//    convergence (1 = success, 0 = failure)

#include <oct.h>
#include <octave/parse.h>
#include <octave/Cell.h>
#include <octave/lo-mappers.h>

// define argument checks
static bool
any_bad_argument(const octave_value_list& args)
{

	if (args.length() > 4)
	{
		error("\nBFGSMin: you have provided more than 4 arguments\n");
		return true;
	}	
	if (args.length() < 3)
	{
		error("\nBFGSMin: you have provided fewer than 3 arguments\n");
		return true;
	}	

	if (!args(0).is_string())
    {
        error("BFGSMin: first argument must be string holding objective function name");
        return true;
    }
	if (!args(1).is_real_matrix())
    {
        error("BFGSMin: second argument must be vector of initial parameter values");
        return true;
    }
	if (!args(2).is_cell())
    {
        error("BFGSMin: third argument must a cell of other arguments\n\
		 (or a placeholder such as {1})");
        return true;
    }
	if (args.length() == 4)
	{
		if (!args(3).is_real_matrix())
    	{
        	error("BFGSMin: fourth argument, if supplied must a 3x1 vector of control options");
	        return true;
    	}
	}	
    return false;
}


DEFUN_DLD(BFGSMin, args, , "Usage:\n\
[m, o, i, c] = BFGSMin(f, x, otherargs, control)\n\
[m, o, i, c] = BFGSMin(f, x, {1}, control)\n\
\n\
BFGS minimization of f() with respect to x\n\
f() should be of the form\n\
[o, possibly more returns...] = f(x, otherargs)\n\
where otherargs is a cell array\n\
\n\
or if there are no other arguments\n\
[o, possibly more returns...] = f(x)\n\
\n\
Inputs:\n\
f: string - the function to be minimized\n\
x: vector - starting values of arg w.r.t. which minimization is done\n\
otherargs: cell array - additional arguments to function\n\
	(e.g., data)\n\
	if there are none, use a place holder such as {1}\n\
control: (optional) 3x1 vector of switches\n\
	control(1): max iters, use -1 for infinity\n\
	control(2): verbosity\n\
		0 (default): no screen output\n\
		1: report at every iteration\n\
		2: report at last iteration\n\
	control(3): strict or weak convergence criterion\n\
		1 (default): function, gradient and parameter\n\
			convergence required\n\
		any other value: only function convergence required\n\
Outputs:\n\
m: the minimizing value of x\n\
o: the value of f() at m\n\
i: number of iterations required\n\
c: convergence message\n\
	c = 1: normal convergence\n\
	c = 2: failure of algorithm\n\
	c = -1: max iters exceeded\n\
\n\
Example:\n\
\n\
octave:1> function obj_value = quadratic(x)\n\
> obj_value = x'*x;\n\
> endfunction\n\
octave:2> x = rand(5,1)\n\
x =\n\
\n\
  0.81472\n\
  0.13548\n\
  0.90579\n\
  0.83501\n\
  0.12699\n\
\n\
octave:3> BFGSMin(\"quadratic\", x, {1}, [10,2,1]);\n\
\n\
BFGSMin final results: Iteration 1\n\
Stepsize 0.0000000\n\
Objective function value     0.0000000000\n\
Function conv 1  Param conv 1  Gradient conv 1\n\
  params  gradient  change\n\
 -0.0000  -0.0000   0.0000\n\
  0.0000   0.0000  -0.0000\n\
 -0.0000  -0.0000   0.0000\n\
 -0.0000  -0.0000   0.0000\n\
 -0.0000  -0.0000  -0.0000\n\
")

{
	int nargin = args.length();
	if (nargin < 1)
    {
      	error("BFGSMin: you must supply at least 2 arguments");
	    return octave_value_list();
    }

    // check the arguments
	if (any_bad_argument(args)) return octave_value_list();


	std::string f (args(0).string_value());
    Matrix theta (args(1).matrix_value());
	Cell otherargs (args(2).cell_value());
	
	octave_value_list f_args(2,1);
	octave_value_list g_args = args;
	octave_value_list step_args(4,1);
	octave_value_list f_return;

	int criterion, max_iters, convergence, verbosity, iter;
	int gradient_failed, i, j, conv1, conv2, conv3;
	int k = theta.rows();

	double func_tol, param_tol, gradient_tol, stepsize, obj_value;
	double last_obj_value, denominator, test, tempscalar;
	Matrix control, thetain, thetanew, H, g, d, g_new, p, q, temp, H1, H2;
	
	f_args(0) = theta;
	f_args(1) = otherargs;

	step_args(0) = f; // the 4 elements are func,  direction, theta, otherargs
	step_args(3) = otherargs;
	
	// tolerances
	func_tol  = 10*sqrt(2.22e-16);
	param_tol = 1e-5;
	gradient_tol = 1e-4;  

	// Default values for controls
	max_iters = INT_MAX; // no limit on iterations
	verbosity = 0; // by default don't report results to screen
	criterion = 1; // strong convergence required

 	// use provided controls, if applicable
	if (args.length() == 4)
	{
		control = args(3).matrix_value();
		max_iters = (int) control(0);
		if (max_iters == -1) max_iters = INT_MAX;
		verbosity = (int) control(1);
		criterion = (int) control(2);
	}	
	
	// initialize things
	convergence = -1; // if this doesn't change, it means that maxiters were exceeded
 	thetain = theta;
	H = identity_matrix(k,k);

	// Initial obj_value
	f_return = feval(f, f_args); 
	last_obj_value = f_return(0).double_value();

	
	// initial gradient
	f_return = feval("NumGradient", args);
	g = f_return(0).matrix_value();
	g = g.transpose();

	// check that gradient is ok
	gradient_failed = 0;  // test = 1 means gradient failed
	for (i=0; i<k;i++)
	{
		gradient_failed = gradient_failed + xisnan(g(i));
	}
	if (gradient_failed)
	{
		error("BFGSMin: Initial gradient could not be calculated: exiting");
		convergence = 2;
	}

	// MAIN LOOP STARTS HERE
 	for (iter = 0; iter < max_iters; iter++)
	{
 		d = -H*g;

		// Regular step
		step_args(1) = d;
		step_args(2) = theta;
		f_return = feval("NewtonStep", step_args);
		stepsize = f_return(0).double_value();
		obj_value = f_return(1).double_value();
		
		// Steepest descent if regular step fails
		if (xisnan(stepsize)) 
		{
			d = -g; // try steepest descent
			step_args(1) = d;
			f_return = feval("NewtonStep", step_args);
			stepsize = f_return(0).double_value();
			obj_value = f_return(1).double_value();

			warning("BFGSMin: Stepsize failure in BFGS direction, trying steepest descent direction");
			// if that didn't work, were out of luck
			if (xisnan(stepsize))
			{
				warning("BFGSMin: unable to find direction of improvement: exiting");
				theta = thetain;
				convergence = 2;
				break;
			}
 		}
      	p = stepsize*d;

 		// check normal convergence: all 3 must be satisfied
		// function convergence
		if (fabs(last_obj_value) > 1.0)
		{
 			conv1 = (fabs(obj_value - last_obj_value)/fabs(last_obj_value)) < func_tol;
		}
		else
		{
			conv1 = fabs(obj_value - last_obj_value) < func_tol;
		}	
		// parameter change convergence
		temp = theta.transpose() * theta;
		test = sqrt(temp(0,0));
		if (test > 1)
		{
			temp = p.transpose() * p;
			conv2 = (temp(0,0) / test) < param_tol ;
			
		}
		else
		{
			temp = p.transpose() * p;
			conv2 = temp(0,0) < param_tol;
		}		
		// gradient convergence
		temp = g.transpose() * g;
		conv3 = sqrt(temp(0,0)) < gradient_tol;


		// Want intermediate results?
		if (verbosity == 1)
		{
			printf("\nBFGSMin intermediate results: Iteration %d\n",iter);
			printf("Stepsize %8.7f\n", stepsize);
			printf("Objective function value %16.10f\n", last_obj_value);
			printf("Function conv %d  Param conv %d  Gradient conv %d\n", conv1, conv2, conv3);	
			printf("  params  gradient  change\n");
	 		for (j = 0; j<k; j++)
			{
				printf("%8.4f %8.4f %8.4f\n",theta(j),g(j),p(j));
			}		
			printf("\n");
		}
	
		// Are we done?
		if (criterion == 1)
		{
			if (conv1 && conv2 && conv3)
			{
		  	convergence = 1;
	     	break;
			}
		}		
		else if (conv1)
		{
				convergence = 1;
				break;
		}
		
			
  	    last_obj_value = obj_value;	
		thetanew = theta + p;
		theta = thetanew;

		// get gradient at new parameter value
		g_args(1) = theta;
		f_return = feval("NumGradient", g_args);
		g_new = f_return(0).matrix_value();
		g_new = g_new.transpose();

		// Check that gradient is ok
		gradient_failed = 0;  // test = 1 means gradient failed
		for (i=0; i<k;i++)
		{
			gradient_failed = gradient_failed + xisnan(g_new(i));
		}

		// Hessian update if gradient ok		
		if (!gradient_failed)
		{
			q = g_new-g;
  	    	g = g_new;
			temp = q.transpose() * p;
			denominator = temp(0,0);
	  	  	if (denominator < 2.2e-16)  // reset Hessian if necessary
			{  	
				// if we already reset Hessian then we're out of luck
				if (H == identity_matrix(k,k))
				{					 
					convergence = 2;
					theta = thetain;

				  	warning("BFGSMin: unable to find direction of improvement: exiting");
					break;
				}	
				// Here's the Hessian re-set
				H = identity_matrix(k,k);
			}	
		  	else // normal update
			{
				temp = (1.0+(q.transpose() * H * q) / denominator) / denominator;
				tempscalar = temp(0,0);
				H1 = tempscalar * (p * p.transpose());
				H2 = (p * q.transpose() * H + H*q*p.transpose());
				H2 = H2 / denominator;
			  	H = H + H1 - H2;
  			}
		}
		else // failed gradient - try Hessian reset
		{
			// if we already tried it, we're out of luck
			if (H == identity_matrix(k,k))
			{
				convergence = 2;
				theta = thetain;
				warning("BFGSMin: failure of gradient: exiting");
				break;
			}	
			H = identity_matrix(k,k);
		}	
	}
	
// 	// Want last iteration results?
	if (verbosity == 2)
	{
		printf("\nBFGSMin final results: Iteration %d\n",iter);
		printf("Stepsize %8.7f\n", stepsize);
		printf("Objective function value %16.10f\n", last_obj_value);
		printf("Function conv %d  Param conv %d  Gradient conv %d\n", conv1, conv2, conv3);	
		printf("  params  gradient  change\n");
		for (j = 0; j<k; j++)
		{
			printf("%8.4f %8.4f %8.4f\n",theta(j),g(j),p(j));
		}		
		printf("\n");
	}
	f_return(0) = theta;
	f_return(1) = obj_value;
	f_return(2) = iter;
	f_return(3) = convergence;
	return octave_value_list(f_return);
}

// ============================ BisectionStep =====================================
// 
//  this is for use by BFGSMin 
//
// Uses bisection to find a stepsize that leads to a decrease, then
// continues until no further improvement
// Michael Creel [email protected]
//
// 
//  usage:
// 	   [a, obj_value] = BisectionStep(f, dx, x, otherargs)
//  inputs: 
//		f: the objective function
// 	   	dx: the direction
//		x: the parameter value
//		args: the other arguments of the function, in a cell array.

#include <oct.h>
#include <octave/parse.h>

DEFUN_DLD(BisectionStep, args, , "BisectionStep.cc")
{
    std::string func = args(0).string_value();
    Matrix dx = args(1).matrix_value();
    Matrix x = args(2).matrix_value();
	octave_value otherargs = args(3);
	octave_value_list f_args(2,1);

	f_args(1) = otherargs;
		
	double obj_0, obj, a;
	octave_value_list f_return;
	octave_value_list stepobj(2,1);

	f_args(0) = x;

	Matrix x_in = x;

	// possibly function returns a cell array
	// obj. value will be in first position
	f_return = feval(func, f_args); 
	obj_0 = f_return(0).double_value();

	a = 1.0;

	// this first loop goes until an improvement is found
  	while (a > 4e-16) // limit iterations
	{
		f_args(0) = x + a*dx;
		f_return = feval(func, f_args); 
		obj = f_return(0).double_value();
 
 		// reduce stepsize if worse, or if function can't be evaluated
		if ((obj > obj_0) || isnan(obj))
		{
			a = 0.5 * a;
		}	
		else
		{
			obj_0 = obj;
			break;
		}
	}
	
	// now keep going until we no longer improve, or reach max trials
	while (a > 4e-16)
	{
	   	a = 0.5*a; 
		f_args(0) = x + a*dx;
		f_return = feval(func, f_args); 
		obj = f_return(0).double_value();
 
		// if improved, record new best and try another step
		if ((obj < obj_0) & !isnan(obj))
		{
			obj_0 = obj;
		}	
		else
		{
			a = a / 0.5; // put it back to best found
			break;
		}				
	}

	stepobj(0) = a;
	stepobj(1) = obj_0;
	return octave_value_list(stepobj);
}

// ============================ NewtonStep =====================================
// 
//  this is for use by BFGSMin
// 
//  Uses a Newton interation to attempt to find a good stepsize
//  falls back to BisectionStep if no improvement found
//  Michael Creel [email protected]
//  13/01/2004
// 
//  usage:
// 	   [a, obj_value] = NewtonStep(f, dx, x, otherargs)
//  inputs: 
//		f: the objective function
// 	   	dx: the direction
//		x: the parameter value
//		args: the other arguments of the function, in a cell array.

#include <oct.h>
#include <octave/parse.h>

DEFUN_DLD(NewtonStep, args, , "NewtonStep.cc")
{
    std::string func = args(0).string_value();
    Matrix dx = args(1).matrix_value();
    Matrix x = args(2).matrix_value();
	octave_value otherargs = args(3);

	octave_value_list f_args(2,1);
	f_args(1) = otherargs;

	double obj, obj_0, obj_left, obj_right, delta, a, gradient, hessian;
	octave_value_list f_return;
	octave_value_list stepobj(2,1);

	f_args(0) = x;

	Matrix x_in = x;
	
	gradient = 1.0;
	
	// possibly function return cell array
	// obj. value will be in first position
	f_return = feval(func, f_args); 
	obj = f_return(0).double_value();

	obj_0 = obj;
	
	delta = 0.001; // experimentation show that this is a good choice
	
	Matrix x_right = x + delta*dx;
	Matrix x_left = x  - delta*dx;

	// possibly function return cell array
	// obj. value will be in first position
	f_args(0) = x_right;
	f_return = feval(func, f_args); 
	obj_right = f_return(0).double_value();
	
	f_args(0) = x_left;
	f_return = feval(func, f_args); 
	obj_left = f_return(0).double_value();
	

  	gradient = (obj_right - obj_left) / (2*delta);  // take central difference
  	hessian =  (obj_right - 2*obj + obj_left) / pow(delta, 2.0);	
	hessian = fabs(hessian); // ensures we're going in a decreasing direction
	if (hessian <= 2e-16) hessian = 1.0; // avoid div by zero

	a = - gradient / hessian;  // hessian inverse gradient: the Newton step
	a = (a < 5.0)*a + 5.0*(a>=5.0); // Let's avoid extreme steps that might cause crashes
	// check that this is improvement
	f_args(0) = x + a*dx;
	f_return = feval(func, f_args); 
	obj = f_return(0).double_value();
 
	// if not, fall back to bisection
	if ((obj > obj_0) || isnan(obj))
	{
		 f_return = feval("BisectionStep", args);
		 a = f_return(0).double_value();
		 obj = f_return(1).double_value();
	}
	
	stepobj(0) = a;
	stepobj(1) = obj;
	return octave_value_list(stepobj);
}

// ============================= NumGradient =============================
//  Central difference gradient of a function of the form
//  f(theta, otherargs), where theta is a vector and otherargs is a cell array.
// 
//  The gradient is respect to the FIRST element of args.
// 
//  * Allows diff. of vector-valued function
//  * Uses systematic finite difference
// 
//  See ExampleNumGradient.m

#include <oct.h>
#include <octave/parse.h>
#include <octave/Cell.h>
#include <octave/lo-mappers.h>


// define argument checks
static bool
any_bad_argument(const octave_value_list& args)
{


	if (!args(0).is_string())
    {
        error("NumGradient: first argument must be string holding objective function name");
        return true;
    }
	if (!args(1).is_real_matrix())
    {
        error("NumGradient: second argument must be vector of initial parameter values");
        return true;
    }
	if (!args(2).is_cell())
    {
        error("NumGradient: third argument must a cell of other arguments\n\
		 (or a placeholder such as {1})");
        return true;
    }
    return false;
}

DEFUN_DLD(NumGradient, args, ,
"Usage:\n\
gradient = NumGradient(f, x, otherargs)\n\
\n\
Numeric central difference gradient of a function f() with respect to x\n\
\n\
f() should be of the form\n\
[f_value, possibly more returns...] = f(x, otherargs)\n\
where otherargs is a cell array\n\
\n\
or if there are no other arguments\n\
[f_value, possibly more returns...] = f(x)\n\
\n\
f_value can be a scalar or column vector.\n\
\n\
gradient will be dimension nxk, where n=dim(f_value) and k=dim(x)\n\
\n\
Inputs:\n\
f: string - the function to be minimized\n\
x: vector - the value at which gradient is taken\n\
otherargs: cell array - additional arguments to function\n\
	(e.g., data)\n\
	if there are none, use a place holder such as {1}\n\
Outputs:\n\
gradient: nxk matrix\n\
\n\
Example:\n\
\n\
function a = myfunc(x)\n\
	a = x'*x;\n\
endfunction\n\
\n\
x = rand(3,1)\n\
x =\n\
\n\
  0.81472\n\
  0.13548\n\
  0.90579\n\
\n\
g = NumGradient(\"myfunc\", x, {1})\n\
g =\n\
\n\
	1.62945  0.27095  1.81158\n\
")
{
	int nargin = args.length();
	if (nargin < 1)
    {
      	error("NumGradient: you must supply 3 arguments\n\
		even if the 3rd is just a placeholder such as {1}");
	    return octave_value_list();
    }

	// check the arguments
	if (any_bad_argument(args)) return octave_value_list();

    std::string f = args(0).string_value();
    Matrix parameter = args(1).matrix_value();
	octave_value otherargs = args(2);
	octave_value_list f_args(2,1);
	octave_value_list fdiff_args(2,1);
	
	f_args(0) = parameter;
	f_args(1) = otherargs;
 
 	Matrix obj_value, obj_left, obj_right;
	octave_value_list f_return;
	double p, d, delta, delta_right, delta_left;
	int i, j;
	
	// possibly function return cell array
	// obj. value will be in first position
	f_return = feval(f, f_args); 
	obj_value = f_return(0).matrix_value();
	
	const int n = obj_value.rows();
    const int k = parameter.rows();

    Matrix derivative(n, k);
	Matrix columnj;

    for (j=0; j<k; j++) // get 1st derivative by central difference
	{
        p = parameter(j);

		fdiff_args(0) = p;
		fdiff_args(1) = 1;
		f_return = feval("FiniteDifference", fdiff_args);
		delta = f_return(0).double_value();

        // right side
		parameter(j) = d = p + delta;
		delta_right = d - p;
		f_args(0) = parameter;
   	  	f_return = feval(f, f_args);
		obj_right = f_return(0).matrix_value();

	  	// left size
		d = p - delta;
		parameter(j) = d;
		delta_left = p - d;
		f_args(0) = parameter;
	  	f_return = feval(f, f_args);
		obj_left = f_return(0).matrix_value();
		
		parameter(j) = p;  // restore original parameter 
		columnj = (obj_right - obj_left) / (delta_right + delta_left);
		for (i=0; i<n; i++)
		{
	        derivative(i, j) = columnj(i);
		}
			
	}
	return octave_value(derivative);

}

// ============================ NewtonMin.m =====================================
// 
//  Minimize a function of the form 
// 
//    [obj_value, ...] = f(arg, otherargs)
// 
//    or
// 
//    [obj_value, ...] = f(arg)
// 
//    with respect to arg. The first return of f() must be the value of f().
//    f() may have other returns, but they are not used here.
// 
//  Calling syntax:
// 
//    [theta, obj_value, iters, convergence] = NewtonMin(f, arg, otherargs, control)
//  or
//    [theta, obj_value, iters, convergence] = NewtonMin(f, arg, {1}, control)
// 
//    NOTE: the last argument (control) is optional
// 
// 
//  Input arguments:
// 
//    REQUIRED
//    f - (string) the function to minimize
//    arg - the arg we min with respect to, a vector
//    otherargs - a cell array of other arguments of func, or a placeholder, e.g., {1}
// 
//    OPTIONAL
//    control - 3x1 vector
// 	   * 1st elem. controls maximum iterations
// 		   scalar > 0 - max. iters
// 		   or -1 for infinity (default)
// 	   * 2nd elem. is verbosity control 
// 		   0 = no results printed (default)
// 		   1 = intermediate results
// 		   2 = only last iteration results
// 	   * 3rd elem. specifies convergence criterion
// 		   1 = function, gradient, and parameter change
// 			   are all tested (default)
// 		   0 = only function conv tested
// 
//  Outputs:
//    theta - the minimizer
//    obj_value - the minimized funtion value
//    iters - number of iterations used
//    convergence (1 = success, 0 = failure)

#include <oct.h>
#include <octave/parse.h>
#include <octave/Cell.h>
#include <octave/lo-mappers.h>

// define argument checks
static bool
any_bad_argument(const octave_value_list& args)
{

	if (args.length() > 4)
	{
		error("\nNewtonMin: you have provided more than 4 arguments\n");
		return true;
	}	
	if (args.length() < 3)
	{
		error("\nNewtonMin: you have provided fewer than 3 arguments\n");
		return true;
	}	

	if (!args(0).is_string())
    {
        error("NewtonMin: first argument must be string holding objective function name");
        return true;
    }
	if (!args(1).is_real_matrix())
    {
        error("NewtonMin: second argument must be vector of initial parameter values");
        return true;
    }
	if (!args(2).is_cell())
    {
        error("NewtonMin: third argument must a cell of other arguments\n\
		 (or a placeholder such as {1})");
        return true;
    }
	if (args.length() == 4)
	{
		if (!args(3).is_real_matrix())
    	{
        	error("NewtonMin: fourth argument, if supplied must a 3x1 vector of control options");
	        return true;
    	}
	}	
    return false;
}


DEFUN_DLD(NewtonMin, args, ,
"Usage:\n\
[m, o, i, c] = NewtonMin(f, x, otherargs, control)\n\
[m, o, i, c] = NewtonMin(f, x, {1}, control)\n\
\n\
Newton minimization of f() with respect to x\n\
f() should be of the form\n\
[o, possibly more returns...] = f(x, otherargs)\n\
where otherargs is a cell array\n\
\n\
or if there are no other arguments\n\
[o, possibly more returns...] = f(x)\n\
\n\
Inputs:\n\
f: string - the function to be minimized\n\
x: vector - starting values of arg w.r.t. which minimization is done\n\
otherargs: cell array - additional arguments to function\n\
	(e.g., data)\n\
	if there are none, use a place holder such as {1}\n\
control: (optional) 3x1 vector of switches\n\
	control(1): max iters, use -1 for infinity\n\
	control(2): verbosity\n\
		0 (default): no screen output\n\
		1: report at every iteration\n\
		2: report at last iteration\n\
	control(3): strict or weak convergence criterion\n\
		1 (default): function, gradient and parameter\n\
			convergence required\n\
		any other value: only function convergence required\n\
Outputs:\n\
m: the minimizing value of x\n\
o: the value of f() at m\n\
i: number of iterations required\n\
c: convergence message\n\
	c = 1: normal convergence\n\
	c = 2: failure of algorithm\n\
	c = -1: max iters exceeded\n\
\n\
Example:\n\
\n\
octave:1> function obj_value = quadratic(x)\n\
> obj_value = x'*x;\n\
> endfunction\n\
octave:2> x = rand(5,1)\n\
x =\n\
\n\
  0.81472\n\
  0.13548\n\
  0.90579\n\
  0.83501\n\
  0.12699\n\
\n\
octave:3> NewtonMin(\"quadratic\", x, {1}, [10,2,1]);\n\
\n\
NewtonMin final results: Iteration 1\n\
Stepsize 0.0000000\n\
Objective function value     0.0000000000\n\
Function conv 1  Param conv 1  Gradient conv 1\n\
  params  gradient  change\n\
 -0.0000  -0.0000   0.0000\n\
  0.0000   0.0000  -0.0000\n\
 -0.0000  -0.0000   0.0000\n\
 -0.0000  -0.0000   0.0000\n\
 -0.0000  -0.0000  -0.0000\n\
")
{
	int nargin = args.length();
	if (nargin < 1)
    {
      	error("NewtonMin: you must supply at least 2 arguments");
	    return octave_value_list();
    }

    // check the arguments
	if (any_bad_argument(args)) return octave_value_list();

	std::string f = args(0).string_value();
    Matrix theta = args(1).matrix_value();
	octave_value otherargs = args(2);
	
	octave_value_list f_args(2,1);
	octave_value_list g_args = args;
	octave_value_list step_args(4,1);
	octave_value_list f_return;

	int criterion, max_iters, convergence, verbosity, iter;
	int gradient_failed, i, j, conv1, conv2, conv3;
	int k = theta.rows();

	double func_tol, param_tol, gradient_tol, stepsize, obj_value;
	double last_obj_value, test, tempscalar;
	Matrix control, thetain, thetanew, H, g, d, g_new, p, q, temp;
	
	f_args(0) = theta;
	f_args(1) = otherargs;

	step_args(0) = f; // the 4 elements are func,  direction, theta, otherargs
	step_args(3) = otherargs;
	
	// Check number of args
	if (args.length() > 4) error("\nNewtonMin: you have provided more than 4 arguments\n");
	if (args.length() < 3) error("\nNewtonMin: you have provided fewer than 3 arguments\n");
	
//	Type checking: not done yet, since I don't know how to do it
// 	// Check types of required arguments
// 	if !ischar(func) error("\nNewtonMin: first argument must be a string that gives name of objective function\n"); endif
// 	if !isvector(theta) error("\nNewtonMin: second argument must be a vector\n"); endif
// 	if !iscell(otherargs) error("\nNewtonMin: third argument must be a cell array\n\
// 	 (even if it's just placeholder such as {1}\n"); endif
//	if (!is_matrix_type(control)) error("\nNewtonMin: If 4 arguments passed, the 4th must be control vector\n"); endif
//	if (control.rows() != 3) error("\nNewtonMin: control (4th argument) must be a 3x1 vector\n"); endif

	// tolerances
	func_tol  = 10*sqrt(2.22e-16);
	param_tol = 1e-5;
	gradient_tol = 1e-4;  

	// Default values for controls
	max_iters = INT_MAX; // no limit on iterations
	verbosity = 0; // by default don't report results to screen
	criterion = 1; // strong convergence required

 	// use provided controls, if applicable
	if (args.length() == 4)
	{
		control = args(3).matrix_value();
		max_iters = (int) control(0);
		if (max_iters == -1) max_iters = INT_MAX;
		verbosity = (int) control(1);
		criterion = (int) control(2);
	}	
	
	// initialize things
	convergence = -1; // if this doesn't change, it means that maxiters were exceeded
 	thetain = theta;
	H = identity_matrix(k,k);

	// Initial obj_value
	f_return = feval(f, f_args); 
	last_obj_value = f_return(0).double_value();

	
	// initial gradient
	f_return = feval("NumGradient", args);
	g = f_return(0).matrix_value();
	g = g.transpose();

	// check that gradient is ok
	gradient_failed = 0;  // test = 1 means gradient failed
	for (i=0; i<k;i++)
	{
		gradient_failed = gradient_failed + isnan(g(i));
	}
	if (gradient_failed)
	{
		error("NewtonMin: Initial gradient could not be calculated: exiting");
		convergence = 2;
	}

	// MAIN LOOP STARTS HERE
 	for (iter = 0; iter < max_iters; iter++)
	{
 		d = -H.inverse() * g;

		// Regular step
		step_args(1) = d;
		step_args(2) = theta;
		f_return = feval("NewtonStep", step_args);
		stepsize = f_return(0).double_value();
		obj_value = f_return(1).double_value();
		
		// Steepest descent if regular step fails
		if (isnan(stepsize)) 
		{
			d = -g; // try steepest descent
			step_args(1) = d;
			f_return = feval("NewtonStep", step_args);
			stepsize = f_return(0).double_value();
			obj_value = f_return(1).double_value();

			warning("NewtonMin: Stepsize failure in Newton direction, trying steepest descent direction");
			// if that didn't work, were out of luck
			if (isnan(stepsize))
			{
				warning("NewtonMin: unable to find direction of improvement: exiting");
				theta = thetain;
				convergence = 2;
				break;
			}
 		}
      	p = stepsize*d;

 		// check normal convergence: all 3 must be satisfied
		// function convergence
		if (fabs(last_obj_value) > 1.0)
		{
 			conv1 = (fabs(obj_value - last_obj_value)/fabs(last_obj_value)) < func_tol;
		}
		else
		{
			conv1 = fabs(obj_value - last_obj_value) < func_tol;
		}	
		// parameter change convergence
		temp = theta.transpose() * theta;
		test = sqrt(temp(0,0));
		if (test > 1)
		{
			temp = p.transpose() * p;
			conv2 = (temp(0,0) / test) < param_tol ;
			
		}
		else
		{
			temp = p.transpose() * p;
			conv2 = temp(0,0) < param_tol;
		}		
		// gradient convergence
		temp = g.transpose() * g;
		conv3 = sqrt(temp(0,0)) < gradient_tol;


		// Want intermediate results?
		if (verbosity == 1)
		{
			printf("\nNewtonMin intermediate results: Iteration %d\n",iter);
			printf("Stepsize %8.7f\n", stepsize);
			printf("Objective function value %16.10f\n", last_obj_value);
			printf("Function conv %d  Param conv %d  Gradient conv %d\n", conv1, conv2, conv3);	
			printf("  params  gradient  change\n");
	 		for (j = 0; j<k; j++)
			{
				printf("%8.4f %8.4f %8.4f\n",theta(j),g(j),p(j));
			}		
			printf("\n");
		}
	
		// Are we done?
		if (criterion == 1)
		{
			if (conv1 && conv2 && conv3)
			{
		  	convergence = 1;
	     	break;
			}
		}		
		else if (conv1)
		{
				convergence = 1;
				break;
		}
		
			
  	    last_obj_value = obj_value;	
		thetanew = theta + p;
		theta = thetanew;

		// get gradient at new parameter value
		g_args(1) = theta;
		f_return = feval("NumGradient", g_args);
		g_new = f_return(0).matrix_value();
		g_new = g_new.transpose();

		// Check that gradient is ok
		gradient_failed = 0;  // test = 1 means gradient failed
		for (i=0; i<k;i++)
		{
			gradient_failed = gradient_failed + isnan(g_new(i));
		}

		// Hessian		
		if (!gradient_failed)
		{
			g = g_new;
			f_return = feval("NumHessian", g_args);
			H = f_return(0).matrix_value();
		}
		else // failed gradient - try Hessian reset
		{
			// if we already tried it, we're out of luck
			if (H == identity_matrix(k,k))
			{
				convergence = 2;
				theta = thetain;
				warning("NewtonMin: failure of gradient: exiting");
				break;
			}	
			H = identity_matrix(k,k);
		}	
	}
	
// 	// Want last iteration results?
	if (verbosity == 2)
	{
		printf("\nNewtonMin final results: Iteration %d\n",iter);
		printf("Stepsize %8.7f\n", stepsize);
		printf("Objective function value %16.10f\n", last_obj_value);
		printf("Function conv %d  Param conv %d  Gradient conv %d\n", conv1, conv2, conv3);	
		printf("  params  gradient  change\n");
		for (j = 0; j<k; j++)
		{
			printf("%8.4f %8.4f %8.4f\n",theta(j),g(j),p(j));
		}		
		printf("\n");
	}
	f_return(0) = theta;
	f_return(1) = obj_value;
	f_return(2) = iter;
	f_return(3) = convergence;
	return octave_value_list(f_return);
}

// ============================= NumHessian =============================
//  Numeric second derivative function of the form
//  f(theta, otherargs), where theta is a vector and otherargs is a cell array.
// 
//  The derivative is respect to the FIRST element of args.
// 

#include <oct.h>
#include <octave/parse.h>
#include <octave/Cell.h>
#include <octave/lo-mappers.h>


// define argument checks
static bool
any_bad_argument(const octave_value_list& args)
{


	if (!args(0).is_string())
    {
        error("NumHessian: first argument must be string holding objective function name");
        return true;
    }
	if (!args(1).is_real_matrix())
    {
        error("NumHessian: second argument must be vector of initial parameter values");
        return true;
    }
	if (!args(2).is_cell())
    {
        error("NumHessian: third argument must a cell of other arguments\n\
		 (or a placeholder such as {1})");
        return true;
    }
    return false;
}

DEFUN_DLD(NumHessian, args, ,
"Usage:\n\
hessian = NumHessian(f, x, otherargs)\n\
\n\
Numeric Hessian matrix of a function f() with respect to x\n\
\n\
f() should be of the form\n\
[f_value, possibly more returns...] = f(x, otherargs)\n\
where otherargs is a cell array\n\
\n\
or if there are no other arguments\n\
[f_value, possibly more returns...] = f(x)\n\
\n\
f_value must be real-valued\n\
\n\
hessian will be dimension kxk, where k=dim(x)\n\
\n\
Inputs:\n\
f: string - the function to be minimized\n\
x: vector - the value to at which gradient is taken\n\
otherargs: cell array - additional arguments to function\n\
	(e.g., data)\n\
	if there are none, use a place holder such as {1}\n\
Outputs:\n\
hessian: kxk matrix\n\
\n\
Example:\n\
\n\
function a = myfunc(x)\n\
	a = x'*x;\n\
endfunction\n\
\n\
x = rand(3,1)\n\
x =\n\
\n\
  0.81472\n\
  0.13548\n\
  0.90579\n\
\n\
h = NumHessian(\"myfunc\", x, {1})\n\
h =\n\
\n\
  2.00000  0.00000  0.00000\n\
  0.00000  2.00000  0.00000\n\
  0.00000  0.00000  2.00000\n\
")
{
	int nargin = args.length();
	if (nargin < 1)
    {
      	error("NumHessian: you must supply 3 arguments\n\
		even if the 3rd is just a placeholder such as {1}");
	    return octave_value_list();
    }
 
    std::string f = args(0).string_value();
    Matrix parameter = args(1).matrix_value();
	octave_value otherargs = args(2);
	octave_value_list f_args(2,1);
	octave_value_list fdiff_args(2,1);
	octave_value_list f_return;

	const int k = parameter.rows();
    Matrix derivative(k, k);
 	
	int i, j;
	double di, hi, pi, dj, hj, pj, hia;
	double hja, fpp, fmm, fmp, fpm, obj_value;
	
	f_args(0) = parameter;
	f_args(1) = otherargs;


	f_return = feval(f, f_args); 
	obj_value = f_return(0).double_value();

 
    for (i = 0; i<k;i++)	// approximate 2nd deriv. by central difference 
    {
	    pi = parameter(i);
		fdiff_args(0) = pi;
		fdiff_args(1) = 2;
		f_return = feval("FiniteDifference", fdiff_args);
		hi = f_return(0).double_value();
        for (j = 0; j < i; j++) // off-diagonal elements
		{
		    pj = parameter(j);
			fdiff_args(0) = pj;
			fdiff_args(1) = 2;
			f_return = feval("FiniteDifference", fdiff_args);
			hj = f_return(0).double_value();
	
            // +1 +1
			parameter(i) = di = pi + hi;
			parameter(j) = dj = pj + hj; 
			hia = di - pi;
			hja = dj - pj;
			f_args(0) = parameter;
 			f_return = feval(f, f_args); 
			fpp = f_return(0).double_value();

			// -1 -1
			parameter(i) = di = pi - hi;
			parameter(j) = dj = pj - hj; 
			hia = hia + pi - di;
			hja = hja + pj - dj;
			f_args(0) = parameter;
 			f_return = feval(f, f_args); 
			fmm = f_return(0).double_value();
			
            // +1 -1
			parameter(i) = pi + hi;
			parameter(j) = pj - hj;
			f_args(0) = parameter;
 			f_return = feval(f, f_args); 
			fpm = f_return(0).double_value();

            // -1 +1 
			parameter(i) = pi - hi;
			parameter(j) = pj + hj;
			f_args(0) = parameter;
 			f_return = feval(f, f_args); 
			fmp = f_return(0).double_value();

            derivative(j,i) = ((fpp - fpm) + (fmm - fmp)) / (hia * hja);
			derivative(i,j) = derivative(j,i);

            parameter(j) = pj;
		}

		// diagonal elements
		// +1 +1  
        parameter(i) = di = pi + 2 * hi;
		f_args(0) = parameter;
 		f_return = feval(f, f_args); 
		fpp = f_return(0).double_value();

		hia = (di - pi) / 2;

		// -1 -1 
        parameter(i) = di = pi - 2 * hi;
		f_args(0) = parameter;
 		f_return = feval(f, f_args); 
		fmm = f_return(0).double_value();
		hia = hia + (pi - di) / 2;

        derivative(i,i) = ((fpp - obj_value) + (fmm - obj_value)) / (hia * hia);

        parameter(i) = pi;
    }

	return octave_value(derivative);
}

# Makefile optimization stuff

.phony: all
all: BisectionStep.oct NewtonStep.oct NumGradient.oct BFGSMin.oct NewtonMin.oct FiniteDifference.oct NumHessian.oct

%.oct: %.cc
	mkoctfile $<

1;
# This shows how to call BFGSMin.m


# This is just used to form objective functions, not important
# Function value and gradient vector of the rosenbrock function
# The minimizer is at the vector (1,1,..,1),
# and the minimized value is 0.
function [obj_value, gradient] = rosenbrock(x);
	dimension = length(x);
	obj_value = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2);
	if nargout > 1
	  gradient = zeros(dimension, 1);
	  gradient(1:dimension-1) = - 400*x(1:dimension-1).*(x(2:dimension)-x(1:dimension-1).^2) - 2*(1-x(1:dimension-1));
	  gradient(2:dimension) = gradient(2:dimension) + 200*(x(2:dimension)-x(1:dimension-1).^2);
	endif
endfunction



# example obj. fn. - illustrates use of "otherargs", and more than 1 return
function [obj_value, junk] = objective(theta, otherargs)
	location = otherargs{1};
	x = theta - location + ones(rows(theta),1); # move minimizer to "location"
	obj_value = rosenbrock(x);
	junk = 1;	
endfunction 


# example obj. fn. - illustrates calling BFGSMin with 2 arguments
function obj_value = objective2(theta)
	# location defined internally, since not passed
	dim = rows(theta) - 1;
	location = 5*(0:dim)/dim;
	location = location';
	x = theta - location + ones(rows(theta),1); # move minimizer to "location"
	obj_value = rosenbrock(x);	
endfunction 



# control options and initial value
control = [1000;2;1];  # max 1000 iterations; only report results of last iteration; strong convergence required
dim = 10; # dimension of Rosenbrock function
theta = zeros(dim+1,1);  # starting values
location = 5*(0:dim)/dim;
location = location';


# do the minimization
printf("If this was successful, the minimizer should be\n");
printf("a vector of even steps from 0 to 5\n\n");

printf("BFGSMin - 4 input args\n");
theta = theta - theta;
t=cputime();
[theta, obj_value, iterations, convergence] = BFGSMin("objective", theta, {location}, control);
t = cputime() - t;
printf("Elapsed time = %f\n",t);


printf("BFGSMin - 4 input args (dummy otherargs, control supplied)\n");
theta = theta - theta;
t=cputime();
[theta, obj_value, iterations, convergence] = BFGSMin("objective2", theta, {1}, control);
t = cputime() - t;
printf("Elapsed time = %f\n",t);

printf("BFGSMin - 3 input args (dummy otherargs, control not supplied)\n");
theta = theta - theta;
t=cputime();
[theta, obj_value, iterations, convergence] = BFGSMin("objective2", theta, {1});
t = cputime() - t;
printf("check that this last worked, since it doesn't print results\n");
theta
printf("Elapsed time = %f\n",t);


printf("NewtonMin - 4 input args\n");
theta = theta - theta;
t=cputime();
[theta, obj_value, iterations, convergence] = NewtonMin("objective", theta, {location}, control);
t = cputime() - t;
printf("Elapsed time = %f\n",t);