Mercurial > forge

changeset 8720:f80846665d3b octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use Fotios' complex step derivative approximation. And a few small fixes.
author i7tiol
date Tue, 01 Nov 2011 13:17:41 +0000
parents b87070d83277
children c159aea4f7df
files main/optim/DESCRIPTION main/optim/PKG_ADD main/optim/inst/nonlin_residmin.m main/optim/inst/private/__collect_constraints__.m main/optim/inst/private/__nonlin_residmin__.m main/optim/inst/private/__residmin_stat__.m main/optim/inst/residmin_stat.m
diffstat 7 files changed, 87 insertions(+), 35 deletions(-) [+]
line wrap: on
line diff
--- a/main/optim/DESCRIPTION	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/DESCRIPTION	Tue Nov 01 13:17:41 2011 +0000
@@ -1,6 +1,6 @@
 Name: Optim
 Version: 1.0.17
-Date: 2011-05-20
+Date: 2011-11-01
 Author: Various Authors
 Maintainer: The Octave Community
 Title: Optimzation.
--- a/main/optim/PKG_ADD	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/PKG_ADD	Tue Nov 01 13:17:41 2011 +0000
@@ -1,6 +1,6 @@
 if (compare_versions (version (), "3.3.55", ">="))
   ## optimset mechanism was fixed for option names with underscores
   ## sometime in 3.3.54+, if I remember right
-  __all_opts__ ("__nonlin_residmin__");
-  __all_opts__ ("__residmin_stat__");
+  __all_opts__ ("nonlin_residmin");
+  __all_opts__ ("residmin_stat");
 endif
--- a/main/optim/inst/nonlin_residmin.m	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/inst/nonlin_residmin.m	Tue Nov 01 13:17:41 2011 +0000
@@ -91,6 +91,18 @@
 ## intervals should be used by jacobian functions performing finite
 ## differencing. Default: @code{false (size (parameters))}.
 ##
+## @code{complex_step_derivative}, @code{complex_step_derivative_inequc},
+## @code{complex_step_derivative_equc}: logical scalars, default: false.
+## Estimate Jacobian of model function, general inequality constraints,
+## and general equality constraints, respectively, with complex step
+## derivative approximation. Use only if you know that your model
+## function, function of general inequality constraints, or function of
+## general equality constraints, respectively, is suitable for this. No
+## user function for the respective Jacobian must be specified.
+##
+## @code{cstep}: scalar step size for complex step derivative
+## approximation. Default: 1e-20.
+##
 ## @code{fixed}: logical column vector indicating which parameters
 ## should not be optimized, but kept to their inital value. Fixing is
 ## done independently of the backend, but the backend may choose to fix
--- a/main/optim/inst/private/__collect_constraints__.m	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/inst/private/__collect_constraints__.m	Tue Nov 01 13:17:41 2011 +0000
@@ -14,7 +14,7 @@
 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
 
 function [mc, vc, f_gencstr, df_gencstr, user_df] = \
-      __collect_constraints__ (cstr)
+      __collect_constraints__ (cstr, do_cstep, context)
 
   mc = vc = f_gencstr = df_gencstr = [];
   user_df = false;
@@ -59,13 +59,21 @@
 	tf_gencstr (f_gencstr, varargin{:});
 
     if (user_df)
+      if (do_cstep)
+	error ("both complex step derivative chosen and user Jacobian function specified for %s", context);
+      endif
       if (ischar (df_gencstr))
 	df_gencstr = str2func (df_gencstr);
       endif
       df_gencstr = @ (p, func, idx, hook) \
 	  df_gencstr (p, idx, hook);
     else
-      df_gencstr = @ (p, func, idx, hook) __dfdp__ (p, func, hook);
+      if (do_cstep)
+	df_gencstr = @ (p, func, idx, hook) jacobs (p, func, hook);
+      else
+	__dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
+	df_gencstr = @ (p, func, idx, hook) __dfdp__ (p, func, hook);
+      endif
     endif
   endif
 
--- a/main/optim/inst/private/__nonlin_residmin__.m	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/inst/private/__nonlin_residmin__.m	Tue Nov 01 13:17:41 2011 +0000
@@ -36,6 +36,7 @@
   ## NA here in the frontend
   diffp_default = .001;
   stol_default = .0001;
