[Octave-cvsupdate] octave-forge/main/optim/mintoolkit Makefile,1.2,1.3 bfgsmin.cc,1.1,1.2 bisectionstep.cc,1.1,1.2 finitedifference.cc,1.1,1.2 newtonstep.cc,1.1,1.2 samin.cc,1.1,1.2 Example2SA.m,1.1,NONE ExampleMin.m,1.1,NONE ExampleSA.m,1.1,NONE newtonmin.cc,1.1,NONE

[Michael C. <mc...@us...>]

Update of /cvsroot/octave/octave-forge/main/optim/mintoolkit
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7445/mintoolkit

Modified Files:
	Makefile bfgsmin.cc bisectionstep.cc finitedifference.cc 
	newtonstep.cc samin.cc 
Removed Files:
	Example2SA.m ExampleMin.m ExampleSA.m newtonmin.cc 
Log Message:
update files

Index: bfgsmin.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/bfgsmin.cc,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- bfgsmin.cc	19 Mar 2004 15:51:23 -0000	1.1
+++ bfgsmin.cc	30 Mar 2004 14:37:23 -0000	1.2
@@ -26,6 +26,7 @@
 #include <octave/parse.h>
 #include <octave/Cell.h>
 #include <octave/lo-mappers.h>
+#include <float.h>
 
 // define argument checks
 static bool
@@ -43,9 +44,9 @@
     }
  	if (args.length() == 3)
 	{
-		if (!args(2).is_real_matrix() || (!(args(2).length() == 4)))
+		if (!args(2).is_cell() || (!(args(2).length() == 4)))
     	{
-        	error("bfgsmin: 3rd argument, if supplied,  must a 4x1 vector of control options");
+        	error("bfgsmin: 3rd argument, if supplied,  must a 4-element cell array of control options");
 	        return true;
     	}
 	}	
