Mercurial > hg > octave-lyh

changeset 8014:44d206ae68c9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve fsolve compatibility
author John W. Eaton <jwe@octave.org>
date Wed, 06 Aug 2008 15:50:44 -0400
parents b3e667f1ab4c
children 30629059b72d
files NEWS src/ChangeLog src/DLD-FUNCTIONS/fsolve.cc
diffstat 3 files changed, 154 insertions(+), 36 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS
+++ b/NEWS
@@ -12,6 +12,11 @@
     removed from Octave. These function were incompatible with the high
     level graphics handle code.
 
+ ** The fsolve function now accepts an option structure argument (see
+    also the optimset function).
+    The INFO values returned from fsolve have changed to be compatible
+    with Matlab's fsolve function.
+
  ** Object Oriented Programming
 
     TO BE WRITTEN
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,10 @@
 2008-08-06  John W. Eaton  <jwe@octave.org>
 
+	* DLD-FUNCTIONS/fsolve.cc (hybrd_info_to_fsolve_info):
+	Update INFO values to be compatible with Matlab's current fsolve.
+	(make_unimplemented_options, override_options): New functions.
+	(Ffsolve): Handle optimset options.  Update doc string.
+
 	* ov-usr-fcn.cc (octave_user_function::do_multi_index_op,
 	octave_user_script::do_multi_index_op):
 	Call octave_call_stack::backtrace_error_message.
--- a/src/DLD-FUNCTIONS/fsolve.cc
+++ b/src/DLD-FUNCTIONS/fsolve.cc
@@ -25,6 +25,8 @@
 #include <config.h>
 #endif
 
+#include <algorithm>
+#include <set>
 #include <string>
 
 #include <iomanip>
@@ -37,6 +39,7 @@
 #include "defun-dld.h"
 #include "error.h"
 #include "gripes.h"
+#include "oct-map.h"
 #include "oct-obj.h"
 #include "ov-fcn.h"
 #include "ov-cell.h"