+  cstep_default = 1e-20;
 
   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
     p = optimset ("param_config", [], \
@@ -55,6 +56,10 @@
 		  "fract_prec", [], \
 		  "diffp", [], \
 		  "diff_onesided", [], \
+		  "complex_step_derivative", false, \
+		  "complex_step_derivative_inequc", false, \
+		  "complex_step_derivative_equc", false, \
+		  "cstep", cstep_default, \
 		  "fixed", [], \
 		  "inequc", [], \
 		  "equc", [], \
@@ -107,6 +112,14 @@
   diffp = optimget (settings, "diffp");
   diff_onesided = optimget (settings, "diff_onesided");
   fixed = optimget (settings, "fixed");
+  do_cstep = optimget (settings, "complex_step_derivative", false);
+  cstep = optimget (settings, "cstep", cstep_default);
+  if (do_cstep && ! isempty (dfdp))
+    error ("both 'complex_step_derivative' and 'dfdp' are set");
+  endif
+  do_cstep_inequc = \
+      optimget (settings, "complex_step_derivative_inequc", false);
+  do_cstep_equc = optimget (settings, "complex_step_derivative_equc", false);
 
   any_vector_conf = ! (isempty (lbound) && isempty (ubound) && \
 		       isempty (max_fract_change) && \
@@ -115,9 +128,11 @@
 
   ## collect constraints
   [mc, vc, f_genicstr, df_gencstr, user_df_gencstr] = \
-      __collect_constraints__ (optimget (settings, "inequc"));
+      __collect_constraints__ (optimget (settings, "inequc"), \
+			       do_cstep_inequc, "inequality constraints");
   [emc, evc, f_genecstr, df_genecstr, user_df_genecstr] = \
-      __collect_constraints__ (optimget (settings, "equc"));
+      __collect_constraints__ (optimget (settings, "equc"), \
+			       do_cstep_equc, "equality constraints");
   mc_struct = isstruct (mc);
   emc_struct = isstruct (emc);
 
@@ -440,8 +455,12 @@
 
   ## jacobian of model function
   if (isempty (dfdp))
-    __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
-    dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+    if (do_cstep)
+      dfdp = @ (p, hook) jacobs (p, f, hook);
+    else
+      __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
+      dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+    endif
   endif
   if (dfdp_pstruct)
     if (pnonscalar)
@@ -730,17 +749,18 @@
 	      ({s_diffp, s_diff_onesided, s_orig_lbound, \
 		s_orig_ubound, s_plabels, \
 		cell2fields(num2cell(hook.fixed), pord(nonfixed), \
-			    1, s_orig_fixed)}, \
+			    1, s_orig_fixed), cstep}, \
 	       {"diffp", "diff_onesided", "lbound", "ubound", \
-		"plabels", "fixed"}, \
+		"plabels", "fixed", "h"}, \
 	       2, hook));
   else
     dfdp = @ (p, hook) \
 	dfdp (p, cell2fields \
 	      ({diffp, diff_onesided, orig_lbound, orig_ubound, \
-		plabels, assign(orig_fixed, nonfixed, hook.fixed)}, \
+		plabels, assign(orig_fixed, nonfixed, hook.fixed), \
+		cstep}, \
 	       {"diffp", "diff_onesided", "lbound", "ubound", \
-		"plabels", "fixed"}, \
+		"plabels", "fixed", "h"}, \
 	       2, hook));
   endif
 
@@ -751,18 +771,18 @@
 		    ({s_diffp, s_diff_onesided, s_orig_lbound, \
 		      s_orig_ubound, s_plabels, \
 		      cell2fields(num2cell(hook.fixed), pord(nonfixed), \
-				  1, s_orig_fixed)}, \
+				  1, s_orig_fixed), cstep}, \
 		     {"diffp", "diff_onesided", "lbound", "ubound", \
-		      "plabels", "fixed"}, \
+		      "plabels", "fixed", "h"}, \
 		     2, hook));
   else
     df_gencstr = @ (p, func, idx, hook) \
 	df_gencstr (p, func, idx, cell2fields \
 		    ({diffp, diff_onesided, orig_lbound, \
 		      orig_ubound, plabels, \
-		      assign(orig_fixed, nonfixed, hook.fixed)}, \
+		      assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
 		     {"diffp", "diff_onesided", "lbound", "ubound", \
-		      "plabels", "fixed"}, \
+		      "plabels", "fixed", "h"}, \
 		     2, hook));
   endif
 
@@ -773,18 +793,18 @@
 		     ({s_diffp, s_diff_onesided, s_orig_lbound, \
 		       s_orig_ubound, s_plabels, \
 		       cell2fields(num2cell(hook.fixed), pord(nonfixed), \
-				   1, s_orig_fixed)}, \
+				   1, s_orig_fixed), cstep}, \
 		      {"diffp", "diff_onesided", "lbound", "ubound", \
-		       "plabels", "fixed"}, \
+		       "plabels", "fixed", "h"}, \
 		      2, hook));
   else
     df_genecstr = @ (p, func, idx, hook) \
 	df_genecstr (p, func, idx, cell2fields \
 		     ({diffp, diff_onesided, orig_lbound, \
 		       orig_ubound, plabels, \
-		       assign(orig_fixed, nonfixed, hook.fixed)}, \
+		       assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
 		      {"diffp", "diff_onesided", "lbound", "ubound", \
-		       "plabels", "fixed"}, \
+		       "plabels", "fixed", "h"}, \
 		      2, hook));
   endif
 