@@ -54,11 +55,14 @@
 
 
 DEFUN_DLD(bfgsmin, args, ,
- " bfgs minimization of a function.\n\
+"bfgsmin: bfgs minimization of a function.\n\
 \n\
-* first argument is function name (string)\n\
-* second is a cell array that holds all arguments of the function,\n\
-* third argument is a optional 4x1 control vector\n\
+[x, obj, convergence] = bfgsmin(\"f\", {args}, {control})\n\
+\n\
+Arguments:\n\
+* \"f\": function name (string)\n\
+* {args}: a cell array that holds all arguments of the function,\n\
+* {control}: an optional cell array of 4 elements\n\
 	* elem 1: maximum iterations  (-1 = infinity (default))\n\
 	* elem 2: verbosity\n\
 		- 0 = no screen output (default)\n\
@@ -68,7 +72,12 @@
 		 - 1 = strict (function, gradient and param change) (default)\n\
 		 - 2 = weak - only function convergence required\n\
 	* elem 4: arg with respect to which minimization is done\n\
-	(default = first)\n\
+		(default = first)\n\
+\n\
+Returns:\n\
+* x: the minimizer\n\
+* obj: the value of f() at x\n\
+* convergence: 1 if normal conv, other values if not\n\
 \n\
 Example:\n\
 function a = f(x,y)\n\
@@ -123,39 +132,42 @@
 	
 	octave_value_list f_return; // holder for feval returns
 
-	int criterion, max_iters, convergence, verbosity, minarg, iter;
+	int max_iters, verbosity, criterion, minarg;
+	int convergence, iter;
 	int gradient_failed, i, j, k, conv1, conv2, conv3;
-
+	int have_gradient = 0;
 	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;
-	
+	double last_obj_value, denominator, test;
+	Matrix H, tempmatrix, H1, H2;
+	ColumnVector thetain, d, g, g_new, p, q;
 	// tolerances
-	func_tol  = 10*sqrt(2.22e-16);
+	func_tol  = 2*DBL_EPSILON;
 	param_tol = 1e-6;
-	gradient_tol = 1e-5;  
+	gradient_tol = 1e-6;  
 
-	// 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
-	minarg = 1; // by default, first arg is one over which we minimize
- 
  	// use provided controls, if applicable
 	if (args.length() == 3)
 	{
-		control = args(2).matrix_value();
-		max_iters = (int) control(0);
+		Cell control (args(2).cell_value());
+		max_iters = control(0).int_value();
 		if (max_iters == -1) max_iters = INT_MAX;
-		verbosity = (int) control(1);
-		criterion = (int) control(2);
-		minarg = (int) control(3);
+		verbosity = control(1).int_value();
+		criterion = control(2).int_value();
+		minarg = control(3).int_value();
 	}	
-
+	else
+	{
+		// 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
+		minarg = 1; // by default, first arg is one over which we minimize
+ 	}
+	
 	step_args(3) = minarg;
 	g_args(2) = minarg;
 	
- 	Matrix theta  = f_args(minarg - 1).matrix_value();
+ 	ColumnVector theta  = f_args(minarg - 1).column_vector_value();
 
 	k = theta.rows();
 	
@@ -167,11 +179,29 @@
 	// Initial obj_value
 	f_return = feval("celleval", c_args); 
 	last_obj_value = f_return(0).double_value();
-
-	// initial gradient
-	f_return = feval("numgradient", g_args);
-	g = f_return(0).matrix_value();
-	g = g.transpose();
+	
+	// maybe we have analytic gradient?
+	// if the function returns more than one item, and the second
+	// is a kx1 vector, it is assumed to be the gradient
+	if (f_return.length() > 1)
+	{
+		if (f_return(1).is_real_matrix())
+		{
+			if ((f_return(1).rows() == k) & (f_return(1).columns() == 1))
+			{
+				g = f_return(1).column_vector_value();
+				have_gradient = 1; // future reference
+			}	
+		}
+		else have_gradient = 0;
+		
+	}
+	if (!have_gradient) // use numeric gradient
+	{
+		f_return = feval("numgradient", g_args);
+		tempmatrix = f_return(0).matrix_value();
+		g = tempmatrix.row(0).transpose();
+	}
 
 	// check that gradient is ok
 	gradient_failed = 0;  // test = 1 means gradient failed
@@ -188,36 +218,47 @@
 	// MAIN LOOP STARTS HERE
  	for (iter = 0; iter < max_iters; iter++)
 	{
- 		d = -H*g;
 
-		// Regular step
-		step_args(2) = d;
-		f_args(minarg - 1) = theta;
-		step_args(1) = f_args;
-		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 = -H*g;
+		// if direction not zero, get stepsize
+		if (!((fabs(d.max()) < DBL_EPSILON) \
+				|| (d.min() < 0 & (fabs(d.min()) < 2*DBL_EPSILON))))
 		{
-			warning("bfgsmin: Stepsize failure in BFGS direction, trying steepest descent direction");
-			d = -g; // try steepest descent
+ 	
+			// stepsize: try bfgs direction, then steepest descent if it fails
 			step_args(2) = d;
+			f_args(minarg - 1) = theta;
+			step_args(1) = f_args;
 			f_return = feval("newtonstep", step_args);
 			stepsize = f_return(0).double_value();
 			obj_value = f_return(1).double_value();
-
-			// if that didn't work, were out of luck
-			if (xisnan(stepsize))
+			if (stepsize == 0.0) // fall back to steepest descent
 			{
-				warning("bfgsmin: unable to find direction of improvement: exiting");
-				theta = thetain;
-				convergence = 2;
-				break;
+				warning("bfgsmin: Stepsize failure in BFGS direction, trying steepest descent");
+				d = -g; // try steepest descent
+				step_args(2) = d;
+				f_return = feval("newtonstep", step_args);
+				stepsize = f_return(0).double_value();
+				obj_value = f_return(1).double_value();
+				if (stepsize == 0.0) // out of luck on stepsize
+				{
+					warning("bfgsmin: unable to find direction of improvement: exiting");
+					theta = thetain;
+					convergence = 2;
+					break;
+				}
 			}
- 		}
-      	p = stepsize*d;
+		}
+		else // if gradient zero, then check function conv.
+		{
+			f_args(minarg - 1) = theta;
+			c_args(1) = f_args;
+			f_return = feval("celleval", c_args); 
+			obj_value = f_return(0).double_value();
+			stepsize = 1;	
+		}	
+
+		p = stepsize*d;
 
  		// check normal convergence: all 3 must be satisfied
 		// function convergence
@@ -230,29 +271,27 @@
 			conv1 = fabs(obj_value - last_obj_value) < func_tol;
 		}	
 		// parameter change convergence
-		temp = theta.transpose() * theta;
-		test = sqrt(temp(0,0));
+		test = sqrt(theta.transpose() * theta);
 		if (test > 1)
 		{
-			temp = p.transpose() * p;
-			conv2 = (temp(0,0) / test) < param_tol ;
+			conv2 = p.transpose() * p / test < param_tol ;
 			
 		}
 		else
 		{
-			temp = p.transpose() * p;
-			conv2 = temp(0,0) < param_tol;
+			conv2 = p.transpose() * p < param_tol;
 		}		
 		// gradient convergence
-		temp = g.transpose() * g;
-		conv3 = sqrt(temp(0,0)) < gradient_tol;
+		conv3 = sqrt(g.transpose() * g) < gradient_tol;
 
 
 		// Want intermediate results?
 		if (verbosity == 1)
 		{
-			printf("\nbfgsmin intermediate results: Iteration %d\n",iter);
+			printf("\nBFGSMIN intermediate results: Iteration %d\n",iter);
 			printf("Stepsize %8.7f\n", stepsize);
+			if (have_gradient) printf("Using analytic gradient\n");
+			else printf("Using numeric gradient\n");
 			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");
@@ -280,15 +319,23 @@
 		
 			
   	    last_obj_value = obj_value;	
-		thetanew = theta + p;
-		theta = thetanew;
+		theta = theta + p;
 
-		// get gradient at new parameter value
+		// new gradient
 		f_args(minarg - 1) = theta;
-		g_args(1) = f_args;
-		f_return = feval("numgradient", g_args);
-		g_new = f_return(0).matrix_value();
-		g_new = g_new.transpose();
+		if (have_gradient)
+		{
+			c_args(1) = f_args;
+			f_return = feval("celleval", c_args);
+			g_new = f_return(1).column_vector_value();
+		}
+		else // use numeric gradient
+		{
+			g_args(1) = f_args;
+			f_return = feval("numgradient", g_args);
+			tempmatrix = f_return(0).matrix_value();
+			g_new = tempmatrix.row(0).transpose();
+		}
 
 		// Check that gradient is ok
 		gradient_failed = 0;  // test = 1 means gradient failed
@@ -297,48 +344,28 @@
 			gradient_failed = gradient_failed + xisnan(g_new(i));
 		}
 
-		// Hessian update if gradient ok		
+	// Hessian update if gradient ok		
 		if (!gradient_failed)
 		{
 			q = g_new-g;
   	    	g = g_new;
-			temp = q.transpose() * p;
-			denominator = temp(0,0);
-	  	  	if (fabs(denominator) < 2.2e-16)  // reset Hessian if necessary
+			denominator = q.transpose()*p;
+	  	  	if ((fabs(denominator) < DBL_EPSILON))  // reset Hessian if necessary
 			{  	
-				// if we already reset Hessian then we're out of luck
-				if ((H == identity_matrix(k,k)) & !(iter==0))
-				{					 
-					convergence = 2;
-					theta = thetain;
-				  	warning("bfgsmin: %d unable to find direction of improvement: exiting",iter);
-					break;
-				}	
-				// Here's the Hessian re-set
+				if (verbosity == 1) printf("bfgsmin: Hessian reset\n");
 				H = identity_matrix(k,k);
 			}	
-		  	else // normal update
+		  	else 
 			{
-				temp = (1.0+(q.transpose() * H * q) / denominator) / denominator;
-				tempscalar = temp(0,0);
-				H1 = tempscalar * (p * p.transpose());
+				H1 = (1.0+(q.transpose() * H * q) / denominator) / denominator \
+					* (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);
-		}	
+		else H = identity_matrix(k,k); // reset hessian if gradient fails
+		// then try to start again with steepest descent
 	}
 	
  	// Want last iteration results?
@@ -350,6 +377,8 @@
 		if ((convergence == 1) & (criterion == 1)) printf("STRONG CONVERGENCE\n");
 		if ((convergence == 1) & !(criterion == 1)) printf("WEAK CONVERGENCE\n");
 		if (convergence == 2) printf("NO CONVERGENCE: algorithm failed\n");
+		if (have_gradient) printf("Used analytic gradient\n");
+		else printf("Used numeric gradient\n");
 		printf("\nObj. fn. value %f     Iteration %d\n", last_obj_value, iter);
 		printf("Function conv %d  Param conv %d  Gradient conv %d\n", conv1, conv2, conv3);	
 		printf("\n  param  gradient  change\n");
@@ -362,7 +391,6 @@
 	}
 	f_return(0) = theta;
 	f_return(1) = obj_value;