@@ -69,29 +72,23 @@
   switch (info)
     {
     case -1:
-      break;
-
-    case 0:
-      info = -2;
-      break;
-
     case 1:
       break;
 
     case 2:
-      info = 4;
+      info = 0;
       break;
 
     case 3:
     case 4:
     case 5:
-      info = 3;
+      info = -2;
       break;
 
     default:
       {
 	std::ostringstream buf;
-	buf << "fsolve: unrecognized value of INFO from MINPACK (= "
+	buf << "fsolve: unexpected value of INFO from MINPACK (= "
 	    << info << ")";
 	std::string msg = buf.str ();
 	warning (msg.c_str ());
@@ -201,6 +198,98 @@
   return retval;
 }
 
+static std::set<std::string>
+make_unimplemented_options (void)
+{
+  static bool initialized = false;
+
+  static std::set<std::string> options;
+
+  if (! initialized)
+    {
+      initialized = true;
+
+      options.insert ("largescale");
+      options.insert ("derivativecheck");
+      options.insert ("diagnostics");
+      options.insert ("diffmaxchange");
+      options.insert ("diffminchange");
+      options.insert ("display");
+      options.insert ("funvalcheck");
+      options.insert ("jacobian");
+      options.insert ("maxfunevals");
+      options.insert ("maxiter");
+      options.insert ("outputfcn");
+      options.insert ("plotfcns");
+      options.insert ("tolfun");
+      options.insert ("tolx");
+      options.insert ("typicalx");
+      options.insert ("jacobmult");
+      options.insert ("jacobpattern");
+      options.insert ("maxpcgiter");
+      options.insert ("precondbandwidth");
+      options.insert ("tolpcg");
+      options.insert ("nonleqnalgorithm");
+      options.insert ("linesearchtype");
+    }
+
+  return options;
+}
+
+static void
+override_options (NLEqn_options& opts, const Octave_map& option_map)
+{
+  string_vector keys = option_map.keys ();
+
+  for (octave_idx_type i = 0; i < keys.length (); i++)
+    {
+      std::string key = keys(i);
+      std::transform (key.begin (), key.end (), key.begin (), tolower);
+
+      if (key == "tolx")
+	{
+	  Cell c = option_map.contents (key);
+
+	  if (c.numel () == 1)
+	    {
+	      octave_value val = c(0);
+
+	      if (! val.is_empty ())
+		{
+		  double dval = val.double_value ();
+
+		  if (! error_state)
+		    opts.set_tolerance (dval);
+		  else
+		    gripe_wrong_type_arg ("fsolve", val);
+		}
+	    }
+	  else
+	    error ("fsolve: invalid value for %s option", key.c_str ());
+	}
+      else
+	{
+	  static std::set<std::string> unimplemented_options
+	    = make_unimplemented_options ();
+
+	  if (unimplemented_options.find (key) != unimplemented_options.end ())
+	    {
+	      Cell c = option_map.contents (key);
+
+	      if (c.numel () == 1)
+		{
+		  octave_value val = c(0);
+
+		  if (! val.is_empty ())
+		    warning_with_id ("Octave:fsolve-unimplemented option",
+				     "fsolve: option `%s' not implemented",
+				     key.c_str ());
+		}
+	    }
+	}
+    }
+}
+
 #define FSOLVE_ABORT() \
   do \
     { \
@@ -227,7 +316,7 @@
 
 DEFUN_DLD (fsolve, args, nargout,
   "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{x}, @var{fval}, @var{info}] =} fsolve (@var{fcn}, @var{x0})\n\
+@deftypefn {Loadable Function} {[@var{x}, @var{fval}, @var{info}] =} fsolve (@var{fcn}, @var{x0}, @var{options})\n\
 Given @var{fcn}, the name of a function of the form @code{f (@var{x})}\n\
 and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\
 equations such that @code{f(@var{x}) == 0}.\n\
@@ -237,17 +326,15 @@
 \n\
 @table @asis\n\
 \n\
-@item -2\n\
-Invalid input parameters.\n\
+@item 1\n\
+Algorithm converged with relative error between two consecutive iterates\n\
+less than or equal to the specified tolerance (see @code{fsolve_options}).\n\
+@item 0\n\
+Iteration limit exceeded.\n\
 @item -1\n\
 Error in user-supplied function.\n\
-@item 1\n\
-Relative error between two consecutive iterates is at most the\n\
-specified tolerance (see @code{fsolve_options}).\n\
-@item 3\n\
+@item -2\n\
 Algorithm failed to converge.\n\
-@item 4\n\
-Limit on number of function calls reached.\n\
 @end table\n\
 \n\
 If @var{fcn} is a two-element string array, or a two element cell array\n\
@@ -269,6 +356,11 @@
 \n\
 You can use the function @code{fsolve_options} to set optional\n\
 parameters for @code{fsolve}.\n\
+\n\
+If the optional argument @var{options} is provided, it is expected to\n\
+be a structure of the form returned by @code{optimset}.  Options\n\
+specified in this structure override any set globally by\n\
+@code{optimset, fsolve_options}.\n\
 @end deftypefn")
 {
   octave_value_list retval;
@@ -286,7 +378,7 @@
 
   int nargin = args.length ();
 
-  if (nargin == 2 && nargout < 4)
+  if ((nargin == 2 || nargin == 3) && nargout < 4)
     {
       std::string fcn_name, fname, jac_name, jname;
       fsolve_fcn = 0;
@@ -339,7 +431,7 @@
 	    FSOLVE_ABORT1 ("incorrect number of elements in cell array");
 	}
 
-      if (!fsolve_fcn && ! f_arg.is_cell())
+      if (! fsolve_fcn && ! f_arg.is_cell())
 	{
 	  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
 	    fsolve_fcn = f_arg.function_value ();
@@ -395,7 +487,7 @@
 		}
 	    }
 	}
-
+      
       if (error_state || ! fsolve_fcn)
 	FSOLVE_ABORT ();
 
@@ -417,26 +509,42 @@
 	nleqn_fcn.set_jacobian_function (fsolve_user_jacobian);
 
       NLEqn nleqn (x, nleqn_fcn);
-      nleqn.set_options (fsolve_opts);
 
-      octave_idx_type info;
-      ColumnVector soln = nleqn.solve (info);
+      NLEqn_options local_fsolve_opts (fsolve_opts);
 
-      if (fcn_name.length())
-	clear_function (fcn_name);
-      if (jac_name.length())
-	clear_function (jac_name);
+      if (nargin > 2)
+	{
+	  Octave_map option_map = args(2).map_value ();
+
+	  if (! error_state)
+	    override_options (local_fsolve_opts, option_map);
+	  else
+	    error ("fsolve: expecting optimset structure as third argument");
+	}
 
       if (! error_state)
 	{
-	  retval(2) = static_cast<double> (hybrd_info_to_fsolve_info (info));
-	  retval(1) = nleqn.function_value ();
-	  retval(0) = NDArray (ArrayN<double> (soln.reshape (x_dims)));
+	  nleqn.set_options (local_fsolve_opts);
+
+	  octave_idx_type info;
+	  ColumnVector soln = nleqn.solve (info);
+
+	  if (fcn_name.length())
+	    clear_function (fcn_name);
+	  if (jac_name.length())
+	    clear_function (jac_name);
 
-	  if (! nleqn.solution_ok () && nargout < 2)
+	  if (! error_state)
 	    {
-	      std::string msg = nleqn.error_message ();
-	      error ("fsolve: %s", msg.c_str ());
+	      retval(2) = static_cast<double> (hybrd_info_to_fsolve_info (info));
+	      retval(1) = nleqn.function_value ();
+	      retval(0) = NDArray (ArrayN<double> (soln.reshape (x_dims)));
+
+	      if (! nleqn.solution_ok () && nargout < 2)
+		{
+		  std::string msg = nleqn.error_message ();
+		  error ("fsolve: %s", msg.c_str ());
+		}
 	    }
 	}
     }
@@ -463,7 +571,7 @@
 %! 2.005014 ];
 %! tol = 1.0e-5;
 %! [x, fval, info] = fsolve ("f", [ 0.5; 2.0; 2.5 ]);
-%! info_bad = (info != 1);
+%! info_bad = (info <= 0);
 %! solution_bad = sum (abs (x - x_opt) > tol);
 %! value_bad = sum (abs (fval) > tol);
 %! if (info_bad)
@@ -497,7 +605,7 @@
 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
 %! tol = 1.0e-5;
 %! [x, fval, info] = fsolve ("f", [-1, 1, 2, -1]);
-%! info_bad = (info != 1);
+%! info_bad = (info <= 0);
 %! solution_bad = sum (abs (x - x_opt) > tol);
 %! value_bad = sum (abs (fval) > tol);
 %! if (info_bad)