--- a/main/optim/inst/private/__residmin_stat__.m	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/inst/private/__residmin_stat__.m	Tue Nov 01 13:17:41 2011 +0000
@@ -44,6 +44,10 @@
     optimget = @ __optimget__;
   endif
 
+  ## scalar defaults
+  diffp_default = .001;
+  cstep_default = 1e-20;
+
   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
     ret = optimset ("param_config", [], \
 		    "param_order", [], \
@@ -53,6 +57,8 @@
 		    "dfdp", [], \
 		    "diffp", [], \
 		    "diff_onesided", [], \
+		    "complex_step_derivative", false, \
+		    "cstep", cstep_default, \
 		    "fixed", [], \
 		    "weights", [], \
 		    "residuals", [], \
@@ -65,9 +71,6 @@
     return;
   endif
 
-  ## scalar defaults
-  diffp_default = .001;
-
   assign = @ assign; # Is this faster in repeated calls?
 
   if (nargin != 4)
@@ -98,6 +101,11 @@
   diff_onesided = optimget (settings, "diff_onesided");
   fixed = optimget (settings, "fixed");
   residuals = optimget (settings, "residuals");
+  do_cstep = optimget (settings, "complex_step_derivative", false);
+  cstep = optimget (settings, "cstep", cstep_default);
+  if (do_cstep && ! isempty (dfdp))
+    error ("both 'complex_step_derivative' and 'dfdp' are set");
+  endif
 
   any_vector_conf = ! (isempty (diffp) && isempty (diff_onesided) && \
 		       isempty (fixed));
@@ -335,8 +343,12 @@
   ## jacobian of model function
   if (isempty (dfdp))
     if (! isempty (f))
-      __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
-      dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+      if (do_cstep)
+	dfdp = @ (p, hook) jacobs (p, f, hook);
+      else
+	__dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
+	dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+      endif
     endif
   elseif (! isa (dfdp, "function_handle"))
     if (ismatrix (dfdp))
@@ -484,15 +496,15 @@
   if (dfdp_pstruct)
     dfdp = @ (p, hook) \
 	dfdp (p, cell2fields \
-	      ({s_diffp, s_diff_onesided, s_plabels, s_fixed}, \
-	       {"diffp", "diff_onesided", "plabels", "fixed"}, \
+	      ({s_diffp, s_diff_onesided, s_plabels, s_fixed, cstep}, \
+	       {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
 	       2, hook));
   else
     if (! isempty (dfdp))
       dfdp = @ (p, hook) \
 	  dfdp (p, cell2fields \
-		({diffp, diff_onesided, plabels, fixed}, \
-		 {"diffp", "diff_onesided", "plabels", "fixed"}, \
+		({diffp, diff_onesided, plabels, fixed, cstep}, \
+		 {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
 		 2, hook));
     endif
   endif
--- a/main/optim/inst/residmin_stat.m	Tue Nov 01 10:45:54 2011 +0000
+++ b/main/optim/inst/residmin_stat.m	Tue Nov 01 13:17:41 2011 +0000
@@ -39,7 +39,7 @@
 ##
 ## @code{dfdp}: Jacobian of model function with respect to parameters.
 ##
-## @code{covr}: Covariance matrix of data (typically guessed by applying
+## @code{covd}: Covariance matrix of data (typically guessed by applying
 ## a factor to the covariance matrix of the residuals).
 ##
 ## @code{covp}: Covariance matrix of final parameters.
@@ -54,10 +54,10 @@
 ## The following settings have the same meaning as in
 ## @code{nonlin_residmin} (please refer to there): @code{param_order},
 ## @code{param_dims}, @code{f_pstruct}, @code{dfdp_pstruct},
-## @code{diffp}, @code{diff_onesided}, @code{fixed}, and @code{weights}.
-## Similarly, @code{param_config} can be used, but only with fields
-## corresponding to the settings @code{fixed}, @code{diffp}, and
-## @code{diff_onesided}.
+## @code{diffp}, @code{diff_onesided}, @code{complex_step_derivative},
+## @code{cstep}, @code{fixed}, and @code{weights}. Similarly,
+## @code{param_config} can be used, but only with fields corresponding
+## to the settings @code{fixed}, @code{diffp}, and @code{diff_onesided}.
 ##
 ## @code{dfdp} can be set in the same way as in @code{nonlin_residmin},
 ## but alternatively may already contain the computed Jacobian of the