-	f_return(2) = iter;
-	f_return(3) = convergence;
+	f_return(2) = convergence;
 	return octave_value_list(f_return);
 }

Index: bisectionstep.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/bisectionstep.cc,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- bisectionstep.cc	19 Mar 2004 15:51:23 -0000	1.1
+++ bisectionstep.cc	30 Mar 2004 14:37:23 -0000	1.2
@@ -15,6 +15,10 @@
 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 
 // bisection step, for use by minimization algorithms
+// 
+// returns optimal stepsize using bisection - first until an improvement
+// is found, then until there is no further improvement
+//
 // since this is not for use by end-users, its documentation
 // is a bit short. See bfgsmin for an example of use
 // args are
@@ -26,19 +30,20 @@
 #include <oct.h>
 #include <octave/parse.h>
 #include <octave/Cell.h>
+#include <float.h>
 
 DEFUN_DLD(bisectionstep, args, , "bisectionstep.cc")
 {
    	std::string f (args(0).string_value());
 	Cell f_args (args(1).cell_value());
-	Matrix dx (args(2).matrix_value());
+	ColumnVector dx (args(2).column_vector_value());
 
 	double obj_0, obj, a;
 	octave_value_list f_return;
 	octave_value_list c_args(2,1); // for cellevall {f, f_args}  
 
 	octave_value_list stepobj(2,1);
-	int minarg;
+	int minarg, found_improvement;
 
 	// Default values for controls
 	minarg = 1; // by default, first arg is one over which we minimize
@@ -49,8 +54,8 @@
 		minarg = args(3).int_value();
 	}	
  
-	Matrix x (f_args(minarg - 1).matrix_value());
-	Matrix x_in = x;
+	ColumnVector x (f_args(minarg - 1).column_vector_value());
+	ColumnVector x_in = x;
 	
 
 	// possibly function returns a cell array
@@ -61,9 +66,9 @@
 	obj_0 = f_return(0).double_value();
 
 	a = 1.0;
-
+	found_improvement = 0;
 	// this first loop goes until an improvement is found
-  	while (a > 4e-16) // limit iterations
+  	while (a > 2*DBL_EPSILON) // limit iterations
 	{
 		f_args(minarg - 1) = x + a*dx;
 		c_args(1) = f_args;
@@ -71,19 +76,28 @@
 		obj = f_return(0).double_value();
  
  		// reduce stepsize if worse, or if function can't be evaluated
-		if ((obj > obj_0) || isnan(obj))
+		if ((obj >= obj_0) || isnan(obj))
 		{
 			a = 0.5 * a;
 		}	
 		else
 		{
 			obj_0 = obj;
+			found_improvement = 1;
 			break;
 		}
 	}
 	
+	// If unable to find any improvement break out with stepsize zero
+	if (!found_improvement)
+	{
+		stepobj(0) = 0.0;
+		stepobj(1) = obj_0;
+		return octave_value_list(stepobj);
+	}	
+	
 	// now keep going until we no longer improve, or reach max trials
-	while (a > 4e-16)
+	while (a > 2*DBL_EPSILON)
 	{
 	   	a = 0.5*a; 
 		f_args(minarg - 1) = x + a*dx;

Index: finitedifference.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/finitedifference.cc,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- finitedifference.cc	19 Mar 2004 15:51:23 -0000	1.1
+++ finitedifference.cc	30 Mar 2004 14:37:23 -0000	1.2
@@ -17,6 +17,7 @@
 // ========================= finitedifference ==========================
 //  finite differences for numeric differentiation
 #include <oct.h>
+#include <float.h>
 DEFUN_DLD(finitedifference, args, ,"finitedifference, C++ version\n\
 differences for numgradient and numhessian")
 {
@@ -25,7 +26,7 @@
 	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
+	eps = DBL_EPSILON; // machine precision
 	SQRT_EPS = sqrt(eps); 
 	DIFF_EPS = exp(log(eps)/2);
 	DIFF_EPS1 = exp(log(eps)/3);

--- newtonmin.cc DELETED ---

--- Example2SA.m DELETED ---

Index: newtonstep.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/newtonstep.cc,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- newtonstep.cc	19 Mar 2004 15:51:23 -0000	1.1
+++ newtonstep.cc	30 Mar 2004 14:37:23 -0000	1.2
@@ -25,6 +25,7 @@
 #include <oct.h>
 #include <octave/parse.h>
 #include <octave/Cell.h>
+#include <float.h>
 
 DEFUN_DLD(newtonstep, args, , "newtonstep.cc")
 {
@@ -83,15 +84,28 @@
   	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
+	if (hessian <= 2*DBL_EPSILON) hessian = 1.0; // avoid div by zero
 
 	a = - gradient / hessian;  // hessian inverse gradient: the Newton step
+
+	if (a < 0) 	// since direction is descending, a must be positive
+	{ 			// if it is not, go to bisection step
+		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);
+ 	}
+
 	a = (a < 5.0)*a + 5.0*(a>=5.0); // Let's avoid extreme steps that might cause crashes
-	// check that this is improvement
+
+	// ensure that this is improvement
 	f_args(minarg - 1) = x + a*dx;
 	c_args(1) = f_args;
 	f_return = feval("celleval", c_args); 
 	obj = f_return(0).double_value();
+
 	// if not, fall back to bisection
 	if ((obj > obj_0) || isnan(obj))
 	{

Index: Makefile
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/Makefile,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -d -r1.2 -r1.3
--- Makefile	22 Mar 2004 08:16:56 -0000	1.2
+++ Makefile	30 Mar 2004 14:37:23 -0000	1.3
@@ -1,5 +1,5 @@
 include ../../../Makeconf
 
-all: bfgsmin.oct bisectionstep.oct celleval.oct finitedifference.oct newtonmin.oct newtonstep.oct numgradient.oct numhessian.oct samin.oct
+all: bfgsmin.oct bisectionstep.oct celleval.oct finitedifference.oct lbfgsmin.oct newtonstep.oct numgradient.oct numhessian.oct samin.oct
 
 clean: ; -$(RM) *.o core octave-core *.oct *~

Index: samin.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/optim/mintoolkit/samin.cc,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- samin.cc	19 Mar 2004 15:51:23 -0000	1.1
+++ samin.cc	30 Mar 2004 14:37:23 -0000	1.2
@@ -65,50 +65,55 @@
         return true;
 	}	
 	Cell control (args(2).cell_value());
-	if (!(control.length()==10))
+	if (!(control.length() == 10))
 	{
 	 	error("samin: third argument must be a cell array with 10 elements");
  		return true;
 	}
 	if (!control(0).is_numeric_type())
 	{
-        error("samin: 1st element of controls must be a vector of lower bounds");
+        error("samin: 1st element of controls must be LB: a vector of lower bounds");
         return true;
 	}
 	if (!control(1).is_numeric_type())
 	{
-        error("samin: 2nd element of controls must be a vector of lower bounds");
+        error("samin: 2nd element of controls must be UB: a vector of upper bounds");
         return true;
 	}
 	if (!control(2).is_real_scalar())
 	{
-        error("samin: 3rd element of controls must be integer NT");
+        error("samin: 3rd element of controls must be NT: positive integer\n\
+		 	loops per temperature reduction");
         return true;
 	}
 	if (!control(3).is_real_scalar())
 	{
-        error("samin: 4th element of controls must integer NS");
+        error("samin: 4th element of controls must be NS: positive integer\n\
+		 	loops per stepsize adjustment");
         return true;
 	}
 	if (!control(4).is_real_scalar())
 	{
-        error("samin: 5th element of controls must be integer neps");
+        error("samin: 5th element of controls must be RT:\n\
+		 	temperature reduction factor, between 0 and 1");
         return true;
 	}
 	if (!control(5).is_real_scalar())
 	{
-        error("samin: 6th element of controls must be number \n\
-		 	temperature reduction factor RT, between 0 and 1");
+        error("samin: 6th element of controls must be integer MAXEVALS ");
         return true;
 	}
 	if (!control(6).is_real_scalar())
 	{
-        error("samin: 7th element of controls must be integer MAXEVALS ");
+        error("samin: 7th element of controls must be NEPS: positive integer\n\
+		 	number of final obj. values that must be within EPS of eachother ");
         return true;
 	}
+
 	if (!control(7).is_real_scalar())
 	{
-        error("samin: 8th element of controls must must be TOL (> 0)");
+        error("samin: 8th element of controls must must be EPS (> 0)\n\
+			used to compare the last NEPS obj values for convergence test");
         return true;
 	}
 	if (!control(8).is_real_scalar())
@@ -128,12 +133,14 @@
 
 //-------------- The annealing algorithm --------------
 DEFUN_DLD(samin, args, ,
-  "simulated annealing minimization of a function.\n\
+"samin: simulated annealing minimization of a function.\n\
 \n\
-* first argument is function name (string)\n\
-* second is a cell array that holds all arguments of the function,\n\
-* third argument is a cell array with 10 elements that controls\n\
-	the algorithm\n\
+[x, obj, convergence] = samin(\"f\", {args}, {control})\n\
+\n\
+Arguments:\n\
+* \"f\": function name (string)\n\
+* {args}: a cell array that holds all arguments of the function,\n\
+* {control}: a cell array with 10 elements\n\
 	* LB  - vector of lower bounds\n\
 	* UB - vector of upper bounds\n\
 	* nt - integer: # of iterations between temperature reductions\n\
@@ -148,6 +155,11 @@
 		* 2 = only final results to screen\n\
 	* minarg - integer: which of function args is minimization over?\n\
 \n\
+Returns:\n\
+* x: the minimizer\n\
+* obj: the value of f() at x\n\
+* convergence: 1 if normal conv, other values if not\n\
+\n\
 Example: A really stupid way to calculate pi\n\
 function a = f(x)\n\
 	a = cos(x) + 0.01*(x-pi)^2;\n\
@@ -156,9 +168,6 @@
 Set up the controls:\n\
 ub = 20;\n\
 lb = -ub;\n\
-nt = 20;\n\
-ns = 5;\n\
-neps = 5;\n\
 rt = 0.5;\n\
 maxevals = 1e10;\n\
 eps = 1e-5;\n\
@@ -201,24 +210,24 @@
 	octave_value_list c_args(2,1); // for cellevall {f, f_args}  
 	octave_value_list f_return; // holder for feval returns
 
-	int m, i, j, h, n, nacc, nobds, func_evals;
+	int m, i, j, h, n, nacc, func_evals;
  	int nup, nrej, nnew, ndown, lnobds;
 	int converge, test;
-
+	
  	// user provided controls
 	const Matrix lb (control(0).matrix_value());
 	const Matrix ub (control(1).matrix_value());
-	const int nt (control(2).int_value());  
+	const int nt (control(2).int_value());
 	const int ns (control(3).int_value());
-	const int neps (control(4).int_value());
-	const double maxevals (control(6).double_value());
+	const double rt (control(4).double_value());
+	const double maxevals (control(5).double_value());
+	const int neps (control(6).int_value());
 	const double eps (control(7).double_value());
 	const int verbosity (control(8).int_value());
 	const int minarg (control(9).int_value());	
 	
 	
 	double f, fp, p, pp, fopt, rand_draw, ratio, t;
-	double rt;
  	Matrix x  = f_args(minarg - 1).matrix_value();
 	Matrix bounds = ub - lb;
 
@@ -229,7 +238,6 @@
     
     //  Set initial values  
     nacc = 0;
-    nobds = 0;
     func_evals = 0;
 	
  	Matrix fstar(neps,1);  
@@ -258,8 +266,7 @@
 
 	// First stage: find initial temperature so that
 	// bounds grow to cover parameter space
-	t = 1;
-	rt = 2; // double temperature to make bounds grow
+	t = 1000;
 	converge = 0;    
 	while(converge==0)
 	{	
@@ -288,8 +295,7 @@
 	                {
 						xp(h) = lb(h) + (ub(h) - lb(h)) * rand_draw;
                         lnobds = lnobds + 1;
-                        nobds = nobds + 1;
-                   	}  
+                    }  
 
 	                // Evaluate function at new point
 					f_args(minarg - 1) = xp;
@@ -304,8 +310,24 @@
 	                if(func_evals >= maxevals)
 	                {
 						warning("samin: NO CONVERGENCE: MAXEVALS exceeded before initial temparature found");
+                		if(verbosity >= 1)
+						{
+							printf("\033[2J");
+							printf("\n================================================\n");
+							printf("SAMIN results\n");
+							printf("NO CONVERGENCE: MAXEVALS exceeded\n");
+							printf("Stage 1, increasing temperature\n");
+							printf("\nObj. fn. value %f\n", fopt);
+  							printf("           parameter        search width\n");
+		    				for(i = 0; i < n; i++)
+							{
+								printf("%20f%20f\n", xopt(i), bounds(i));
+							}
+							printf("================================================\n");
+	    				} 
                 		f_return(0) = xopt;
 						f_return(1) = fopt;
+						f_return(2) = 0;
 						return octave_value_list(f_return);
 					}
                 
@@ -366,7 +388,7 @@
 	                bounds(i) = bounds(i) / (1.0 + 2.0 * ((0.4 - ratio) / 0.4));
 	            }
 				// keep within initial bounds
-	            if(bounds(i) > (ub(i) - lb(i)))
+	            if(bounds(i) >= (ub(i) - lb(i)))
 	            {
 	                test = test + 1; // when this gets to n, we're done with fist stage
 					bounds(i) = ub(i) - lb(i);
@@ -378,11 +400,12 @@
     
 	    if(verbosity == 1)
 	    {
-	    	feval("clc");
-			printf("First stage: Increasing temperature to cover parameter space\n");
+	    	printf("\033[2J");
+			printf("\nFirst stage: Increasing temperature to cover parameter space\n");
 		    printf("\nTemperature  %e", t);
 			printf("\nmin function value so far %f", fopt);
-		    printf("\ntotal moves  %d", nup + ndown + nrej);
+			printf("\ntotal evaluations so far %d", func_evals);
+		    printf("\ntotal moves since temp change %d", nup + ndown + nrej);
 		    printf("\ndownhill  %d", nup);
 		    printf("\naccepted uphill %d", ndown);
 		    printf("\nrejected uphill %d", nrej);
@@ -397,7 +420,7 @@
 	    }
     
 	    // Increase temperature
-	    t = rt * t;
+	    t = 10 * t;
 		for(i = neps-1; i > 0; i--)
 		{
 		fstar(i) = fstar(i-1);
@@ -409,7 +432,6 @@
 
 	
 	// Second stage: temperature reduction loop
-	rt = control(5).double_value(); // set to provided value, not that t has been found
 	converge = 0;    
 	while(converge==0)
 	{
@@ -438,8 +460,8 @@
 	                {
 						xp(h) = lb(h) + (ub(h) - lb(h)) * rand_draw;
                         lnobds = lnobds + 1;
-                        nobds = nobds + 1;
-                   	}  
+					}
+                 
 
 	                // Evaluate function at new point
 					f_args(minarg - 1) = xp;
@@ -454,8 +476,24 @@
 	                if(func_evals >= maxevals)
 	                {
 						warning("samin: NO CONVERGENCE to tolerance %f: maxevals exceeded", eps);
-                		f_return(0) = xopt;
+                		if(verbosity >= 1)
+						{
+							printf("\033[2J");
+							printf("\n================================================\n");
+							printf("SAMIN results\n");
+							printf("NO CONVERGENCE: MAXEVALS exceeded\n");
+							printf("Stage 2, decreasing temperature\n");
+							printf("\nObj. fn. value %f\n", fopt);
+  							printf("           parameter        search width\n");
+		    				for(i = 0; i < n; i++)
+							{
+								printf("%20f%20f\n", xopt(i), bounds(i));
+							}
+							printf("================================================\n");
+	    				}						
+						f_return(0) = xopt;
 						f_return(1) = fopt;
+						f_return(2) = 0;
 						return octave_value_list(f_return);
 					}
                 
@@ -525,11 +563,12 @@
     
 	    if(verbosity == 1)
 	    {
-	    	feval("clc");
-			printf("Intermediate results before next temperature reduction\n");
+	    	printf("\033[2J");
+			printf("\nIntermediate results before next temperature reduction\n");
 		    printf("\nTemperature  %e", t);
 			printf("\nmin function value so far %f", fopt);
-		    printf("\ntotal moves  %d", nup + ndown + nrej);
+			printf("\ntotal evaluations so far %d", func_evals);
+		    printf("\ntotal moves since last temp reduction  %d", nup + ndown + nrej);
 		    printf("\ndownhill  %d", nup);
 		    printf("\naccepted uphill %d", ndown);
 		    printf("\nrejected uphill %d", nrej);
@@ -557,18 +596,35 @@
 		test = (test > 0);
 	    if( ((fopt - fstar(0)) <= eps) && (!test))
 		{
-			converge = 1;
+	        // check for bound narrow enough for parameter convergence
+	        for(i = 0;i < n;i++)
+	        {
+	            if(bounds(i) > 1e-3) test = 0; // no conv. if bounds too wide
+	        }
 		}
-    
+		else test = 0;
+		
+		if (test)
+		{
+			converge = 1;
+			if (lnobds > 0) converge = 2;
+    	}
 	    // Are we done yet?    
-	    if(converge)
+	    if(converge>0)
 	    {
 	       	if(verbosity >= 1)
 			{
-				feval("clc");
+				printf("\033[2J");
 				printf("\n================================================\n");
 				printf("SAMIN final results\n");
-				printf("Sucessful convergence to tolerance %f\n", eps);
+				if (converge == 1) printf("NORMAL CONVERGENCE to tolerance %f\n", eps);
+				if (converge == 2)
+				{
+					printf("CONVERGENCE to bound of parameter space\n");
+					printf("WARNING: %d recent trials out of bounds\n", lnobds);
+					printf("consider expanding bounds and re-running\n");
+				}	
+				printf("\ntotal evaluations %d", func_evals);
 				printf("\nObj. fn. value %f\n", fopt);
   				printf("           parameter        search width\n");
 		    	for(i = 0; i < n; i++)
@@ -579,6 +635,8 @@
 	    	}
 			f_return(0) = xopt;
 			f_return(1) = fopt;
+			if (lnobds > 0) converge = 2;
+			f_return(2) = converge;
 			return octave_value_list(f_return); 
 		}
     
@@ -594,6 +652,7 @@
 	}
 	f_return(0) = xopt;
 	f_return(1) = fopt;
+	f_return(2) = converge;
 	return octave_value_list(f_return);
 }
 

--- ExampleSA.m DELETED ---

--- ExampleMin.m DELETED ---