Menu
▾
▴
octave-cvsupdate
|
From: <i7...@us...> - 2011-05-20 09:52:46
|
Revision: 8274
http://octave.svn.sourceforge.net/octave/?rev=8274&view=rev
Author: i7tiol
Date: 2011-05-20 09:52:37 +0000 (Fri, 20 May 2011)
Log Message:
-----------
Replace function calls: cell2cell -> num2cell, partarray -> mat2cell.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/inst/__dfdp__.m
trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2011-05-20 07:45:05 UTC (rev 8273)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2011-05-20 09:52:37 UTC (rev 8274)
@@ -1,11 +1,11 @@
Name: Optim
Version: 1.0.17
-Date: 2011-03-07
+Date: 2011-05-20
Author: Various Authors
Maintainer: The Octave Community
Title: Optimzation.
Description: Non-linear optimization toolkit.
-Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.11), struct (>= 1.0.9)
+Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.9)
Autoload: yes
License: GPL version 2 or later and GFDL
Url: http://octave.sf.net
Modified: trunk/octave-forge/main/optim/inst/__dfdp__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/__dfdp__.m 2011-05-20 07:45:05 UTC (rev 8273)
+++ trunk/octave-forge/main/optim/inst/__dfdp__.m 2011-05-20 09:52:37 UTC (rev 8274)
@@ -66,7 +66,7 @@
if (isfield (hook, 'plabels'))
plabels = hook.plabels;
else
- plabels = cell2cell (num2cell ((1:n).'), 1);
+ plabels = num2cell (num2cell ((1:n).'));
end
else
@@ -75,7 +75,7 @@
diffp = .001 * ones (n, 1);
lbound = - Inf (n, 1);
ubound = Inf (n, 1);
- plabels = cell2cell (num2cell ((1:n).'), 1);
+ plabels = num2cell (num2cell ((1:n).'));
end
prt = zeros (m, n); % initialise Jacobian to Zero
Modified: trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m 2011-05-20 07:45:05 UTC (rev 8273)
+++ trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m 2011-05-20 09:52:37 UTC (rev 8274)
@@ -178,8 +178,7 @@
error ("lengths of 'param_order' and 'param_dims' not equal");
endif
pnel = cellfun (@ prod, pdims);
- ppartidx = cumsum (pnel);
- ppartidx = [[1; ppartidx(1:end-1)+1], ppartidx];
+ ppartidx = pnel;
if (any (pnel > 1))
pnonscalar = true;
cpnel = num2cell (pnel);
@@ -210,10 +209,10 @@
endif
endif
- plabels = cell2cell (num2cell ((1:np).'), 1);
+ plabels = num2cell (num2cell ((1:np).'));
if (! isempty (pord))
- plabels = cat (2, plabels, cell2cell (epord, 1), \
- cell2cell (num2cell (psubidx), 1));
+ plabels = cat (2, plabels, num2cell (epord), \
+ num2cell (num2cell (psubidx)));
endif
## some useful vectors
@@ -422,7 +421,7 @@
if (pnonscalar)
f = @ (p, varargin) \
f (cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), \
pord, 1), varargin{:});
else
@@ -449,7 +448,7 @@
cat (2, \
fields2cell \
(dfdp (cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), \
pord, 1), hook), \
pord){:});
@@ -467,7 +466,7 @@
if (pnonscalar)
f_genicstr = @ (p, varargin) \
f_genicstr (cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), \
pord, 1), varargin{:});
else
@@ -488,7 +487,7 @@
fields2cell \
(df_gencstr \
(cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), pord, 1), \
func, idx, hook), \
pord){:});
@@ -507,7 +506,7 @@
if (pnonscalar)
f_genecstr = @ (p, varargin) \
f_genecstr (cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), \
pord, 1), varargin{:});
else
@@ -528,7 +527,7 @@
fields2cell \
(df_genecstr \
(cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), pord, 1), \
func, idx, hook), \
pord){:});
@@ -564,40 +563,39 @@
emc(idx(prepidx), :) = cat (1, fields2cell (semc, pord(idx)){:});
endif
- ## parameter-related configuration for jacobi functions
+ ## parameter-related configuration for jacobian functions
if (dfdp_pstruct || df_inequc_pstruct || df_equc_pstruct)
if(pnonscalar)
s_diffp = cell2struct \
- (cellfun (@ reshape, partarray (diffp, ppartidx), \
+ (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
pdims, "UniformOutput", false), pord, 1);
s_diff_onesided = cell2struct \
- (cellfun (@ reshape, partarray (diff_onesided, ppartidx), \
+ (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
pdims, "UniformOutput", false), pord, 1);
s_orig_lbound = cell2struct \
- (cellfun (@ reshape, partarray (lbound, ppartidx), \
+ (cellfun (@ reshape, mat2cell (lbound, ppartidx), \
pdims, "UniformOutput", false), pord, 1);
s_orig_ubound = cell2struct \
- (cellfun (@ reshape, partarray (ubound, ppartidx), \
+ (cellfun (@ reshape, mat2cell (ubound, ppartidx), \
pdims, "UniformOutput", false), pord, 1);
s_plabels = cell2struct \
- (cell2cell \
+ (num2cell \
(cat (2, cellfun \
(@ (x) cellfun \
- (@ reshape, partarray (cat (1, x{:}), ppartidx), \
+ (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
pdims, "UniformOutput", false), \
- cell2cell (plabels, 2), "UniformOutput", false){:}), \
- 1), \
+ num2cell (plabels, 1), "UniformOutput", false){:}), \
+ 2), \
pord, 1);
- s_plabels = cell2struct (cell2cell (plabels, 1), pord, 1);
s_orig_fixed = cell2struct \
- (cellfun (@ reshape, partarray (fixed, ppartidx), \
+ (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
pdims, "UniformOutput", false), pord, 1);
else
s_diffp = cell2struct (num2cell (diffp), pord, 1);
s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
s_orig_lbound = cell2struct (num2cell (lbound), pord, 1);
s_orig_ubound = cell2struct (num2cell (ubound), pord, 1);
- s_plabels = cell2struct (cell2cell (plabels, 1), pord, 1);
+ s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
s_orig_fixed = cell2struct (num2cell (fixed), pord, 1);
endif
endif
@@ -913,7 +911,7 @@
if (pin_struct)
if (pnonscalar)
p = cell2struct \
- (cellfun (@ reshape, partarray (p, ppartidx), \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
pdims, "UniformOutput", false), \
pord, 1);
else
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2011-05-20 07:45:05 UTC (rev 8273)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2011-05-20 09:52:37 UTC (rev 8274)
@@ -174,9 +174,10 @@
##
## The setting @code{param_order} is a cell-array with names of the
## optimized parameters. If not given, and initial parameters are a
-## structure, all parameters in the structure are optimized. It is an
-## error if @code{param_order} is not given and there are any
-## non-structure-based configuration items or functions.
+## structure, all parameters in the structure are optimized. If initial
+## parameters are a structure, it is an error if @code{param_order} is
+## not given and there are any non-structure-based configuration items
+## or functions.
##
## The initial parameters @var{pin} can be given as a structure
## containing at least all fields named in @code{param_order}. In this
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2011-07-04 13:02:11
|
Revision: 8374
http://octave.svn.sourceforge.net/octave/?rev=8374&view=rev
Author: i7tiol
Date: 2011-07-04 13:02:03 +0000 (Mon, 04 Jul 2011)
Log Message:
-----------
Add general statistics framework for residual-based minimization.
Modified Paths:
--------------
trunk/octave-forge/main/optim/PKG_ADD
trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m
trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
trunk/octave-forge/main/optim/src/Makefile
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/__residmin_stat__.m
trunk/octave-forge/main/optim/inst/curvefit_stat.m
trunk/octave-forge/main/optim/inst/private/__covd_wls__.m
trunk/octave-forge/main/optim/inst/private/__covp_corp_wls__.m
trunk/octave-forge/main/optim/inst/private/__null_optim__.m
trunk/octave-forge/main/optim/inst/residmin_stat.m
trunk/octave-forge/main/optim/src/__disna_optim__.cc
Modified: trunk/octave-forge/main/optim/PKG_ADD
===================================================================
--- trunk/octave-forge/main/optim/PKG_ADD 2011-07-03 17:55:56 UTC (rev 8373)
+++ trunk/octave-forge/main/optim/PKG_ADD 2011-07-04 13:02:03 UTC (rev 8374)
@@ -2,4 +2,5 @@
## 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__");
endif
Modified: trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m 2011-07-03 17:55:56 UTC (rev 8373)
+++ trunk/octave-forge/main/optim/inst/__nonlin_residmin__.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -129,10 +129,10 @@
if (! user_df_genecstr) df_equc_pstruct = false; endif
## some settings require a parameter order
- if ((pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
- f_equc_pstruct || f_pstruct || dfdp_pstruct || \
- df_inequc_pstruct || df_equc_pstruct || mc_struct || \
- emc_struct))
+ if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
+ f_equc_pstruct || f_pstruct || dfdp_pstruct || \
+ df_inequc_pstruct || df_equc_pstruct || mc_struct || \
+ emc_struct)
if (isempty (pord))
if (pin_struct)
if (any_vector_conf || \
Added: trunk/octave-forge/main/optim/inst/__residmin_stat__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/__residmin_stat__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/__residmin_stat__.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,608 @@
+## Copyright (C) 2011 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## Internal function, called by residmin_stat --- see there --- and
+## others. Calling __residmin_stat__ indirectly hides the argument
+## "hook", usable by wrappers, from users. Currently, hook can contain
+## the field "observations". Since much uf the interface code is taken
+## from __nonlin_residmin__, it may be that not everything is ideal for
+## the present case; but I think it's allright to leave it so.
+##
+## Some general considerations while making this function:
+##
+## Different Functions for optimization statistics should be made for
+## mere objective function optimization (to be made yet) and
+## residual-derived minimization (this function), since there are
+## different computing aspects. Don't put the contained functionality
+## (statistics) into the respective optimization functions (or
+## backends), since different optimization algorithms can share a way to
+## compute statistics (e.g. even stochastic optimizers can mimize
+## (weighted) squares of residuals). Also, don't use the same frontend
+## for optimization and statistics, since the differences in the
+## interface for both uses may be confusing otherwise, also the optimset
+## options only partially overlap.
+
+## disabled PKG_ADD: __all_opts__ ("__residmin_stat__");
+
+function ret = __residmin_stat__ (f, pfin, settings, hook)
+
+ 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
+ optimget = @ __optimget__;
+ endif
+
+ if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
+ ret = optimset ("param_config", [], \
+ "param_order", [], \
+ "param_dims", [], \
+ "f_pstruct", false, \
+ "dfdp_pstruct", false, \
+ "dfdp", [], \
+ "diffp", [], \
+ "diff_onesided", [], \
+ "fixed", [], \
+ "weights", [], \
+ "residuals", [], \
+ "covd", [], \
+ "objf", [], \ # no default, e.g. "wls"
+ "ret_dfdp", false, \
+ "ret_covd", false, \
+ "ret_covp", false, \
+ "ret_corp", false);
+ return;
+ endif
+
+ ## scalar defaults
+ diffp_default = .001;
+
+ assign = @ assign; # Is this faster in repeated calls?
+
+ if (nargin != 4)
+ error ("incorrect number of arguments");
+ endif
+
+ if (ischar (f))
+ f = str2func (f);
+ endif
+
+ if (! (p_struct = isstruct (pfin)))
+ if (! isempty (pfin) && (! isvector (pfin) || columns (pfin) > 1))
+ error ("parameters must be either a structure or a column vector");
+ endif
+ endif
+
+ #### processing of settings and consistency checks
+
+ pconf = optimget (settings, "param_config");
+ pord = optimget (settings, "param_order");
+ pdims = optimget (settings, "param_dims");
+ f_pstruct = optimget (settings, "f_pstruct", false);
+ dfdp_pstruct = optimget (settings, "dfdp_pstruct", f_pstruct);
+ dfdp = optimget (settings, "dfdp");
+ if (ischar (dfdp)) dfdp = str2func (dfdp); endif
+ if (isstruct (dfdp)) dfdp_pstruct = true; endif
+ diffp = optimget (settings, "diffp");
+ diff_onesided = optimget (settings, "diff_onesided");
+ fixed = optimget (settings, "fixed");
+ residuals = optimget (settings, "residuals");
+
+ any_vector_conf = ! (isempty (diffp) && isempty (diff_onesided) && \
+ isempty (fixed));
+
+ ## correct "_pstruct" settings if functions are not supplied
+ if (isempty (dfdp)) dfdp_pstruct = false; endif
+ if (isempty (f)) f_pstruct = false; endif
+
+ ## some settings require a parameter order
+ if (p_struct || ! isempty (pconf) || f_pstruct || dfdp_pstruct)
+ if (isempty (pord))
+ if (p_struct)
+ if (any_vector_conf || \
+ ! ((f_pstruct || isempty (f)) && \
+ (dfdp_pstruct || isempty (dfdp))))
+ error ("no parameter order specified and constructing a parameter order from the structure of parameters can not be done since not all configuration or given functions are structure based");
+ else
+ pord = fieldnames (pfin);
+ endif
+ else
+ error ("given settings require specification of parameter order or parameters in the form of a structure");
+ endif
+ endif
+ pord = pord(:);
+ if (p_struct && ! all (arefields (pfin, pord)))
+ error ("some parameters lacking");
+ endif
+ if ((nnames = rows (unique (pord))) < rows (pord))
+ error ("duplicate parameter names in 'param_order'");
+ endif
+ if (isempty (pdims))
+ if (p_struct)
+ pdims = cellfun \
+ (@ size, fields2cell (pfin, pord), "UniformOutput", false);
+ else
+ pdims = num2cell (ones (nnames, 2), 2);
+ endif
+ else
+ pdims = pdims(:);
+ if (p_struct && \
+ ! all (cellfun (@ (x, y) prod (size (x)) == prod (y), \
+ struct2cell (pfin), pdims)))
+ error ("given param_dims and dimensions of parameters do not match");
+ endif
+ endif
+ if (nnames != rows (pdims))
+ error ("lengths of 'param_order' and 'param_dims' not equal");
+ endif
+ pnel = cellfun (@ prod, pdims);
+ ppartidx = pnel;
+ if (any (pnel > 1))
+ pnonscalar = true;
+ cpnel = num2cell (pnel);
+ prepidx = cat (1, cellfun \
+ (@ (x, n) x(ones (1, n), 1), \
+ num2cell ((1:nnames).'), cpnel, \
+ "UniformOutput", false){:});
+ epord = pord(prepidx, 1);
+ psubidx = cat (1, cellfun \
+ (@ (n) (1:n).', cpnel, \
+ "UniformOutput", false){:});
+ else
+ pnonscalar = false; # some less expensive interfaces later
+ prepidx = (1:nnames).';
+ epord = pord;
+ psubidx = ones (nnames, 1);
+ endif
+ else
+ pord = []; # spares checks for given but not needed
+ endif
+
+ if (p_struct)
+ np = sum (pnel);
+ else
+ np = length (pfin);
+ if (! isempty (pord) && np != sum (pnel))
+ error ("number of initial parameters not correct");
+ endif
+ endif
+ if (ismatrix (dfdp) && ! ischar (dfdp) && ! isempty (dfdp) && \
+ np == 0)
+ np = columns (dfdp);
+ endif
+
+ plabels = num2cell (num2cell ((1:np).'));
+ if (! isempty (pord))
+ plabels = cat (2, plabels, num2cell (epord), \
+ num2cell (num2cell (psubidx)));
+ endif
+
+ ## some useful vectors
+ zerosvec = zeros (np, 1);
+ falsevec = false (np, 1);
+ sizevec = [np, 1];
+
+ ## collect parameter-related configuration
+ if (! isempty (pconf))
+ ## use supplied configuration structure
+
+ ## parameter-related configuration is either allowed by a structure
+ ## or by vectors
+ if (any_vector_conf)
+ error ("if param_config is given, its potential items must not \
+ be configured in another way");
+ endif
+
+ ## supplement parameter names lacking in param_config
+ nidx = ! arefields (pconf, pord);
+ pconf = cell2fields ({struct()}(ones (1, sum (nidx))), \
+ pord(nidx), 2, pconf);
+
+ pconf = structcat (1, fields2cell (pconf, pord){:});
+
+ ## in the following, use reshape with explicit dimensions (instead
+ ## of x(:)) so that errors are thrown if a configuration item has
+ ## incorrect number of elements
+
+ diffp = zerosvec;
+ diffp(:) = diffp_default;
+ if (isfield (pconf, "diffp"))
+ idx = ! fieldempty (pconf, "diffp");
+ if (pnonscalar)
+ diffp(idx(prepidx)) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).diffp}.', cpnel(idx), \
+ "UniformOutput", false){:});
+ else
+ diffp(idx) = [pconf.diffp];
+ endif
+ endif
+
+ diff_onesided = fixed = falsevec;
+
+ if (isfield (pconf, "diff_onesided"))
+ idx = ! fieldempty (pconf, "diff_onesided");
+ if (pnonscalar)
+ diff_onesided(idx(prepidx)) = \
+ logical \
+ (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).diff_onesided}.', cpnel(idx), \
+ "UniformOutput", false){:}));
+ else
+ diff_onesided(idx) = logical ([pconf.diff_onesided]);
+ endif
+ endif
+
+ if (isfield (pconf, "fixed"))
+ idx = ! fieldempty (pconf, "fixed");
+ if (pnonscalar)
+ fixed(idx(prepidx)) = \
+ logical \
+ (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).fixed}.', cpnel(idx), \
+ "UniformOutput", false){:}));
+ else
+ fixed(idx) = logical ([pconf.fixed]);
+ endif
+ endif
+ else
+ ## use supplied configuration vectors
+
+ if (isempty (diffp))
+ diffp = zerosvec;
+ diffp(:) = diffp_default;
+ else
+ if (any (size (diffp) != sizevec))
+ error ("diffp: wrong dimensions");
+ endif
+ diffp(isna (diffp)) = diffp_default;
+ endif
+
+ if (isempty (diff_onesided))
+ diff_onesided = falsevec;
+ else
+ if (any (size (diff_onesided) != sizevec))
+ error ("diff_onesided: wrong dimensions")
+ endif
+ diff_onesided(isna (diff_onesided)) = false;
+ diff_onesided = logical (diff_onesided);
+ endif
+
+ if (isempty (fixed))
+ fixed = falsevec;
+ else
+ if (any (size (fixed) != sizevec))
+ error ("fixed: wrong dimensions");
+ endif
+ fixed(isna (fixed)) = false;
+ fixed = logical (fixed);
+ endif
+ endif
+
+ #### consider whether parameters and functions are based on parameter
+ #### structures or parameter vectors; wrappers for call to default
+ #### function for jacobians
+
+ ## parameters
+ if (p_struct)
+ if (pnonscalar)
+ pfin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ fields2cell (pfin, pord), cpnel, \
+ "UniformOutput", false){:});
+ else
+ pfin = cat (1, fields2cell (pfin, pord){:});
+ endif
+ endif
+
+ ## model function
+ if (f_pstruct)
+ if (pnonscalar)
+ f = @ (p, varargin) \
+ f (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), varargin{:});
+ else
+ f = @ (p, varargin) \
+ f (cell2struct (num2cell (p), pord, 1), varargin{:});
+ endif
+ endif
+ if (isempty (residuals))
+ if (isempty (f))
+ error ("neither model function nor residuals given");
+ endif
+ residuals = f (pfin);
+ endif
+ if (isfield (hook, "observations"))
+ if (any (size (residuals) != size (obs = hook.observations)))
+ error ("dimensions of observations and values of model function must match");
+ endif
+ f = @ (p) f (p) - obs;
+ residuals -= obs;
+ endif
+
+ ## jacobian of model function
+ if (isempty (dfdp))
+ if (! isempty (f))
+ dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+ endif
+ elseif (! isa (dfdp, "function_handle"))
+ if (ismatrix (dfdp))
+ if (any (size (dfdp) != [prod(size(residuals)), np]))
+ error ("jacobian has wrong size");
+ endif
+ elseif (! dfdp_pstruct)
+ error ("jacobian has wrong type");
+ endif
+ dfdp = @ (varargin) dfdp; # simply make a function returning it
+ endif
+ if (dfdp_pstruct)
+ if (pnonscalar)
+ dfdp = @ (p, hook) \
+ cat (2, \
+ fields2cell \
+ (dfdp (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), hook), \
+ pord){:});
+ else
+ dfdp = @ (p, hook) \
+ cat (2, \
+ fields2cell \
+ (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
+ pord){:});
+ endif
+ endif
+
+ ## parameter-related configuration for jacobian function
+ if (dfdp_pstruct)
+ if(pnonscalar)
+ s_diffp = cell2struct \
+ (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_diff_onesided = cell2struct \
+ (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_plabels = cell2struct \
+ (num2cell \
+ (cat (2, cellfun \
+ (@ (x) cellfun \
+ (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
+ pdims, "UniformOutput", false), \
+ num2cell (plabels, 1), "UniformOutput", false){:}), \
+ 2), \
+ pord, 1);
+ s_orig_fixed = cell2struct \
+ (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ else
+ s_diffp = cell2struct (num2cell (diffp), pord, 1);
+ s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
+ s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
+ s_fixed = cell2struct (num2cell (fixed), pord, 1);
+ endif
+ endif
+
+ #### further values and checks
+
+ ## check weights dimensions
+ weights = optimget (settings, "weights", ones (size (residuals)));
+ if (any (size (weights) != size (residuals)))
+ error ("dimension of weights and residuals must match");
+ endif
+
+
+ #### collect remaining settings
+ need_dfdp = false;
+ covd = optimget (settings, "covd");
+ need_objf_label = false;
+ if ((ret_dfdp = optimget (settings, "ret_dfdp", false)))
+ need_dfdp = true;
+ endif
+ if ((ret_covd = optimget (settings, "ret_covd", false)))
+ need_objf_label = true;
+ if (np == 0)
+ error ("number of parameters must be known for 'covd', specify either parameters or a jacobian matrix");
+ endif
+ endif
+ if ((ret_covp = optimget (settings, "ret_covp", false)))
+ need_objf_label = true;
+ need_dfdp = true;
+ endif
+ if ((ret_corp = optimget (settings, "ret_corp", false)))
+ need_objf_label = true;
+ need_dfdp = true;
+ endif
+ if (need_objf_label)
+ if (isempty (objf = optimget (settings, "objf")))
+ error ("label of objective function must be specified");
+ else
+ funs = map_objf (objf);
+ endif
+ else
+ funs = struct ();
+ endif
+ if (isempty (dfdp) && need_dfdp)
+ error ("jacobian required and default function for jacobian requires a model function");
+ endif
+
+ ####
+
+ ## Everything which is computed is stored in a hook structure which is
+ ## passed to and returned by every backend function. This hook is not
+ ## identical to the returned structure, since some more results could
+ ## be computed by the way.
+
+ #### handle fixing of parameters
+
+ orig_p = pfin;
+
+ if (all (fixed) && ! isempty (fixed))
+ error ("no free parameters");
+ endif
+
+ ## The policy should be that everything which is computed is left as
+ ## it is up to the end --- since other computations might need it in
+ ## this form --- and supplemented with values corresponding to fixed
+ ## parameters (mostly NA, probably) not until then.
+
+ nonfixed = ! fixed;
+
+ np_after_fixing = sum (nonfixed);
+
+ if (any (fixed))
+
+ if (! isempty (pfin))
+ pfin = pfin(nonfixed);
+ endif
+
+ ## model function
+ f = @ (p, varargin) f (assign (pfin, nonfixed, p), varargin{:});
+
+ ## jacobian of model function
+ if (! isempty (dfdp))
+ dfdp = @ (p, hook) \
+ dfdp (assign (pfin, nonfixed, p), hook)(:, nonfixed);
+ endif
+
+ endif
+
+ #### supplement constants to jacobian function
+ if (dfdp_pstruct)
+ dfdp = @ (p, hook) \
+ dfdp (p, cell2fields \
+ ({s_diffp, s_diff_onesided, s_plabels, s_fixed}, \
+ {"diffp", "diff_onesided", "plabels", "fixed"}, \
+ 2, hook));
+ else
+ if (! isempty (dfdp))
+ dfdp = @ (p, hook) \
+ dfdp (p, cell2fields \
+ ({diffp, diff_onesided, plabels, fixed}, \
+ {"diffp", "diff_onesided", "plabels", "fixed"}, \
+ 2, hook));
+ endif
+ endif
+
+ #### prepare interface hook
+
+ ## passed final parameters of an optimization
+ hook.pfin = pfin;
+
+ ## passed function for derivative of model function
+ hook.dfdp = dfdp;
+
+ ## passed function for complementary pivoting
+ ## hook.cpiv = cpiv; # set before
+
+ ## passed value of residual function for initial parameters
+ hook.residuals = residuals;
+
+ ## passed weights
+ hook.weights = weights;
+
+ ## passed dimensions
+ hook.np = np_after_fixing;
+ hook.nm = prod (size (residuals));
+
+ ## passed statistics functions
+ hook.funs = funs;
+
+ ## passed covariance matrix of data (if given by user)
+ if (! isempty (covd))
+ covd_dims = size (covd);
+ if (length (covd_dims) != 2 || any (covd_dims != hook.nm))
+ error ("wrong dimensions of covariance matrix of data");
+ endif
+ hook.covd = covd;
+ endif
+
+ #### do the actual work
+
+ if (ret_dfdp)
+ hook.jac = hook.dfdp (hook.pfin, hook);
+ endif
+
+ if (ret_covd)
+ hook = funs.covd (hook);
+ endif
+
+ if (ret_covp || ret_corp)
+ hook = funs.covp_corp (hook);
+ endif
+
+ #### convert (consider fixing ...) and return results
+
+ ret = struct ();
+
+ if (ret_dfdp)
+ ret.dfdp = zeros (hook.nm, np);
+ ret.dfdp(:, nonfixed) = hook.jac;
+ endif
+
+ if (ret_covd)
+ ret.covd = hook.covd;
+ endif
+
+ if (ret_covp)
+ if (any (fixed))
+ ret.covp = NA (np);
+ ret.covp(nonfixed, nonfixed) = hook.covp;
+ else
+ ret.covp = hook.covp;
+ endif
+ endif
+
+ if (ret_corp)
+ if (any (fixed))
+ ret.corp = NA (np);
+ ret.corp(nonfixed, nonfixed) = hook.corp;
+ else
+ ret.corp = hook.corp;
+ endif
+ endif
+
+endfunction
+
+function funs = map_objf (objf)
+
+ switch (objf)
+ case "wls" # weighted least squares
+ funs.covd = str2func ("__covd_wls__");
+ funs.covp_corp = str2func ("__covp_corp_wls__");
+ otherwise
+ error ("no statistics implemented for objective function '%s'", \
+ objf);
+ endswitch
+
+endfunction
+
+function lval = assign (lval, lidx, rval)
+
+ lval(lidx) = rval;
+
+endfunction
+
+function ret = __optimget__ (s, name, default)
+
+ if (isfield (s, name))
+ ret = s.(name);
+ elseif (nargin > 2)
+ ret = default;
+ else
+ ret = [];
+ endif
+
+endfunction
Added: trunk/octave-forge/main/optim/inst/curvefit_stat.m
===================================================================
--- trunk/octave-forge/main/optim/inst/curvefit_stat.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/curvefit_stat.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,78 @@
+## Copyright (C) 2011 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{info} =} residmin_stat (@var{f}, @var{p}, @var{x}, @var{y}, @var{settings})
+##
+## Frontend for computation of statistics for fitting of values,
+## computed by a model function, to observed values.
+##
+## Please refer to the description of @code{residmin_stat}. The only
+## differences to @code{residmin_stat} are the additional arguments
+## @var{x} (independent values) and @var{y} (observations), that the
+## model function @var{f}, if provided, has a second obligatory argument
+## which will be set to @var{x} and is supposed to return guesses for
+## the observations (with the same dimensions), and that the possibly
+## user-supplied function for the jacobian of the model function has
+## also a second obligatory argument which will be set to @var{x}.
+##
+## @seealso {residmin_stat}
+##
+## @end deftypefn
+
+function ret = curvefit_stat (f, pfin, x, y, settings)
+
+ if (nargin == 1)
+ ret = __residmin_stat__ (f);
+ return;
+ endif
+
+ if (nargin != 5)
+ print_usage ()
+ endif
+
+ 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
+ optimget = @ __optimget__;
+ endif
+ if (! isempty (dfdp = optimget (settings, "dfdp")) && \
+ ! (ismatrix (dfdp) && ! ischar (dfdp)))
+ if (ischar (dfdp))
+ dfdp = str2func (dfdp);
+ endif
+ settings.dfdp = @ (p, varargin) dfdp (p, x, varargin{:});
+ endif
+ if (! isempty (f))
+ f = @ (p) f (p, x);
+ endif
+
+ ret = __residmin_stat__ \
+ (f, pfin, settings, struct ("observations", y));
+
+endfunction
+
+function ret = __optimget__ (s, name, default)
+
+ if (isfield (s, name))
+ ret = s.(name);
+ elseif (nargin > 2)
+ ret = default;
+ else
+ ret = [];
+ endif
+
+endfunction
Modified: trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2011-07-03 17:55:56 UTC (rev 8373)
+++ trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -53,8 +53,11 @@
settings = struct ();
endif
- ## optimset mechanism is broken in Octave 3.2.4
- optimget = @ __optimget__;
+ 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
+ optimget = @ __optimget__;
+ endif
if (! isempty (dfdp = optimget (settings, "dfdp")))
if (ischar (dfdp))
dfdp = str2func (dfdp);
Added: trunk/octave-forge/main/optim/inst/private/__covd_wls__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__covd_wls__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__covd_wls__.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,32 @@
+## Copyright (C) 2011 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+function hook = __covd_wls__ (hook)
+
+ m = hook.nm;
+ n = hook.np;
+
+ if (m <= n)
+ error ("guessing covariance-matrix of residuals for weighted least squares requires at least one more residual than free parameters");
+ endif
+
+ w = hook.weights(:);
+ res = hook.residuals(:);
+
+ w2 = w .^ 2;
+
+ hook.covd = diag (res.' * diag (w2) * res / (m - n) ./ w2);
+
+endfunction
Added: trunk/octave-forge/main/optim/inst/private/__covp_corp_wls__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__covp_corp_wls__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__covp_corp_wls__.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,105 @@
+## Copyright (C) 2011 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## This is based on Bard, Nonlinear Parameter Estimation, Academic
+## Press, 1974, section 7-5, but also on own re-calculations for the
+## specific case of weighted least squares. The part with only certain
+## elements of covp or corp being defined is self-made, I don't know a
+## reference.
+
+function hook = __covp_corp_wls__ (hook)
+
+ ## compute jacobian, if not already done
+ if (! isfield (hook, "jac"))
+ hook.jac = hook.dfdp (hook.pfin, hook);
+ endif
+ jact = (jac = hook.jac).';
+
+ ## compute guessed covariance matrix of data, if not already done
+ if (! isfield (hook, "covd"))
+ hook = hook.funs.covd (hook);
+ endif
+ covd_inv = inv (covd = hook.covd);
+
+ if (rcond (A = jact * covd_inv * jac) > eps)
+
+ covp = hook.covp = inv (A);
+
+ d = sqrt (diag (covp));
+
+ hook.corp = covp ./ (d * d.');
+
+ else
+
+ n = hook.np;
+
+ covp = NA (n);
+
+ ## Now we have the equation "A * covp * A.' == A".
+
+ ## Find a particular solution for "covp * A.'".
+ part_covp_At = A \ A;
+
+ ## Find a particular solution for "covp". Only uniquely defined
+ ## elements (identified later) will be further used.
+ part_covp = A \ part_covp_At.';
+
+ ## Find a basis for the nullspace of A.
+ if (true) # test for Octave version once submitted patch is applied
+ # to Octave (bug #33503)
+ null = @ __null_optim__;
+ endif
+ if (isempty (basis = null (A)))
+ error ("internal error, singularity assumed, but null-space computed to be zero-dimensional");
+ endif
+
+ ## Find an index (applied to both row and column) of uniquely
+ ## defined elements of covp.
+ idun = all (basis == 0, 2);
+
+ ## Fill in these elements.
+ covp(idun, idun) = part_covp(idun, idun);
+
+ ## Compute corp as far as possible at the moment.
+ d = sqrt (diag (covp));
+ corp = covp ./ (d * d.');
+
+ ## All diagonal elements of corp should be one, even those as yet
+ ## NA.
+ corp(1 : n + 1 : n * n) = 1;
+
+ ## If there are indices, applied to both row and column, so that
+ ## indexed elements within one row or one column are determined up
+ ## to a multiple of a vector, find both these vectors and the
+ ## respective indices. In the same run, use them to further fill in
+ ## corp as described below.
+ for id = 1 : (cb = columns (basis))
+ if (any (idx = \
+ all (basis(:, (1 : cb) != id) == 0, 2) & \
+ basis(:, id) != 0))
+ vec = sign (basis(idx, id));
+
+ ## Depending on "vec", single coefficients of correlation
+ ## indexed by "idx" are either +1 or -1.
+ corp(idx, idx) = vec * vec.';
+ endif
+ endfor
+
+ hook.covp = covp;
+ hook.corp = corp;
+
+ endif
+
+endfunction
Added: trunk/octave-forge/main/optim/inst/private/__null_optim__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__null_optim__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__null_optim__.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,120 @@
+## Copyright (C) 1994-2011 John W. Eaton
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} null (@var{A})
+## @deftypefnx {Function File} {} null (@var{A}, @var{tol})
+## Return an orthonormal basis of the null space of @var{A}.
+##
+## The dimension of the null space is taken as the number of singular
+## values of @var{A} not greater than @var{tol}. If the argument @var{tol}
+## is missing, it is computed as
+##
+## @example
+## max (size (@var{A})) * max (svd (@var{A})) * eps
+## @end example
+## @seealso{orth}
+## @end deftypefn
+
+## Author: KH <Kur...@wu...>
+## Created: 24 December 1993.
+## Adapted-By: jwe
+## Adapted-By: Olaf Till <ola...@un...>
+
+## This function has also been submitted to Octave (bug #33503).
+
+function retval = __null_optim__ (A, tol)
+
+ if (isempty (A))
+ retval = [];
+ else
+ [U, S, V] = svd (A);
+
+ [rows, cols] = size (A);
+
+ [S_nr, S_nc] = size (S);
+
+ if (S_nr == 1 || S_nc == 1)
+ s = S(1);
+ else
+ s = diag (S);
+ endif
+
+ if (nargin == 1)
+ if (isa (A, "single"))
+ tol = max (size (A)) * (vtol = s (1) * (meps = eps ("single")));
+ else
+ tol = max (size (A)) * (vtol = s (1) * (meps = eps));
+ endif
+ elseif (nargin != 2)
+ print_usage ();
+ endif
+
+ rank = sum (s > tol);
+
+ if (rank < cols)
+ retval = V (:, rank+1:cols);
+
+ if (rows >= cols)
+ cb = columns (retval);
+
+ ## Set those elements of each vector to zero whose absolute
+ ## values are smallest and which together could be zero without
+ ## making the angle to the originally computed vector larger
+ ## than given by the error bound. Do this in an approximative
+ ## but numerically feasible way.
+
+ ## error bounds of basis vectors in radians, see LAPACK user
+ ## guide, http://www.netlib.org/lapack/lug/node96.html
+ if (true) # test for Octave version once submitted patch is applied
+ # to Octave (bug #33503)
+ __disna__ = @ __disna_optim__;
+ endif
+ ebnd = vtol ./ (__disna__ ("R", s, rows, cols)(rank+1:cols));
+
+ ## sort elements by magnitude
+ sb = conj (retval) .* retval;
+ [sb, idx] = sort (sb);
+ idx += repmat (0:cols:cols*(cb-1), cols, 1); # for un-sorting
+
+ ## norms of vectors made by all elements up to this
+ sb = sqrt (cumsum (sb));
+
+ ## The norm of the vectors made up by elements settable to zero
+ ## is small enough to be approximately equal to the angle
+ ## between the full vectors before and after setting these
+ ## elements to zero (considering the norms of the full vectors
+ ## being 1). Index of approximated angles not exceeding error
+ ## bound.
+ zidx = sb <= repmat (ebnd, cols, 1);
+
+ ## set indexed elements to zero in original basis
+ zidx = zidx(idx);
+ retval(zidx) = 0;
+
+ else
+ ## no error bounds computable with LAPACK
+
+ retval(abs (retval) < meps) = 0;
+ endif
+ else
+ retval = zeros (cols, 0);
+ endif
+ endif
+
+endfunction
Added: trunk/octave-forge/main/optim/inst/residmin_stat.m
===================================================================
--- trunk/octave-forge/main/optim/inst/residmin_stat.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/residmin_stat.m 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,103 @@
+## Copyright (C) 2011 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{info} =} residmin_stat (@var{f}, @var{p}, @var{settings})
+##
+## Frontend for computation of statistics for a residual-based
+## minimization.
+##
+## @var{settings} is a structure whose fields can be set by
+## @code{optimset} with Octave versions 3.3.55 or greater; with older
+## Octave versions, the fields must be set directly and in the correct
+## case. With @var{settings} the computation of certain statistics is
+## requested by setting the fields @code{ret_<name_of_statistic>} to
+## @code{true}. The respective statistics will be returned in a
+## structure as fields with name @code{<name_of_statistic>}. Depending
+## on the requested statistic and on the additional information provided
+## in @var{settings}, @var{f} and @var{p} may be empty. Otherwise,
+## @var{f} is the model function of an optimization (the interface of
+## @var{f} is described e.g. in @code{nonlin_residmin}, please see
+## there), and @var{p} is a real column vector with parameters resulting
+## from the same optimization.
+##
+## Currently, the following statistics (or general information) can be
+## requested:
+##
+## @code{dfdp}: Jacobian of model function with respect to parameters.
+##
+## @code{covr}: Covariance matrix of data (typically guessed by applying
+## a factor to the covariance matrix of the residuals).
+##
+## @code{covp}: Covariance matrix of final parameters.
+##
+## @code{corp}: Correlation matrix of final parameters.
+##
+## Further @var{settings}
+##
+## The functionality of the interface is similar to
+## @code{nonlin_residmin}. In particular, structure-based, possibly
+## non-scalar, parameters and flagging parameters as fixed are possible.
+## 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{dfdp} can be set in the same way as in @code{nonlin_residmin},
+## but alternatively may already contain the computed Jacobian of the
+## model function at the final parameters in matrix- or structure-form.
+## Users may pass information on the result of the optimization in
+## @code{residuals} (self-explaining) and @code{covd} (covariance matrix
+## of data). In many cases the type of objective function of the
+## optimization must be specified in @code{objf}; currently, there is
+## only a backend for the type "wls" (weighted least squares).
+##
+## Backend-specific information
+##
+## The backend for @code{objf == "wls"} (currently the only backend)
+## computes @code{cord} (due to user request or as a prerequisite for
+## @code{covp} and @code{corp}) as a diagonal matrix by assuming that
+## the variances of data points are proportional to the reciprocal of
+## the squared @code{weights} and guessing the factor of proportionality
+## from the residuals. If @code{covp} is not defined (e.g. because the
+## Jacobian has no full rank), it makes an attempt to still compute its
+## uniquely defined elements, if any, and to find the additional defined
+## elements (being @code{1} or @code{-1}), if any, in @code{corp}.
+##
+## @seealso {curvefit_stat}
+##
+## @end deftypefn
+
+
+function ret = residmin_stat (varargin)
+
+ if (nargin == 1)
+ ret = __residmin_stat__ (varargin{1});
+ return;
+ endif
+
+ if (nargin != 3)
+ print_usage ();
+ endif
+
+ varargin{4} = struct ();
+
+ ret = __residmin_stat__ (varargin{:});
+
+endfunction
\ No newline at end of file
Modified: trunk/octave-forge/main/optim/src/Makefile
===================================================================
--- trunk/octave-forge/main/optim/src/Makefile 2011-07-03 17:55:56 UTC (rev 8373)
+++ trunk/octave-forge/main/optim/src/Makefile 2011-07-04 13:02:03 UTC (rev 8374)
@@ -2,7 +2,7 @@
MKOCTFILE := mkoctfile
endif
-all: __bfgsmin.oct numgradient.oct numhessian.oct samin.oct
+all: __bfgsmin.oct numgradient.oct numhessian.oct samin.oct __disna_optim__.oct
%.oct: %.cc
$(MKOCTFILE) -s $<
Added: trunk/octave-forge/main/optim/src/__disna_optim__.cc
===================================================================
--- trunk/octave-forge/main/optim/src/__disna_optim__.cc (rev 0)
+++ trunk/octave-forge/main/optim/src/__disna_optim__.cc 2011-07-04 13:02:03 UTC (rev 8374)
@@ -0,0 +1,165 @@
+/*
+
+Copyright (C) 2011 Olaf Till
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING. If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+// This function has also been submitted to Octave (bug #33503).
+
+#include <octave/oct.h>
+#include "f77-fcn.h"
+
+extern "C"
+{
+ F77_RET_T
+ F77_FUNC (ddisna, DDISNA) (F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&,
+ const octave_idx_type&,
+ const double*,
+ double*,
+ octave_idx_type&);
+
+ F77_RET_T
+ F77_FUNC (sdisna, SDISNA) (F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&,
+ const octave_idx_type&,
+ const float*,
+ float*,
+ octave_idx_type&);
+}
+
+DEFUN_DLD (__disna_optim__, args, ,
+ "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{rcond} =} __disna__ (@var{job}, @var{d})\n\
+@deftypefnx {Loadable Function} {@var{rcond} =} __disna__ (@var{job}, @var{d}, @var{m}, @var{n})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+ /*
+ Interface to DDISNA and SDISNA of LAPACK.
+
+ If job is 'E', no third or fourth argument are given. If job is 'L'
+ or 'R', M and N are given.
+ */
+
+ std::string fname ("__disna__");
+
+ octave_value retval;
+
+ if (args.length () != 2 && args.length () != 4)
+ print_usage ();
+
+ std::string job_str (args(0).string_value ());
+
+ char job;
+
+ if (job_str.length () != 1)
+ error ("%s: invalid job label", fname.c_str ());
+ else
+ job = job_str[0];
+
+ octave_idx_type m, n, l;
+ bool single;
+ octave_value d;
+
+ if (args(1).is_single_type ())
+ {
+ single = true;
+ d = args(1).float_column_vector_value ();
+ }
+ else
+ {
+ single = false;
+ d = args(1).column_vector_value ();
+ }
+
+ if (! error_state)
+ {
+ l = d.length ();
+ switch (job)
+ {
+ case 'E' :
+ if (args.length () != 2)
+ error ("%s: with job label 'E' only two arguments are allowed",
+ fname.c_str ());
+ else
+ m = l;
+ break;
+ case 'L' :
+ case 'R' :
+ if (args.length () != 4)
+ error ("%s: with job labels 'L' or 'R', four arguments must be given",
+ fname.c_str ());
+ else
+ {
+ m = args(2).idx_type_value ();
+ n = args(3).idx_type_value ();
+ if (! error_state)
+ {
+ octave_idx_type md = m < n ? m : n;
+ if (l != md)
+ error ("%s: given dimensions don't match length of second argument",
+ fname.c_str ());
+ }
+ }
+ break;
+ default :
+ error ("%s: job label not correct", fname.c_str ());
+ }
+ }
+
+ if (error_state)
+ {
+ error ("%s: invalid arguments", fname.c_str ());
+ return retval;
+ }
+
+ octave_idx_type info;
+
+ if (single)
+ {
+ FloatColumnVector srcond (l);
+
+ F77_XFCN (sdisna, SDISNA, (F77_CONST_CHAR_ARG2 (&job, 1),
+ m, n,
+ d.float_column_vector_value ().fortran_vec (),
+ srcond.fortran_vec (),
+ info));
+
+ retval = srcond;
+ }
+ else
+ {
+ ColumnVector drcond (l);
+
+ F77_XFCN (ddisna, DDISNA, (F77_CONST_CHAR_ARG2 (&job, 1),
+ m, n,
+ d.column_vector_value ().fortran_vec (),
+ drcond.fortran_vec (),
+ info));
+
+ retval = drcond;
+ }
+
+ if (info < 0)
+ error ("%s: LAPACK routine says %i-th argument had an illegal value",
+ fname.c_str (), -info);
+
+ return retval;
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2011-11-01 13:17:51
|
Revision: 8906
http://octave.svn.sourceforge.net/octave/?rev=8906&view=rev
Author: i7tiol
Date: 2011-11-01 13:17:41 +0000 (Tue, 01 Nov 2011)
Log Message:
-----------
Use Fotios' complex step derivative approximation. And a few small fixes.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/PKG_ADD
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/private/__collect_constraints__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
trunk/octave-forge/main/optim/inst/private/__residmin_stat__.m
trunk/octave-forge/main/optim/inst/residmin_stat.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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.
Modified: trunk/octave-forge/main/optim/PKG_ADD
===================================================================
--- trunk/octave-forge/main/optim/PKG_ADD 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/PKG_ADD 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
Modified: trunk/octave-forge/main/optim/inst/private/__collect_constraints__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__collect_constraints__.m 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/inst/private/__collect_constraints__.m 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
Modified: trunk/octave-forge/main/optim/inst/private/__residmin_stat__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__residmin_stat__.m 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/inst/private/__residmin_stat__.m 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
Modified: trunk/octave-forge/main/optim/inst/residmin_stat.m
===================================================================
--- trunk/octave-forge/main/optim/inst/residmin_stat.m 2011-11-01 10:45:54 UTC (rev 8905)
+++ trunk/octave-forge/main/optim/inst/residmin_stat.m 2011-11-01 13:17:41 UTC (rev 8906)
@@ -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
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2011-11-01 14:09:10
|
Revision: 8907
http://octave.svn.sourceforge.net/octave/?rev=8907&view=rev
Author: i7tiol
Date: 2011-11-01 14:08:59 +0000 (Tue, 01 Nov 2011)
Log Message:
-----------
Fix INDEX.
Modified Paths:
--------------
trunk/octave-forge/main/optim/INDEX
trunk/octave-forge/main/optim/inst/test_fminunc_1.m
Modified: trunk/octave-forge/main/optim/INDEX
===================================================================
--- trunk/octave-forge/main/optim/INDEX 2011-11-01 13:17:41 UTC (rev 8906)
+++ trunk/octave-forge/main/optim/INDEX 2011-11-01 14:08:59 UTC (rev 8907)
@@ -18,16 +18,18 @@
nonlin_residmin
nonlin_curvefit
LinearRegression
+Optimization statistics
+ residmin_stat
+ curvefit_stat
Zero finding
vfzero
Compatibility
fminunc_compat
- optimset_compat
quadprog= try Yinyu Ye's <a href="http://dollar.biz.uiowa.edu/col/ye/matlab.html">code</a>
linprog
Numerical derivatives
dfdp dcdp dfpdp dfxpdp cdiff deriv
- numgradient numhessian
+ numgradient numhessian jacobs
Pivoting
cpiv_bard
gjp
Modified: trunk/octave-forge/main/optim/inst/test_fminunc_1.m
===================================================================
--- trunk/octave-forge/main/optim/inst/test_fminunc_1.m 2011-11-01 13:17:41 UTC (rev 8906)
+++ trunk/octave-forge/main/optim/inst/test_fminunc_1.m 2011-11-01 14:08:59 UTC (rev 8907)
@@ -10,7 +10,7 @@
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
## for more details.
-## test_fminunc_compat_1 - Test that fminunc_compat and optimset_compat work
+## test_fminunc_compat_1 - Test that fminunc_compat works
##
## A quadratic function is fminunc_compatd. Various options are tested. Options
## are passed incomplete (to see if properly completed) and
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2011-11-01 14:11:49
|
Revision: 8908
http://octave.svn.sourceforge.net/octave/?rev=8908&view=rev
Author: i7tiol
Date: 2011-11-01 14:11:37 +0000 (Tue, 01 Nov 2011)
Log Message:
-----------
Revert a mistake in previous commit.
Modified Paths:
--------------
trunk/octave-forge/main/optim/INDEX
trunk/octave-forge/main/optim/inst/test_fminunc_1.m
Modified: trunk/octave-forge/main/optim/INDEX
===================================================================
--- trunk/octave-forge/main/optim/INDEX 2011-11-01 14:08:59 UTC (rev 8907)
+++ trunk/octave-forge/main/optim/INDEX 2011-11-01 14:11:37 UTC (rev 8908)
@@ -25,6 +25,7 @@
vfzero
Compatibility
fminunc_compat
+ optimset_compat
quadprog= try Yinyu Ye's <a href="http://dollar.biz.uiowa.edu/col/ye/matlab.html">code</a>
linprog
Numerical derivatives
Modified: trunk/octave-forge/main/optim/inst/test_fminunc_1.m
===================================================================
--- trunk/octave-forge/main/optim/inst/test_fminunc_1.m 2011-11-01 14:08:59 UTC (rev 8907)
+++ trunk/octave-forge/main/optim/inst/test_fminunc_1.m 2011-11-01 14:11:37 UTC (rev 8908)
@@ -10,7 +10,7 @@
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
## for more details.
-## test_fminunc_compat_1 - Test that fminunc_compat works
+## test_fminunc_compat_1 - Test that fminunc_compat and optimset_compat work
##
## A quadratic function is fminunc_compatd. Various options are tested. Options
## are passed incomplete (to see if properly completed) and
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2011-11-30 08:36:30
|
Revision: 9223
http://octave.svn.sourceforge.net/octave/?rev=9223&view=rev
Author: i7tiol
Date: 2011-11-30 08:36:23 +0000 (Wed, 30 Nov 2011)
Log Message:
-----------
Nir Krakauer added function powell.m.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/INDEX
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/powell.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2011-11-29 17:48:56 UTC (rev 9222)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2011-11-30 08:36:23 UTC (rev 9223)
@@ -1,6 +1,6 @@
Name: Optim
Version: 1.0.18
-Date: 2011-11-01
+Date: 2011-11-30
Author: Various Authors
Maintainer: The Octave Community
Title: Optimzation.
Modified: trunk/octave-forge/main/optim/INDEX
===================================================================
--- trunk/octave-forge/main/optim/INDEX 2011-11-29 17:48:56 UTC (rev 9222)
+++ trunk/octave-forge/main/optim/INDEX 2011-11-30 08:36:23 UTC (rev 9223)
@@ -6,6 +6,7 @@
nrm
fmin
line_min
+ powell
fmins adsmax mdsmax nmsmax
bfgsmin samin battery
fminsearch
Added: trunk/octave-forge/main/optim/inst/powell.m
===================================================================
--- trunk/octave-forge/main/optim/inst/powell.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/powell.m 2011-11-30 08:36:23 UTC (rev 9223)
@@ -0,0 +1,196 @@
+## Copyright (C) 2011 Nir Krakauer
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} [@var{p}, @var{obj_value}, @var{convergence}, @var{iters}, @var{nevs}] = powell (@var{f}, @var{p0}, @var{control})
+##powell: implements a direction-set (Powell's) method for multidimensional minimization of a function without calculation of the gradient [1, 2]
+##
+## @subheading Arguments
+##
+## @itemize @bullet
+## @item
+## @var{f}: name of function to minimize (string or handle), which should accept one input variable (see example for how to pass on additional input arguments)
+##
+## @item
+## @var{p0}: An initial value of the function argument to minimize
+##
+## @item
+## @var{options}: an optional structure, which can be generated by optimset, with some or all of the following fields:
+## @itemize @minus
+## @item
+## MaxIter: maximum iterations (positive integer, or -1 or Inf for unlimited (default))
+## @item
+## TolFun: minimum amount by which function value must decrease in each iteration to continue (default is 1E-8)
+## @item
+## MaxFunEvals: maximum function evaluations (positive integer, or -1 or Inf for unlimited (default))
+## @item
+## SearchDirections: an n*n matrix whose columns contain the initial set of (presumably orthogonal) directions to minimize along, where n is the number of elements in the argument to be minimized for; or an n*1 vector of magnitudes for the initial directions (defaults to the set of unit direction vectors)
+## @end itemize
+## @end itemize
+##
+## @subheading Examples
+##
+## @example
+## @group
+## y = @@(x, s) x(1) ^ 2 + x(2) ^ 2 + s;
+## o = optimset('MaxIter', 100, 'TolFun', 1E-10);
+## s = 1;
+## [x_optim, y_min, conv, iters, nevs] = powell(@@(x) y(x, s), [1 0.5], o); %pass y wrapped in an anonymous function so that all other arguments to y, which are held constant, are set
+## %should return something like x_optim = [4E-14 3E-14], y_min = 1, conv = 1, iters = 2, nevs = 24
+## @end group
+##
+## @end example
+##
+## @subheading Returns:
+##
+## @itemize @bullet
+## @item
+## @var{p}: the minimizing value of the function argument
+## @item
+## @var{obj_value}: the value of @var{f}() at @var{p}
+## @item
+## @var{convergence}: 1 if normal convergence, 0 if not
+## @item
+## @var{iters}: number of iterations performed
+## @item
+## @var{nevs}: number of function evaluations
+## @end itemize
+##
+## @subheading References
+##
+## @enumerate
+## @item
+## Powell MJD (1964), An efficient method for finding the minimum of a function of several variables without calculating derivatives, @cite{Computer Journal}, 7 :155-162
+##
+## @item
+## Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (1992). @cite{Numerical Recipes in Fortran: The Art of Scientific Computing} (2nd Ed.). New York: Cambridge University Press (Section 10.5)
+## @end enumerate
+## @end deftypefn
+
+## Author: Nir Krakauer <nkr...@cc...>
+## Description: Multidimensional minimization (direction-set method)
+
+## PKG_ADD: __all_opts__ ("powell");
+
+function [p, obj_value, convergence, iters, nevs] = powell (f, p0, options);
+
+ if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
+ p = optimset ("MaxIter", Inf, \
+ "TolFun", 1e-8, \
+ "MaxFunEvals", Inf, \
+ "SearchDirections", []);
+ return;
+ endif
+
+ ## check number of arguments
+ if ((nargin < 2) || (nargin > 3))
+ usage('powell: you must supply 2 or 3 arguments');
+ endif
+
+
+ ## default or input values
+
+ if (nargin < 3)
+ options = struct ();
+ endif
+
+ xi_set = 0;
+ xi = optimget (options, 'SearchDirections');
+ if (! isempty (xi))
+ if (isvector (xi)) # assume that xi is is n*1 or 1*n
+ xi = diag (x);
+ endif
+ xi_set = 1;
+ endif
+
+
+ MaxIter = optimget (options, 'MaxIter', Inf);
+ if (MaxIter < 0) MaxIter = Inf; endif
+ MaxFunEvals = optimget (options, 'MaxFunEvals', Inf);
+ TolFun = optimget (options, 'TolFun', 1E-8);
+
+ nevs = 0;
+ iters = 0;
+ convergence = 0;
+
+ p = p0; # initial value of the argument being minimized
+
+ try
+ obj_value = f(p);
+ catch
+ error ("function does not exist or cannot be evaluated");
+ end_try_catch
+
+ nevs++;
+
+ n = numel (p); # number of dimensions to minimize over
+
+ xit = zeros (n, 1);
+ if (! xi_set)
+ xi = eye(n);
+ endif
+
+
+
+ ## do an iteration
+ while (iters <= MaxIter && nevs <= MaxFunEvals && ! convergence)
+ iters++;
+ pt = p; # best point as iteration begins
+ fp = obj_value; # value of the objective function as iteration begins
+ ibig = 0; # will hold direction along which the objective function decreased the most in this iteration
+ dlt = 0; # will hold decrease in objective function value in this iteration
+ for i = 1:n
+ xit = reshape (xi(:, i), size(p));
+ fptt = obj_value;
+ [a, obj_value, nev] = line_min (f, xit, {p});
+ nevs = nevs + nev;
+ p = p + a*xit;
+ change = fptt - obj_value;
+ if (change > dlt)
+ dlt = change;
+ ibig = i;
+ endif
+ endfor
+
+ if ( 2*abs(fp-obj_value) <= TolFun*(abs(fp) + abs(obj_value)) )
+ convergence = 1;
+ return
+ endif
+
+ if (iters == MaxIter)
+ disp ("iteration maximum exceeded");
+ return
+ endif
+
+ ## attempt parabolic extrapolation
+ ptt = 2*p - pt;
+ xit = p - pt;
+ fptt = f(ptt);
+ nevs++;
+ if (fptt < fp) # check whether the extrapolation actually makes the objective function smaller
+ t = 2 * (fp - 2*obj_value + fptt) * (fp-obj_value-dlt)^2 - dlt * (fp-fptt)^2;
+ if (t < 0)
+ p = ptt;
+ [a, obj_value, nev] = line_min (f, xit, {p});
+ nevs = nevs + nev;
+ p = p + a*xit;
+
+ ## add the net direction from this iteration to the direction set
+ xi(:, ibig) = xi(:, n);
+ xi(:, n) = xit(:);
+ endif
+ endif
+ endwhile
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2011-12-29 23:10:32
|
Revision: 9485
http://octave.svn.sourceforge.net/octave/?rev=9485&view=rev
Author: paramaniac
Date: 2011-12-29 23:10:26 +0000 (Thu, 29 Dec 2011)
Log Message:
-----------
optim: apply patch from Marco Atzeri, slightly modified to improve compatibility with older versions of make
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/src/Makefile
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2011-12-29 22:49:10 UTC (rev 9484)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2011-12-29 23:10:26 UTC (rev 9485)
@@ -3,7 +3,7 @@
Date: 2011-11-30
Author: Various Authors
Maintainer: The Octave Community
-Title: Optimzation.
+Title: Optimization.
Description: Non-linear optimization toolkit.
Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.9)
Autoload: yes
Modified: trunk/octave-forge/main/optim/src/Makefile
===================================================================
--- trunk/octave-forge/main/optim/src/Makefile 2011-12-29 22:49:10 UTC (rev 9484)
+++ trunk/octave-forge/main/optim/src/Makefile 2011-12-29 23:10:26 UTC (rev 9485)
@@ -2,10 +2,12 @@
MKOCTFILE := mkoctfile
endif
+LAPACK_LIBS = $$(mkoctfile -p LAPACK_LIBS)
+
all: __bfgsmin.oct numgradient.oct numhessian.oct samin.oct __disna_optim__.oct
%.oct: %.cc
- $(MKOCTFILE) -s $<
+ $(MKOCTFILE) -s $< ${LAPACK_LIBS}
clean:
-rm *.o core octave-core *.oct *~
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2012-03-30 12:29:23
|
Revision: 10113
http://octave.svn.sourceforge.net/octave/?rev=10113&view=rev
Author: carandraug
Date: 2012-03-30 12:29:17 +0000 (Fri, 30 Mar 2012)
Log Message:
-----------
optim: package is no longer automatically loaded
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/NEWS
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-03-30 12:27:14 UTC (rev 10112)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-03-30 12:29:17 UTC (rev 10113)
@@ -1,11 +1,11 @@
Name: Optim
Version: 1.0.18
Date: 2011-11-30
-Author: Various Authors
-Maintainer: The Octave Community
+Author: various authors
+Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Description: Non-linear optimization toolkit.
Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.9)
-Autoload: yes
+Autoload: no
License: GPL version 2 or later and GFDL
Url: http://octave.sf.net
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-03-30 12:27:14 UTC (rev 10112)
+++ trunk/octave-forge/main/optim/NEWS 2012-03-30 12:29:17 UTC (rev 10113)
@@ -17,3 +17,5 @@
deriv linprog
** The function `line_min' has a configurable setpesize and max evals
+
+ ** Package is no longer automatically loaded.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-04-18 07:15:18
|
Revision: 10274
http://octave.svn.sourceforge.net/octave/?rev=10274&view=rev
Author: jpicarbajal
Date: 2012-04-18 07:15:07 +0000 (Wed, 18 Apr 2012)
Log Message:
-----------
optim: Numerical derivatives using the Cauchy integral.
Modified Paths:
--------------
trunk/octave-forge/main/optim/INDEX
trunk/octave-forge/main/optim/NEWS
Added Paths:
-----------
trunk/octave-forge/main/optim/cauchy.m
Modified: trunk/octave-forge/main/optim/INDEX
===================================================================
--- trunk/octave-forge/main/optim/INDEX 2012-04-17 15:54:44 UTC (rev 10273)
+++ trunk/octave-forge/main/optim/INDEX 2012-04-18 07:15:07 UTC (rev 10274)
@@ -1,6 +1,6 @@
optimization >> Optimization
Minimization
- minimize
+ minimize
nelder_mead_min
d2_min
nrm
@@ -31,7 +31,7 @@
linprog
Numerical derivatives
dfdp dcdp dfpdp dfxpdp cdiff deriv
- numgradient numhessian jacobs
+ numgradient numhessian jacobs cauchy
Pivoting
cpiv_bard
gjp
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-04-17 15:54:44 UTC (rev 10273)
+++ trunk/octave-forge/main/optim/NEWS 2012-04-18 07:15:07 UTC (rev 10274)
@@ -3,7 +3,7 @@
** The following functions are new optim 1.1.0:
- powell
+ powell cauchy
** The following functions have been deprecated since they have been
implemented in Octave core:
Added: trunk/octave-forge/main/optim/cauchy.m
===================================================================
--- trunk/octave-forge/main/optim/cauchy.m (rev 0)
+++ trunk/octave-forge/main/optim/cauchy.m 2012-04-18 07:15:07 UTC (rev 10274)
@@ -0,0 +1,123 @@
+## Copyright (C) 2011 Fernando Damian Nieuwveldt <fdn...@gm...>
+## 2012 Adapted by Juan Pablo Carbajal <car...@if...>
+##
+## This program is free software; you can redistribute it and/or
+## modify it under the terms of the GNU General Public License
+## as published by the Free Software Foundation; either version 3
+## of the License, or (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} cauchy (@var{N}, @var{r}, @var{x}, @var{f} )
+## Return the Taylor coefficients and numerical differentiation of a function
+## @var{f} for the first @var{N-1} coefficients or derivatives using the fft.
+## @var{N} is the number of points to evaluate,
+## @var{r} is the radius of convergence, needs to be chosen less then the smallest singularity,
+## @var{x} is point to evaluate the Taylor expansion or differentiation. For example,
+##
+## If @var{x} is a scalar, the function @var{f} is evaluated in a row vector
+## of length @var{N}. If @var{x} is a column vector, @var{f} is evaluated in a
+## matrix of length(x)-by-N elements and must return a matrix of the same size.
+##
+## @example
+## @group
+## d = cauchy(16,1.5,0,@(x) exp(x));
+## @result{} d(2) = 1.0000 # first (2-1) derivative of function f (index starts from zero)
+## @end group
+## @end example
+## @end deftypefn
+
+function deriv = cauchy(N, r, x, f)
+
+ if nargin != 4
+ print_usage ();
+ end
+
+ [nx m] = size (x);
+ if m > 1
+ error('cauchy:InvalidArgument', 'The 3rd argument must be a column vector');
+ end
+
+ n = 0:N-1;
+ th = 2*pi*n/N;
+
+ f_p = f (bsxfun (@plus, x, r * exp (i * th) ) );
+
+ evalfft = real(fft (f_p, [], 2));
+
+ deriv = bsxfun (@times, evalfft, 1./(N*(r.^n)).* factorial(n)) ;
+
+endfunction
+
+%!demo
+%! # Cauchy integral formula: Application to Hermite polynomials
+%! # Author: Fernando Damian Nieuwveldt
+%! # Edited by: Juan Pablo Carbajal
+%!
+%! function g = hermite(order,x)
+%!
+%! ## N should be bigger than order+1
+%! N = 32;
+%! r = 0.5;
+%! Hnx = @(t) exp ( bsxfun (@minus, kron(t(:).', x(:)) , t(:).'.^2/2) );
+%! Hnxfft = cauchy(N, r, 0, Hnx);
+%! g = Hnxfft(:, order+1);
+%!
+%! endfunction
+%!
+%! t = linspace(-1,1,30);
+%! he2 = hermite(2,t);
+%! he2_ = t.^2-1;
+%!
+%! plot(t,he2,'bo;Contour integral representation;', t,he2_,'r;Exact;');
+%! grid
+%! clear all
+%!
+%! % --------------------------------------------------------------------------
+%! % The plots compares the approximation of the Hermite polynomial using the
+%! % Cauchy integral (circles) and the corresposind polynomial H_2(x) = x.^2 - 1.
+%! % See http://en.wikipedia.org/wiki/Hermite_polynomials#Contour_integral_representation
+
+%!demo
+%! # Cauchy integral formula: Application to Hermite polynomials
+%! # Author: Fernando Damian Nieuwveldt
+%! # Edited by: Juan Pablo Carbajal
+%!
+%! xx = sort (rand (100,1));
+%! yy = sin (3*2*pi*xx);
+%!
+%! # Exact first derivative derivative
+%! diffy = 6*pi*cos (3*2*pi*xx);
+%!
+%! np = [10 15 30 100];
+%!
+%! for i =1:4
+%! idx = sort(randperm (100,np(i)));
+%! x = xx(idx);
+%! y = yy(idx);
+%!
+%! p = spline (x,y);
+%! yval = ppval (ppder(p),x);
+%! # Use the cauchy formula for computing the derivatives
+%! deriv = cauchy (fix (np(i)/4), .1, x, @(x) sin (3*2*pi*x));
+%!
+%! subplot(2,2,i)
+%! h = plot(xx,diffy,'-b;Exact;',...
+%! x,yval,'-or;ppder solution;',...
+%! x,deriv(:,2),'-og;Cauchy formula;');
+%! set(h(1),'linewidth',2);
+%! set(h(2:3),'markersize',3);
+%!
+%! legend('Location','Northoutside','Orientation','horizontal');
+%!
+%! if i!=1
+%! legend('hide','Orientation','horizontal')
+%! end
+%! end
+%!
+%! % --------------------------------------------------------------------------
+%! % The plots compares the derivatives calculated with Cauchy and with ppder.
+%! % Each subplot shows the results with increasing number of samples.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-04-19 09:00:58
|
Revision: 10290
http://octave.svn.sourceforge.net/octave/?rev=10290&view=rev
Author: jpicarbajal
Date: 2012-04-19 09:00:48 +0000 (Thu, 19 Apr 2012)
Log Message:
-----------
optim: Numerical derivatives using the Cauchy integral. Was in worng place
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
Removed Paths:
-------------
trunk/octave-forge/main/optim/cauchy.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-04-19 08:28:21 UTC (rev 10289)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-04-19 09:00:48 UTC (rev 10290)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.0.18
-Date: 2011-11-30
+Version: 1.1.0
+Date: 2012-XX-XX
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Deleted: trunk/octave-forge/main/optim/cauchy.m
===================================================================
--- trunk/octave-forge/main/optim/cauchy.m 2012-04-19 08:28:21 UTC (rev 10289)
+++ trunk/octave-forge/main/optim/cauchy.m 2012-04-19 09:00:48 UTC (rev 10290)
@@ -1,123 +0,0 @@
-## Copyright (C) 2011 Fernando Damian Nieuwveldt <fdn...@gm...>
-## 2012 Adapted by Juan Pablo Carbajal <car...@if...>
-##
-## This program is free software; you can redistribute it and/or
-## modify it under the terms of the GNU General Public License
-## as published by the Free Software Foundation; either version 3
-## of the License, or (at your option) any later version.
-##
-## This program is distributed in the hope that it will be useful,
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {} cauchy (@var{N}, @var{r}, @var{x}, @var{f} )
-## Return the Taylor coefficients and numerical differentiation of a function
-## @var{f} for the first @var{N-1} coefficients or derivatives using the fft.
-## @var{N} is the number of points to evaluate,
-## @var{r} is the radius of convergence, needs to be chosen less then the smallest singularity,
-## @var{x} is point to evaluate the Taylor expansion or differentiation. For example,
-##
-## If @var{x} is a scalar, the function @var{f} is evaluated in a row vector
-## of length @var{N}. If @var{x} is a column vector, @var{f} is evaluated in a
-## matrix of length(x)-by-N elements and must return a matrix of the same size.
-##
-## @example
-## @group
-## d = cauchy(16,1.5,0,@(x) exp(x));
-## @result{} d(2) = 1.0000 # first (2-1) derivative of function f (index starts from zero)
-## @end group
-## @end example
-## @end deftypefn
-
-function deriv = cauchy(N, r, x, f)
-
- if nargin != 4
- print_usage ();
- end
-
- [nx m] = size (x);
- if m > 1
- error('cauchy:InvalidArgument', 'The 3rd argument must be a column vector');
- end
-
- n = 0:N-1;
- th = 2*pi*n/N;
-
- f_p = f (bsxfun (@plus, x, r * exp (i * th) ) );
-
- evalfft = real(fft (f_p, [], 2));
-
- deriv = bsxfun (@times, evalfft, 1./(N*(r.^n)).* factorial(n)) ;
-
-endfunction
-
-%!demo
-%! # Cauchy integral formula: Application to Hermite polynomials
-%! # Author: Fernando Damian Nieuwveldt
-%! # Edited by: Juan Pablo Carbajal
-%!
-%! function g = hermite(order,x)
-%!
-%! ## N should be bigger than order+1
-%! N = 32;
-%! r = 0.5;
-%! Hnx = @(t) exp ( bsxfun (@minus, kron(t(:).', x(:)) , t(:).'.^2/2) );
-%! Hnxfft = cauchy(N, r, 0, Hnx);
-%! g = Hnxfft(:, order+1);
-%!
-%! endfunction
-%!
-%! t = linspace(-1,1,30);
-%! he2 = hermite(2,t);
-%! he2_ = t.^2-1;
-%!
-%! plot(t,he2,'bo;Contour integral representation;', t,he2_,'r;Exact;');
-%! grid
-%! clear all
-%!
-%! % --------------------------------------------------------------------------
-%! % The plots compares the approximation of the Hermite polynomial using the
-%! % Cauchy integral (circles) and the corresposind polynomial H_2(x) = x.^2 - 1.
-%! % See http://en.wikipedia.org/wiki/Hermite_polynomials#Contour_integral_representation
-
-%!demo
-%! # Cauchy integral formula: Application to Hermite polynomials
-%! # Author: Fernando Damian Nieuwveldt
-%! # Edited by: Juan Pablo Carbajal
-%!
-%! xx = sort (rand (100,1));
-%! yy = sin (3*2*pi*xx);
-%!
-%! # Exact first derivative derivative
-%! diffy = 6*pi*cos (3*2*pi*xx);
-%!
-%! np = [10 15 30 100];
-%!
-%! for i =1:4
-%! idx = sort(randperm (100,np(i)));
-%! x = xx(idx);
-%! y = yy(idx);
-%!
-%! p = spline (x,y);
-%! yval = ppval (ppder(p),x);
-%! # Use the cauchy formula for computing the derivatives
-%! deriv = cauchy (fix (np(i)/4), .1, x, @(x) sin (3*2*pi*x));
-%!
-%! subplot(2,2,i)
-%! h = plot(xx,diffy,'-b;Exact;',...
-%! x,yval,'-or;ppder solution;',...
-%! x,deriv(:,2),'-og;Cauchy formula;');
-%! set(h(1),'linewidth',2);
-%! set(h(2:3),'markersize',3);
-%!
-%! legend('Location','Northoutside','Orientation','horizontal');
-%!
-%! if i!=1
-%! legend('hide','Orientation','horizontal')
-%! end
-%! end
-%!
-%! % --------------------------------------------------------------------------
-%! % The plots compares the derivatives calculated with Cauchy and with ppder.
-%! % Each subplot shows the results with increasing number of samples.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-05-24 05:50:50
|
Revision: 10510
http://octave.svn.sourceforge.net/octave/?rev=10510&view=rev
Author: i7tiol
Date: 2012-05-24 05:50:43 +0000 (Thu, 24 May 2012)
Log Message:
-----------
Untabify and some cleanup.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/PKG_ADD
trunk/octave-forge/main/optim/doc/development/interfaces.txt
trunk/octave-forge/main/optim/inst/cpiv_bard.m
trunk/octave-forge/main/optim/inst/dcdp.m
trunk/octave-forge/main/optim/inst/dfdp.m
trunk/octave-forge/main/optim/inst/dfpdp.m
trunk/octave-forge/main/optim/inst/dfxpdp.m
trunk/octave-forge/main/optim/inst/gjp.m
trunk/octave-forge/main/optim/inst/leasqr.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/optim_problems.m
trunk/octave-forge/main/optim/inst/private/__dfdp__.m
trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,11 +1,11 @@
Name: Optim
Version: 1.1.0
-Date: 2012-XX-XX
+Date: 2012-05-24
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Description: Non-linear optimization toolkit.
-Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.9)
+Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.10)
Autoload: no
License: GFDL, GPLv3+, modified BSD, public domain
Url: http://octave.sf.net
Modified: trunk/octave-forge/main/optim/PKG_ADD
===================================================================
--- trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 05:50:43 UTC (rev 10510)
@@ -3,4 +3,5 @@
## sometime in 3.3.54+, if I remember right
__all_opts__ ("nonlin_residmin");
__all_opts__ ("residmin_stat");
+ ## __all_opts__ ("nonlin_min");
endif
Modified: trunk/octave-forge/main/optim/doc/development/interfaces.txt
===================================================================
--- trunk/octave-forge/main/optim/doc/development/interfaces.txt 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/doc/development/interfaces.txt 2012-05-24 05:50:43 UTC (rev 10510)
@@ -80,6 +80,9 @@
'vc' (required): vector (possibly empty) of the function "mc.' *
parameters + vc" of linear constraints,
+If bounds have been specified, they are contained in 'mc' and 'vc'
+_before_ all other linear constraints.
+
'n_gencstr' (required): number of general constraints (except the
linear constraints given by 'mc' and 'vc',
Modified: trunk/octave-forge/main/optim/inst/cpiv_bard.m
===================================================================
--- trunk/octave-forge/main/optim/inst/cpiv_bard.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/cpiv_bard.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
Modified: trunk/octave-forge/main/optim/inst/dcdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dcdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dcdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
@@ -40,4 +40,3 @@
hook.diff_onesided = dp < 0;
prt = __dfdp__ (p, func, hook);
-end
Modified: trunk/octave-forge/main/optim/inst/dfdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,7 +1,7 @@
%% Copyright (C) 1992-1994 Richard Shrager
%% Copyright (C) 1992-1994 Arthur Jutan
%% Copyright (C) 1992-1994 Ray Muzic
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
Modified: trunk/octave-forge/main/optim/inst/dfpdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfpdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfpdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -48,5 +48,3 @@
end
ret = __dfdp__ (varargin{:});
-
-end
Modified: trunk/octave-forge/main/optim/inst/dfxpdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfxpdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfxpdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
@@ -51,5 +51,3 @@
end
ret = __dfdp__ (varargin{2:end});
-
-end
Modified: trunk/octave-forge/main/optim/inst/gjp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/gjp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/gjp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
@@ -50,4 +50,3 @@
m([1:k-1, k+1:end], l) * m(k, [1:l-1, l+1:end]);
m([1:k-1, k+1:end], l) = - m([1:k-1, k+1:end], l) / p; % pivot column
m(k, l) = 1 / p;
-end
Modified: trunk/octave-forge/main/optim/inst/leasqr.m
===================================================================
--- trunk/octave-forge/main/optim/inst/leasqr.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/leasqr.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,7 +1,7 @@
%% Copyright (C) 1992-1994 Richard Shrager
%% Copyright (C) 1992-1994 Arthur Jutan <ju...@ch...>
%% Copyright (C) 1992-1994 Ray Muzic <rf...@ds...>
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -90,7 +90,8 @@
## 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_f},
+## @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
@@ -176,7 +177,7 @@
##
## @code{plot_cmd}: Function enabling backend to plot results or
## intermediate results. Will be called with current computed
-## residualse. Default: plot nothing.
+## residuals. Default: plot nothing.
##
## @code{debug}: Logical scalar, default: @code{false}. Will be passed
## to the backend, which might print debugging information if true.
Modified: trunk/octave-forge/main/optim/inst/optim_problems.m
===================================================================
--- trunk/octave-forge/main/optim/inst/optim_problems.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/optim_problems.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,6 +1,6 @@
%% Copyright (C) 2007 Paul Kienzle (sort-based lookup in ODE solver)
-%% Copyright (C) 2009 Thomas Treichl <tho...@gm...>
-%% Copyright (C) 2010 Olaf Till <ola...@un...>
+%% Copyright (C) 2009 Thomas Treichl <tho...@gm...> (ode23 code)
+%% Copyright (C) 2010 Olaf Till <i7...@t-...>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
Modified: trunk/octave-forge/main/optim/inst/private/__dfdp__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -83,7 +83,7 @@
idxa = p == 0;
del(idxa) = diffp(idxa);
del(diff_onesided) = - del(diff_onesided); % keep course of
- % optimization of previous versions
+ % optimization of previous versions
absdel = abs (del);
idxd = ~(diff_onesided | fixed); % double sided interval
p1 = zeros (n, 1);
@@ -107,9 +107,9 @@
idxnvd = idxd & idxnv; % within bounds, double sided interval
%% remaining single sided intervals
p1(idxnvs) = p(idxnvs) + del(idxnvs); % don't take absdel, this could
- % change course of optimization without
- % bounds with respect to previous
- % versions
+ % change course of optimization without
+ % bounds with respect to previous
+ % versions
%% remaining single sided intervals, violating a bound -> take largest
%% possible direction of single sided interval
idxvs(idxnvs) = p1(idxnvs, 1) < lbound(idxnvs, 1) | ...
@@ -120,14 +120,14 @@
idx1g2w(idxvs) = idx1g2;
idx1le2w(idxvs) = ~idx1g2;
p1(idx1g2w) = max (p(idx1g2w, 1) - absdel(idx1g2w, 1), ...
- lbound(idx1g2w, 1));
+ lbound(idx1g2w, 1));
p1(idx1le2w) = min (p(idx1le2w, 1) + absdel(idx1le2w, 1), ...
- ubound(idx1le2w, 1));
+ ubound(idx1le2w, 1));
%% double sided interval
p1(idxnvd) = min (p(idxnvd, 1) + absdel(idxnvd, 1), ...
- ubound(idxnvd, 1));
+ ubound(idxnvd, 1));
p2(idxnvd) = max (p(idxnvd, 1) - absdel(idxnvd, 1), ...
- lbound(idxnvd, 1));
+ lbound(idxnvd, 1));
del(idxs) = p1(idxs) - p(idxs);
del(idxd) = p1(idxd) - p2(idxd);
@@ -141,16 +141,16 @@
ps = p;
ps(j) = p1(j);
if (idxs(j))
- info.side = 0; % onesided interval
- tp1 = func (ps, info);
- prt(:, j) = (tp1(:) - f) / del(j);
+ info.side = 0; % onesided interval
+ tp1 = func (ps, info);
+ prt(:, j) = (tp1(:) - f) / del(j);
else
- info.side = 1; % centered interval, side 1
- tp1 = func (ps, info);
- ps(j) = p2(j);
- info.side = 2; % centered interval, side 2
- tp2 = func (ps, info);
- prt(:, j) = (tp1(:) - tp2(:)) / del(j);
+ info.side = 1; % centered interval, side 1
+ tp1 = func (ps, info);
+ ps(j) = p2(j);
+ info.side = 2; % centered interval, side 2
+ tp2 = func (ps, info);
+ prt(:, j) = (tp1(:) - tp2(:)) / del(j);
end
end
end
Modified: trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -37,7 +37,7 @@
df_cstr = hook.df_cstr; % function of derivatives of all constraints
n_gencstr = hook.n_gencstr; % number of non-linear constraints
eq_idx = hook.eq_idx; % logical index of equality constraints in all
- % constraints
+ % constraints
lbound = hook.lbound; % bounds, subset of linear inequality
ubound = hook.ubound; % constraints in mc and vc
@@ -86,16 +86,16 @@
wtl = wt(:);
nz = 20 * eps; % This is arbitrary. Constraint function will be
- % regarded as <= zero if less than nz.
+ % regarded as <= zero if less than nz.
%% backend-specific checking of options and constraints
if (have_constraints_except_bounds)
if (any (pin_cstr.inequ.lin_except_bounds < 0) || ...
- (n_gencstr > 0 && any (pin_cstr.inequ.gen < 0)))
+ (n_gencstr > 0 && any (pin_cstr.inequ.gen < 0)))
warning ('initial parameters violate inequality constraints');
end
if (any (abs (pin_cstr.equ.lin) >= nz) || ...
- (n_gencstr > 0 && any (abs (pin_cstr.equ.gen) >= nz)))
+ (n_gencstr > 0 && any (abs (pin_cstr.equ.gen) >= nz)))
warning ('initial parameters violate equality constraints');
end
end
@@ -128,7 +128,7 @@
%% fill constant fields of hook for derivative-functions; some fields
%% may be backend-specific
dfdp_hook.fixed = fixed; % this may be handled by the frontend, but
- % the backend still may add to it
+ % the backend still may add to it
%% set up for iterations
%%
@@ -156,18 +156,18 @@
v_cstr = f_cstr (p, ac_idx);
%% index of active constraints
c_act = v_cstr < nz | eq_idx; # equality constraints might be
- # violated at start
+ # violated at start
if (any (c_act))
if (n_gencstr > 0)
- dct = df_cstr (p, ac_idx, ...
- setfield (dfdp_hook, 'f', v_cstr));
- dct(:, fixed) = 0; % for user supplied dfdp; necessary?
- dc = dct.';
- dcat = dct(c_act, :);
+ dct = df_cstr (p, ac_idx, ...
+ setfield (dfdp_hook, 'f', v_cstr));
+ dct(:, fixed) = 0; % for user supplied dfdp; necessary?
+ dc = dct.';
+ dcat = dct(c_act, :);
else
- dcat = df_cstr (p, c_act, ...
- setfield (dfdp_hook, 'f', v_cstr));
- dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
+ dcat = df_cstr (p, c_act, ...
+ setfield (dfdp_hook, 'f', v_cstr));
+ dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
end
dca = dcat.';
end
@@ -207,46 +207,46 @@
%%% ser = 1 ./ sqrt((s.*s)+epsL);
se = sqrt ((s.*s) + epsL);
if (new_s)
- %% for testing
- ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
+ %% for testing
+ ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
else
- ser = 1 ./ se;
+ ser = 1 ./ se;
end
tp1 = (v * (g .* ser)) .* nrm;
if (any (c_act))
- deb_printf (testing, 'constraints are active:\n');
- deb_printf (testing, '%i\n', c_act);
- %% calculate chg by 'quadratic programming'
- nrme= diag (nrm);
- ser2 = diag (ser .* ser);
- mfc1 = nrme * v * ser2 * v.' * nrme;
- tp2 = mfc1 * dca;
- a_eq_idx = eq_idx(c_act);
- [lb, bidx, ridx, tbl] = cpiv (dcat * tp1, dcat * tp2, a_eq_idx);
- chg = tp1 + tp2(:, bidx) * lb; % if a parameter is 'fixed',
- % the respective component of chg should
- % be zero too, even here (with active
- % constraints)
- deb_printf (testing, 'change:\n');
- deb_printf (testing, '%e\n', chg);
- deb_printf (testing, '\n');
- %% indices for different types of constraints
- c_inact = ~c_act; % inactive constraints
- c_binding = nc_idx;
- c_binding(c_act) = bidx; % constraints selected binding
- c_unbinding = nc_idx;
- c_unbinding(c_act) = ridx; % constraints unselected binding
- c_nonbinding = c_act & ~(c_binding | c_unbinding); % constraints
- % selected non-binding
+ deb_printf (testing, 'constraints are active:\n');
+ deb_printf (testing, '%i\n', c_act);
+ %% calculate chg by 'quadratic programming'
+ nrme= diag (nrm);
+ ser2 = diag (ser .* ser);
+ mfc1 = nrme * v * ser2 * v.' * nrme;
+ tp2 = mfc1 * dca;
+ a_eq_idx = eq_idx(c_act);
+ [lb, bidx, ridx, tbl] = cpiv (dcat * tp1, dcat * tp2, a_eq_idx);
+ chg = tp1 + tp2(:, bidx) * lb; % if a parameter is 'fixed',
+ % the respective component of chg should
+ % be zero too, even here (with active
+ % constraints)
+ deb_printf (testing, 'change:\n');
+ deb_printf (testing, '%e\n', chg);
+ deb_printf (testing, '\n');
+ %% indices for different types of constraints
+ c_inact = ~c_act; % inactive constraints
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx; % constraints selected binding
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx; % constraints unselected binding
+ c_nonbinding = c_act & ~(c_binding | c_unbinding); % constraints
+ % selected non-binding
else
- %% chg is the Levenberg/Marquardt step
- chg = tp1;
- %% indices for different types of constraints
- c_inact = ac_idx; % inactive constraints consist of all
- % constraints
- c_binding = nc_idx;
- c_unbinding = nc_idx;
- c_nonbinding = nc_idx;
+ %% chg is the Levenberg/Marquardt step
+ chg = tp1;
+ %% indices for different types of constraints
+ c_inact = ac_idx; % inactive constraints consist of all
+ % constraints
+ c_binding = nc_idx;
+ c_unbinding = nc_idx;
+ c_nonbinding = nc_idx;
end
%% apply constraints to step width (since this is a
%% Levenberg/Marquardt algorithm, no line-search is performed
@@ -258,200 +258,200 @@
hstep = mcit * chg;
idx = hstep < 0;
if (any (idx))
- k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
- hstep(idx)));
+ k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
+ hstep(idx)));
end
if (k < 1)
- deb_printf (testing, 'stepwidth: linear constraints\n');
+ deb_printf (testing, 'stepwidth: linear constraints\n');
end
if (n_gencstr > 0)
- c_tp = gc_idx & (c_nonbinding | c_inact);
- if (any (c_tp) && any (f_cstr (pprev + k * chg, c_tp) < 0))
- [k, fval, info] = ...
- fzero (@ (x) min (cat (1, ...
- f_cstr (pprev + x * chg, c_tp), ...
- k - x, ...
- ifelse (x < 0, -Inf, Inf))), ...
- 0);
- if (info ~= 1 || abs (fval) >= nz)
- error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
- end
- deb_printf (testing, 'general constraints limit stepwidth\n');
- end
+ c_tp = gc_idx & (c_nonbinding | c_inact);
+ if (any (c_tp) && any (f_cstr (pprev + k * chg, c_tp) < 0))
+ [k, fval, info] = ...
+ fzero (@ (x) min (cat (1, ...
+ f_cstr (pprev + x * chg, c_tp), ...
+ k - x, ...
+ ifelse (x < 0, -Inf, Inf))), ...
+ 0);
+ if (info ~= 1 || abs (fval) >= nz)
+ error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
+ end
+ deb_printf (testing, 'general constraints limit stepwidth\n');
+ end
end
chg = k * chg;
if (any (gc_idx & c_binding)) % none selected binding =>
- % none unselected binding
- deb_printf (testing, 'general binding constraints must be regained:\n');
- %% regain binding constraints and one of the possibly active
- %% previously inactive or non-binding constraints
- ptp1 = pprev + chg;
+ % none unselected binding
+ deb_printf (testing, 'general binding constraints must be regained:\n');
+ %% regain binding constraints and one of the possibly active
+ %% previously inactive or non-binding constraints
+ ptp1 = pprev + chg;
- tp = true;
- nt_nosuc = true;
- lim = 20;
- while (nt_nosuc && lim >= 0)
- deb_printf (testing, 'starting from new value of p in regaining:\n');
- deb_printf (testing, '%e\n', ptp1);
- %% we keep d_p.' * inv (mfc1) * d_p minimal in each step of
- %% the inner loop; this is both sensible (this metric
- %% considers a guess of curvature of sum of squared residuals)
- %% and convenient (we have useful matrices available for it)
- c_tp0 = c_inact | c_nonbinding;
- c_tp1 = c_inact | (gc_idx & c_nonbinding);
- btbl = tbl(bidx, bidx);
- c_tp2 = c_binding;
- if (any (tp)) % if none before, does not get true again
- tp = f_cstr (ptp1, c_tp1) < nz;
- if (any (tp)) % could be less clumsy, but ml-compatibility..
- %% keep only the first true entry in tp
- tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
- %% supplement binding index with one (the first) getting
- %% binding in c_tp1
- c_tp2(c_tp1) = tp;
- %% gradient of this added constraint
- caddt = dct(c_tp2 & ~c_binding, :);
- cadd = caddt.';
- C = dct(c_binding, :) * mfc1 * cadd;
- Ct = C.';
- G = [btbl, btbl * C; ...
- -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
- btbl = gjp (G, size (G, 1));
- end
- end
- dcbt = dct(c_tp2, :);
- mfc = - mfc1 * dcbt.' * btbl;
- deb_printf (testing, 'constraints to regain:\n');
- deb_printf (testing, '%i\n', c_tp2);
+ tp = true;
+ nt_nosuc = true;
+ lim = 20;
+ while (nt_nosuc && lim >= 0)
+ deb_printf (testing, 'starting from new value of p in regaining:\n');
+ deb_printf (testing, '%e\n', ptp1);
+ %% we keep d_p.' * inv (mfc1) * d_p minimal in each step of
+ %% the inner loop; this is both sensible (this metric
+ %% considers a guess of curvature of sum of squared residuals)
+ %% and convenient (we have useful matrices available for it)
+ c_tp0 = c_inact | c_nonbinding;
+ c_tp1 = c_inact | (gc_idx & c_nonbinding);
+ btbl = tbl(bidx, bidx);
+ c_tp2 = c_binding;
+ if (any (tp)) % if none before, does not get true again
+ tp = f_cstr (ptp1, c_tp1) < nz;
+ if (any (tp)) % could be less clumsy, but ml-compatibility..
+ %% keep only the first true entry in tp
+ tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
+ %% supplement binding index with one (the first) getting
+ %% binding in c_tp1
+ c_tp2(c_tp1) = tp;
+ %% gradient of this added constraint
+ caddt = dct(c_tp2 & ~c_binding, :);
+ cadd = caddt.';
+ C = dct(c_binding, :) * mfc1 * cadd;
+ Ct = C.';
+ G = [btbl, btbl * C; ...
+ -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
+ btbl = gjp (G, size (G, 1));
+ end
+ end
+ dcbt = dct(c_tp2, :);
+ mfc = - mfc1 * dcbt.' * btbl;
+ deb_printf (testing, 'constraints to regain:\n');
+ deb_printf (testing, '%i\n', c_tp2);
- ptp2 = ptp1;
- nt_niter_start = 100;
- nt_niter = nt_niter_start;
- while (nt_nosuc && nt_niter >= 0)
- hv = f_cstr (ptp2, c_tp2);
- if (all (abs (hv) < nz))
- nt_nosuc = false;
- chg = ptp2 - pprev;
- else
- ptp2 = ptp2 + mfc * hv; % step should be zero for each
- % component for which the parameter is
- % 'fixed'
- end
- nt_niter = nt_niter - 1;
- end
- deb_printf (testing, 'constraints after regaining:\n');
- deb_printf (testing, '%e\n', hv);
- if (nt_nosuc || ...
- any (abs (chg) > abs (pprev .* maxstep)) || ...
- any (f_cstr (ptp2, c_tp0) < -nz))
- if (nt_nosuc)
- deb_printf (testing, 'regaining did not converge\n');
- else
- deb_printf (testing, 'regaining violated type 3 and 4\n');
- end
- nt_nosuc = true;
- ptp1 = (pprev + ptp1) / 2;
- end
- if (~nt_nosuc)
- tp = f_cstr (ptp2, c_unbinding);
- if (any (tp) < 0) % again ml-compatibility clumsyness..
- [discarded, id] = min(tp);
- tid = find (ridx);
- id = tid(id); % index within active constraints
- unsuccessful_exchange = false;
- if (abs (tbl(id, id)) < nz) % Bard: not absolute value
- %% exchange this unselected binding constraint against a
- %% binding constraint, but not against an equality
- %% constraint
- tbidx = bidx & ~a_eq_idx;
- if (~any (tbidx))
- unsuccessful_exchange = true;
- else
- [discarded, idm] = max (abs (tbl(tbidx, id)));
- tid = find (tbidx);
- idm = tid(idm); % -> index within active constraints
- tbl = gjp (tbl, idm);
- bidx(idm) = false;
- ridx(idm) = true;
- end
- end
- if (unsuccessful_exchange)
- %% It probably doesn't look good now; this desperate
- %% last attempt is not in the original algortithm, since
- %% that didn't account for equality constraints.
- ptp1 = (pprev + ptp1) / 2;
- else
- tbl = gjp (tbl, id);
- bidx(id) = true;
- ridx(id) = false;
- c_binding = nc_idx;
- c_binding(c_act) = bidx;
- c_unbinding = nc_idx;
- c_unbinding(c_act) = ridx;
- end
- nt_nosuc = true;
- deb_printf (testing, 'regaining violated type 2\n');
- end
- end
- if (~nt_nosuc)
- deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
- nt_niter_start - nt_niter);
- deb_printf (testing, '%e\n', ptp2);
- end
- lim = lim - 1;
- end
- if (nt_nosuc)
- error ('could not regain binding constraints');
- end
+ ptp2 = ptp1;
+ nt_niter_start = 100;
+ nt_niter = nt_niter_start;
+ while (nt_nosuc && nt_niter >= 0)
+ hv = f_cstr (ptp2, c_tp2);
+ if (all (abs (hv) < nz))
+ nt_nosuc = false;
+ chg = ptp2 - pprev;
+ else
+ ptp2 = ptp2 + mfc * hv; % step should be zero for each
+ % component for which the parameter is
+ % 'fixed'
+ end
+ nt_niter = nt_niter - 1;
+ end
+ deb_printf (testing, 'constraints after regaining:\n');
+ deb_printf (testing, '%e\n', hv);
+ if (nt_nosuc || ...
+ any (abs (chg) > abs (pprev .* maxstep)) || ...
+ any (f_cstr (ptp2, c_tp0) < -nz))
+ if (nt_nosuc)
+ deb_printf (testing, 'regaining did not converge\n');
+ else
+ deb_printf (testing, 'regaining violated type 3 and 4\n');
+ end
+ nt_nosuc = true;
+ ptp1 = (pprev + ptp1) / 2;
+ end
+ if (~nt_nosuc)
+ tp = f_cstr (ptp2, c_unbinding);
+ if (any (tp) < 0) % again ml-compatibility clumsyness..
+ [discarded, id] = min(tp);
+ tid = find (ridx);
+ id = tid(id); % index within active constraints
+ unsuccessful_exchange = false;
+ if (abs (tbl(id, id)) < nz) % Bard: not absolute value
+ %% exchange this unselected binding constraint against a
+ %% binding constraint, but not against an equality
+ %% constraint
+ tbidx = bidx & ~a_eq_idx;
+ if (~any (tbidx))
+ unsuccessful_exchange = true;
+ else
+ [discarded, idm] = max (abs (tbl(tbidx, id)));
+ tid = find (tbidx);
+ idm = tid(idm); % -> index within active constraints
+ tbl = gjp (tbl, idm);
+ bidx(idm) = false;
+ ridx(idm) = true;
+ end
+ end
+ if (unsuccessful_exchange)
+ %% It probably doesn't look good now; this desperate
+ %% last attempt is not in the original algortithm, since
+ %% that didn't account for equality constraints.
+ ptp1 = (pprev + ptp1) / 2;
+ else
+ tbl = gjp (tbl, id);
+ bidx(id) = true;
+ ridx(id) = false;
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx;
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx;
+ end
+ nt_nosuc = true;
+ deb_printf (testing, 'regaining violated type 2\n');
+ end
+ end
+ if (~nt_nosuc)
+ deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
+ nt_niter_start - nt_niter);
+ deb_printf (testing, '%e\n', ptp2);
+ end
+ lim = lim - 1;
+ end
+ if (nt_nosuc)
+ error ('could not regain binding constraints');
+ end
else
- %% check the maximal stepwidth and apply as necessary
- ochg=chg;
- idx = ~isinf(maxstep);
- limit = abs(maxstep(idx).*pprev(idx));
- chg(idx) = min(max(chg(idx),-limit),limit);
- if (verbose && any(ochg ~= chg))
- disp(['Change in parameter(s): ', ...
- sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']);
- end
+ %% check the maximal stepwidth and apply as necessary
+ ochg=chg;
+ idx = ~isinf(maxstep);
+ limit = abs(maxstep(idx).*pprev(idx));
+ chg(idx) = min(max(chg(idx),-limit),limit);
+ if (verbose && any(ochg ~= chg))
+ disp(['Change in parameter(s): ', ...
+ sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']);
+ end
end
aprec=abs(pprec.*pbest); %---
%% ss=scalar sum of squares=sum((wt.*f)^2).
if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating
- % function if there is some non-miniscule
- % change
- p=chg+pprev;
- %% since the projection method may have slightly violated
- %% constraints due to inaccuracy, correct parameters to bounds
- %% --- but only if no further constraints are given, otherwise
- %% the inaccuracy in honoring them might increase by this
- if (~have_constraints_except_bounds)
- lidx = p < lbound;
- uidx = p > ubound;
- p(lidx, 1) = lbound(lidx, 1);
- p(uidx, 1) = ubound(uidx, 1);
- chg(lidx, 1) = p(lidx, 1) - pprev(lidx, 1);
- chg(uidx, 1) = p(uidx, 1) - pprev(uidx, 1);
- end
- %%
- f = F (p);
- r = wt .* f;
- r = r(:);
- if (~isreal (r))
- error ('weighted residuals are not real');
- end
- ss = r.' * r;
- deb_printf (testing, 'sbest: %.16e\n', sbest);
- deb_printf (testing, 'sgoal: %.16e\n', sgoal);
- deb_printf (testing, ' ss: %.16e\n', ss);
- if (ss<sbest)
+ % function if there is some non-miniscule
+ % change
+ p=chg+pprev;
+ %% since the projection method may have slightly violated
+ %% constraints due to inaccuracy, correct parameters to bounds
+ %% --- but only if no further constraints are given, otherwise
+ %% the inaccuracy in honoring them might increase by this
+ if (~have_constraints_except_bounds)
+ lidx = p < lbound;
+ uidx = p > ubound;
+ p(lidx, 1) = lbound(lidx, 1);
+ p(uidx, 1) = ubound(uidx, 1);
+ chg(lidx, 1) = p(lidx, 1) - pprev(lidx, 1);
+ chg(uidx, 1) = p(uidx, 1) - pprev(uidx, 1);
+ end
+ %%
+ f = F (p);
+ r = wt .* f;
+ r = r(:);
+ if (~isreal (r))
+ error ('weighted residuals are not real');
+ end
+ ss = r.' * r;
+ deb_printf (testing, 'sbest: %.16e\n', sbest);
+ deb_printf (testing, 'sgoal: %.16e\n', sgoal);
+ deb_printf (testing, ' ss: %.16e\n', ss);
+ if (ss<sbest)
pbest=p;
fbest=f;
sbest=ss;
- end
- if (ss<=sgoal)
+ end
+ if (ss<=sgoal)
break;
- end
+ end
end %---
end
%% printf ('epsL no.: %i\n', jjj); % for testing
@@ -472,7 +472,7 @@
if (all(abs(chg) <= aprec) && all(abs(chgprev) <= aprec))
cvg = 2;
if (verbose)
- fprintf('Parameter changes converged to specified precision\n');
+ fprintf('Parameter changes converged to specified precision\n');
end
break;
else
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -31,6 +31,12 @@
optimget = @ __optimget__;
endif
+ if (compare_versions (version (), "3.2.4", "<="))
+ ## For bug #31484; but Octave 3.6... shows bug #36288 due to this
+ ## workaround. Octave 3.7... seems to be all right.
+ __dfdp__ = @ __dfdp__;
+ endif
+
## some scalar defaults; some defaults are backend specific, so
## lacking elements in respective constructed vectors will be set to
## NA here in the frontend
@@ -56,7 +62,7 @@
"fract_prec", [], \
"diffp", [], \
"diff_onesided", [], \
- "complex_step_derivative", false, \
+ "complex_step_derivative_f", false, \
"complex_step_derivative_inequc", false, \
"complex_step_derivative_equc", false, \
"cstep", cstep_default, \
@@ -112,10 +118,10 @@
diffp = optimget (settings, "diffp");
diff_onesided = optimget (settings, "diff_onesided");
fixed = optimget (settings, "fixed");
- do_cstep = optimget (settings, "complex_step_derivative", false);
+ do_cstep = optimget (settings, "complex_step_derivative_f", false);
cstep = optimget (settings, "cstep", cstep_default);
if (do_cstep && ! isempty (dfdp))
- error ("both 'complex_step_derivative' and 'dfdp' are set");
+ error ("both 'complex_step_derivative_f' and 'dfdp' are set");
endif
do_cstep_inequc = \
optimget (settings, "complex_step_derivative_inequc", false);
@@ -458,7 +464,6 @@
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
@@ -814,8 +819,8 @@
tp = eye (sum (nonfixed));
lidx = lbound != - Inf;
uidx = ubound != Inf;
- mc = cat (2, mc, tp(:, lidx), - tp(:, uidx));
- vc = cat (1, vc, - lbound(lidx, 1), ubound(uidx, 1));
+ mc = cat (2, tp(:, lidx), - tp(:, uidx), mc);
+ vc = cat (1, - lbound(lidx, 1), ubound(uidx, 1), vc);
## concatenate linear inequality and equality constraints
mc = cat (2, mc, emc);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-05-24 06:03:58
|
Revision: 10512
http://octave.svn.sourceforge.net/octave/?rev=10512&view=rev
Author: i7tiol
Date: 2012-05-24 06:03:47 +0000 (Thu, 24 May 2012)
Log Message:
-----------
Add functions for constrained optimization of scalar-valued objective function.
Modified Paths:
--------------
trunk/octave-forge/main/optim/PKG_ADD
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/nonlin_min.m
trunk/octave-forge/main/optim/inst/private/__s2mat__.m
trunk/octave-forge/main/optim/inst/private/__siman__.m
trunk/octave-forge/main/optim/inst/private/__sqp__.m
Modified: trunk/octave-forge/main/optim/PKG_ADD
===================================================================
--- trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 05:54:36 UTC (rev 10511)
+++ trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 06:03:47 UTC (rev 10512)
@@ -3,5 +3,5 @@
## sometime in 3.3.54+, if I remember right
__all_opts__ ("nonlin_residmin");
__all_opts__ ("residmin_stat");
- ## __all_opts__ ("nonlin_min");
+ __all_opts__ ("nonlin_min");
endif
Added: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-05-24 06:03:47 UTC (rev 10512)
@@ -0,0 +1,1354 @@
+## Copyright (C) 2012 Olaf Till <i7...@t-...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin})
+## @deftypefnx {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin}, @var{settings})
+##
+## Frontend for nonlinear minimization of a scalar objective function.
+## The functions supplied by the user have a minimal interface; any
+## additionally needed constants can be supplied by wrapping the user
+## functions into anonymous functions.
+##
+## At the moment there is only a stochastic backend implemented, so that
+## some settings in the frontend have no use.
+##
+## The following description applies to usage with vector-based
+## parameter handling. Differences in usage for structure-based
+## parameter handling will be explained in a separate section below.
+##
+## @var{f}: objective function. It gets a column vector of real
+## parameters as argument. In gradient determination, this function may
+## be called with an informational second argument, whose content
+## depends on the function for gradient determination.
+##
+## @var{pin}: real column vector of initial parameters.
+##
+## @var{settings}: structure whose fields stand for optional settings
+## referred to below. The fields can be set by @code{optimset()} with
+## Octave versions 3.3.55 or greater; with older Octave versions, the
+## fields must be set directly as structure-fields in the correct case.
+##
+## The returned values are the column vector of final parameters
+## @var{p}, the final value of the objective function @var{objf}, an
+## integer @var{cvg} indicating if and how optimization succeeded or
+## failed, and a structure @var{outp} with additional information,
+## curently with only one field: @var{niter}, the number of iterations.
+## @var{cvg} is greater than zero for success and less than or equal to
+## zero for failure; its possible values depend on the used backend and
+## currently can be @code{0} (maximum number of iterations exceeded),
+## @code{1} (fixed number of iterations completed, e.g. in stochastic
+## optimizers), @code{2} (parameter change less than specified precision
+## in two consecutive iterations), or @code{3} (improvement in objective
+## function less than specified).
+##
+## @var{settings}:
+##
+## @code{Algorithm}: String specifying the backend. Currently there is
+## no default, since only a stochastic backend is implemented
+## (@code{"siman"}). It is described in a separate section below.
+##
+## @code{objf_grad}: Function computing the gradient of the objective
+## function with respect to the parameters, assuming residuals are
+## reshaped to a vector. Default: finite differences. Will be called
+## with the column vector of parameters and an informational structure
+## as arguments. The structure has the fields @code{f}: value of
+## objective function for current parameters, @code{fixed}: logical
+## vector indicating which parameters are not optimized, so these
+## partial derivatives need not be computed and can be set to zero,
+## @code{diffp}, @code{diff_onesided}, @code{lbound}, @code{ubound}:
+## identical to the user settings of this name, @code{plabels}:
+## 1-dimensional cell-array of column-cell-arrays, each column with
+## labels for all parameters, the first column contains the numerical
+## indices of the parameters. The default gradient function will call
+## the objective function with the second argument set with fields
+## @code{f}: as the @code{f} passed to the gradient function,
+## @code{plabels}: cell-array of 1x1 cell-arrays with the entries of the
+## column-cell-arrays of @code{plabels} as passed to the jacobian
+## function corresponding to current parameter, @code{side}: @code{0}
+## for one-sided interval, @code{1} or @code{2}, respectively, for the
+## sides of a two-sided interval, and @code{parallel}: logical scalar
+## indicating parallel computation of partial derivatives.
+##
+## @code{objf_hessian}: Function computing the Hessian of the objective
+## function with respect to the parameters. The default is backend
+## specific. Will be called with the column vector of parameters as
+## argument.
+##
+## @code{diffp}: column vector of fractional intervals (doubled for
+## central intervals) supposed to be used by gradient functions
+## performing finite differencing. Default: @code{.001 * ones (size
+## (parameters))}. The default gradient function will use these as
+## absolute intervals for parameters with value zero.
+##
+## @code{diff_onesided}: logical column vector indicating that one-sided
+## intervals should be used by gradient functions performing finite
+## differencing. Default: @code{false (size (parameters))}.
+##
+## @code{complex_step_derivative_objf},
+## @code{complex_step_derivative_inequc},
+## @code{complex_step_derivative_equc}: logical scalars, default: false.
+## Estimate gradient of objective function, general inequality
+## constraints, and general equality constraints, respectively, with
+## complex step derivative approximation. Use only if you know that your
+## objective function, function of general inequality constraints, or
+## function of general equality constraints, respectively, is suitable
+## for this. No user function for the respective gradient 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
+## additional parameters under certain conditions.
+##
+## @code{lbound}, @code{ubound}: column vectors of lower and upper
+## bounds for parameters. Default: @code{-Inf} and @code{+Inf},
+## respectively. The bounds are non-strict, i.e. parameters are allowed
+## to be exactly equal to a bound. The default gradient function will
+## respect bounds (but no further inequality constraints) in finite
+## differencing.
+##
+## @code{inequc}: Further inequality constraints. Cell-array containing
+## up to four entries, two entries for linear inequality constraints
+## and/or one or two entries for general inequality constraints. Either
+## linear or general constraints may be the first entries, but the two
+## entries for linear constraints must be adjacent and, if two entries
+## are given for general constraints, they also must be adjacent. The
+## two entries for linear constraints are a matrix (say @code{m}) and a
+## vector (say @code{v}), specifying linear inequality constraints of
+## the form @code{m.' * parameters + v >= 0}. The first entry for
+## general constraints must be a differentiable vector valued function
+## (say @code{h}), specifying general inequality constraints of the form
+## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
+## paraters and the optional argument @code{idx} is a logical index.
+## @code{h} has to return the values of all constraints if @code{idx} is
+## not given, and has to return only the indexed constraints if
+## @code{idx} is given (so computation of the other constraints can be
+## spared). In gradient determination, this function may be called with
+## an informational third argument, whose content depends on the
+## function for gradient determination. If a second entry for general
+## inequality constraints is given, it must be a function computing the
+## jacobian of the constraints with respect to the parameters. For this
+## function, the description of @code{dfdp} above applies, except that
+## it is called with 3 arguments since it has an additional argument
+## @code{idx} --- a logical index --- at second position, indicating
+## which rows of the jacobian must be returned, and except that the
+## default function calls @code{h} with 3 arguments, since the argument
+## @code{idx} is also supplied. Note that specifying linear constraints
+## as general constraints will generally waste performance, even if
+## further, non-linear, general constraints are also specified.
+##
+## @code{equc}: Equality constraints. Specified the same way as
+## inequality constraints (see @code{inequc}).
+##
+## @code{cpiv}: Function for complementary pivoting, usable in
+## algorithms for constraints. Default: @ cpiv_bard. Only the default
+## function is supplied with the package.
+##
+## @code{TolFun}: Minimum fractional improvement in objective function
+## in an iteration (abortion criterium). Default: .0001.
+##
+## @code{MaxIter}: Maximum number of iterations (abortion criterium).
+## Default: backend-specific.
+##
+## @code{fract_prec}: Column Vector, minimum fractional change of
+## parameters in an iteration (abortion criterium if violated in two
+## consecutive iterations). Default: backend-specific.
+##
+## @code{max_fract_change}: Column Vector, enforced maximum fractional
+## change in parameters in an iteration. Default: backend-specific.
+##
+## @code{Display}: String indicating the degree of verbosity. Default:
+## @code{"off"}. Possible values are currently @code{"off"} (no
+## messages) and @code{"iter"} (some messages after each iteration).
+## Support of this setting and its exact interpretation are
+## backend-specific.
+##
+## @code{debug}: Logical scalar, default: @code{false}. Will be passed
+## to the backend, which might print debugging information if true.
+##
+## Structure-based parameter handling
+##
+## The setting @code{param_order} is a cell-array with names of the
+## optimized parameters. If not given, and initial parameters are a
+## structure, all parameters in the structure are optimized. If initial
+## parameters are a structure, it is an error if @code{param_order} is
+## not given and there are any non-structure-based configuration items
+## or functions.
+##
+## The initial parameters @var{pin} can be given as a structure
+## containing at least all fields named in @code{param_order}. In this
+## case the returned parameters @var{p} will also be a structure.
+##
+## Each user-supplied function can be called with the argument
+## containing the current parameters being a structure instead of a
+## column vector. For this, a corresponding setting must be set to
+## @code{true}: @code{objf_pstruct} (objective function),
+## @code{objf_grad_pstruct} (gradient of objective function),
+## @code{objf_hessian_pstruct} (hessian of objective function),
+## @code{f_inequc_pstruct} (general inequality constraints),
+## @code{df_inequc_pstruct} (jacobian of general inequality
+## constraints), @code{f_equc_pstruct} (general equality constraints),
+## and @code{df_equc_pstruct} (jacobian of general equality
+## constraints). If a gradient (jacobian) function is configured in such
+## a way, it must return the entries (columns) of the gradient
+## (jacobian) as fields of a structure under the respective parameter
+## names. If the hessian function is configured in such a way, it must
+## return a structure (say @code{h}) with fields e.g. as @code{h.a.b =
+## value} for @{value} being the 2nd partial derivative with respect to
+## @code{a} and @code{b}. There is no need to also specify the field
+## @code{h.b.a} in this example.
+##
+## Similarly, for specifying linear constraints, instead of the matrix
+## (called @code{m} above), a structure containing the rows of the
+## matrix in fields under the respective parameter names can be given.
+## In this case, rows containing only zeros need not be given.
+##
+## The vector-based settings @code{lbound}, @code{ubound},
+## @code{fixed}, @code{diffp}, @code{diff_onesided}, @code{fract_prec},
+## and @code{max_fract_change} can be replaced by the setting
+## @code{param_config}. It is a structure that can contain fields named
+## in @code{param_order}. For each such field, there may be subfields
+## with the same names as the above vector-based settings, but
+## containing a scalar value for the respective parameter. If
+## @code{param_config} is specified, none of the above
+## vector/matrix-based settings may be used.
+##
+## Additionally, named parameters are allowed to be non-scalar real
+## arrays. In this case, their dimensions are given by the setting
+## @code{param_dims}, a cell-array of dimension vectors, each containing
+## at least two dimensions; if not given, dimensions are taken from the
+## initial parameters, if these are given in a structure. Any
+## vector-based settings or not structure-based linear constraints then
+## must correspond to an order of parameters with all parameters
+## reshaped to vectors and concatenated in the user-given order of
+## parameter names. Structure-based settings or structure-based initial
+## parameters must contain arrays with dimensions reshapable to those of
+## the respective parameters.
+##
+## Description of backends (currently only one)
+##
+## "siman"
+##
+## A simulated annealing (stochastic) optimizer, changing all parameters
+## at once in a single step, so being suitable for non-bound
+## constraints.
+##
+## No gradient or hessian of the objective function is used. The
+## settings @code{MaxIter}, @code{fract_prec}, @code{TolFun}, and
+## @code{max_fract_change} are not honoured.
+##
+## Accepts the additional settings @code{T_init} (initial temperature,
+## default 0.01), @code{T_min} (final temperature, default 1.0e-5),
+## @code{mu_T} (factor of temperature decrease, default 1.005),
+## @code{iters_fixed_T} (iterations within one temperature step, default
+## 10), @code{max_rand_step} (column vector or structure-based
+## configuration of maximum random steps for each parameter, default
+## 0.005 * @var{pin}), @code{stoch_regain_constr} (if @code{true},
+## regain constraints after a random step, otherwise take new random
+## value until constraints are met, default false), @code{trace_steps}
+## (set field @code{trace} of @var{outp} with a matrix with a row for
+## each step, first column iteration number, second column repeat number
+## within iteration, third column value of objective function, rest
+## columns parameter values, default false), and @code{siman_log} (set
+## field @code{log} of @var{outp} with a matrix with a row for each
+## iteration, first column temperature, second column value of objective
+## function, rest columns numbers of tries with decrease, no decrease
+## but accepted, and no decrease and rejected.
+##
+## Steps with increase @code{diff} of objective function are accepted if
+## @code{rand (1) < exp (- diff / T)}, where @code{T} is the temperature
+## of the current iteration.
+##
+## If regaining of constraints failed, optimization will be aborted and
+## returned value of @var{cvg} will be @code{0}. Otherwise, @var{cvg}
+## will be @code{1}.
+##
+## Interpretation of @code{Display}: if set to @code{"iter"}, an
+## informational line is printed after each iteration.
+##
+## @end deftypefn
+
+## disabled PKG_ADD: __all_opts__ ("nonlin_min");
+
+function [p, objf, cvg, outp] = nonlin_min (f, pin, settings)
+
+ 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
+ optimget = @ __optimget__;
+ endif
+
+ if (compare_versions (version (), "3.2.4", "<="))
+ ## For bug #31484; but Octave 3.6... shows bug #36288 due to this
+ ## workaround. Octave 3.7... seems to be all right.
+ __dfdp__ = @ __dfdp__;
+ endif
+
+ ## some scalar defaults; some defaults are backend specific, so
+ ## lacking elements in respective constructed vectors will be set to
+ ## 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", [], \
+ "param_order", [], \
+ "param_dims", [], \
+ "f_inequc_pstruct", false, \
+ "f_equc_pstruct", false, \
+ "objf_pstruct", false, \
+ "df_inequc_pstruct", false, \
+ "df_equc_pstruct", false, \
+ "objf_grad_pstruct", false, \
+ "objf_hessian_pstruct", false, \
+ "lbound", [], \
+ "ubound", [], \
+ "objf_grad", [], \
+ "objf_hessian", [], \
+ "cpiv", @ cpiv_bard, \
+ "max_fract_change", [], \
+ "fract_prec", [], \
+ "diffp", [], \
+ "diff_onesided", [], \
+ "complex_step_derivative_objf", false, \
+ "complex_step_derivative_inequc", false, \
+ "complex_step_derivative_equc", false, \
+ "cstep", cstep_default, \
+ "fixed", [], \
+ "inequc", [], \
+ "equc", [], \
+ "TolFun", stol_default, \
+ "MaxIter", [], \
+ "Display", "off", \
+ "Algorithm", [], \ # set reasonable default, once available
+ "T_init", .01, \
+ "T_min", 1.0e-5, \
+ "mu_T", 1.005, \
+ "iters_fixed_T", 10, \
+ "max_rand_step", [], \
+ "stoch_regain_constr", false, \
+ "trace_steps", false, \
+ "siman_log", false, \
+ "debug", false);
+ return;
+ endif
+
+ assign = @ assign; # Is this faster in repeated calls?
+
+ if (nargin != 3)
+ error ("incorrect number of arguments");
+ endif
+
+ if (ischar (f))
+ f = str2func (f);
+ endif
+
+ if (! (pin_struct = isstruct (pin)))
+ if (! isvector (pin) || columns (pin) > 1)
+ error ("initial parameters must be either a structure or a column vector");
+ endif
+ endif
+
+ #### processing of settings and consistency checks
+
+ pconf = optimget (settings, "param_config");
+ pord = optimget (settings, "param_order");
+ pdims = optimget (settings, "param_dims");
+ f_inequc_pstruct = optimget (settings, "f_inequc_pstruct", false);
+ f_equc_pstruct = optimget (settings, "f_equc_pstruct", false);
+ f_pstruct = optimget (settings, "objf_pstruct", false);
+ dfdp_pstruct = optimget (settings, "objf_grad_pstruct", f_pstruct);
+ hessian_pstruct = optimget (settings, "objf_hessian_pstruct", f_pstruct);
+ df_inequc_pstruct = optimget (settings, "df_inequc_pstruct", \
+ f_inequc_pstruct);
+ df_equc_pstruct = optimget (settings, "df_equc_pstruct", \
+ f_equc_pstruct);
+ lbound = optimget (settings, "lbound");
+ ubound = optimget (settings, "ubound");
+ dfdp = optimget (settings, "objf_grad");
+ if (ischar (dfdp)) dfdp = str2func (dfdp); endif
+ hessian = optimget (settings, "objf_hessian");
+ max_fract_change = optimget (settings, "max_fract_change");
+ fract_prec = optimget (settings, "fract_prec");
+ diffp = optimget (settings, "diffp");
+ diff_onesided = optimget (settings, "diff_onesided");
+ fixed = optimget (settings, "fixed");
+ do_cstep = optimget (settings, "complex_step_derivative_objf", false);
+ cstep = optimget (settings, "cstep", cstep_default);
+ if (do_cstep && ! isempty (dfdp))
+ error ("both 'complex_step_derivative_objf' and 'objf_grad' are set");
+ endif
+ do_cstep_inequc = \
+ optimget (settings, "complex_step_derivative_inequc", false);
+ do_cstep_equc = optimget (settings, "complex_step_derivative_equc", \
+ false);
+ max_rand_step = optimget (settings, "max_rand_step");
+
+ any_vector_conf = ! (isempty (lbound) && isempty (ubound) && \
+ isempty (max_fract_change) && \
+ isempty (fract_prec) && isempty (diffp) && \
+ isempty (diff_onesided) && isempty (fixed) && \
+ isempty (max_rand_step));
+
+ ## collect constraints
+ [mc, vc, f_genicstr, df_gencstr, user_df_gencstr] = \
+ __collect_constraints__ (optimget (settings, "inequc"), \
+ do_cstep_inequc, "inequality constraints");
+ [emc, evc, f_genecstr, df_genecstr, user_df_genecstr] = \
+ __collect_constraints__ (optimget (settings, "equc"), \
+ do_cstep_equc, "equality constraints");
+ mc_struct = isstruct (mc);
+ emc_struct = isstruct (emc);
+
+ ## correct "_pstruct" settings if functions are not supplied
+ if (isempty (dfdp)) dfdp_pstruct = false; endif
+ if (isempty (hessian)) hessian_pstruct = false; endif
+ if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
+ if (isempty (f_genecstr)) f_equc_pstruct = false; endif
+ if (! user_df_gencstr) df_inequc_pstruct = false; endif
+ if (! user_df_genecstr) df_equc_pstruct = false; endif
+
+ ## some settings require a parameter order
+ if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
+ f_equc_pstruct || f_pstruct || dfdp_pstruct || \
+ hessian_pstruct || df_inequc_pstruct || df_equc_pstruct || \
+ mc_struct || emc_struct)
+ if (isempty (pord))
+ if (pin_struct)
+ if (any_vector_conf || \
+ ! (f_pstruct && \
+ (f_inequc_pstruct || isempty (f_genicstr)) && \
+ (f_equc_pstruct || isempty (f_genecstr)) && \
+ (dfdp_pstruct || isempty (dfdp)) && \
+ (hessian_pstruct || isempty (hessian)) && \
+ (df_inequc_pstruct || ! user_df_gencstr) && \
+ (df_equc_pstruct || ! user_df_genecstr) && \
+ (mc_struct || isempty (mc)) && \
+ (emc_struct || isempty (emc))))
+ error ("no parameter order specified and constructing a parameter order from the structure of initial parameters can not be done since not all configuration or given functions are structure based");
+ else
+ pord = fieldnames (pin);
+ endif
+ else
+ error ("given settings require specification of parameter order or initial parameters in the form of a structure");
+ endif
+ endif
+ pord = pord(:);
+ if (pin_struct && ! all (isfield (pin, pord)))
+ error ("some initial parameters lacking");
+ endif
+ if ((nnames = rows (unique (pord))) < rows (pord))
+ error ("duplicate parameter names in 'param_order'");
+ endif
+ if (isempty (pdims))
+ if (pin_struct)
+ pdims = cellfun \
+ (@ size, fields2cell (pin, pord), "UniformOutput", false);
+ else
+ pdims = num2cell (ones (nnames, 2), 2);
+ endif
+ else
+ pdims = pdims(:);
+ if (pin_struct && \
+ ! all (cellfun (@ (x, y) prod (size (x)) == prod (y), \
+ struct2cell (pin), pdims)))
+ error ("given param_dims and dimensions of initial parameters do not match");
+ endif
+ endif
+ if (nnames != rows (pdims))
+ error ("lengths of 'param_order' and 'param_dims' not equal");
+ endif
+ pnel = cellfun (@ prod, pdims);
+ ppartidx = pnel;
+ if (any (pnel > 1))
+ pnonscalar = true;
+ cpnel = num2cell (pnel);
+ prepidx = cat (1, cellfun \
+ (@ (x, n) x(ones (1, n), 1), \
+ num2cell ((1:nnames).'), cpnel, \
+ "UniformOutput", false){:});
+ epord = pord(prepidx, 1);
+ psubidx = cat (1, cellfun \
+ (@ (n) (1:n).', cpnel, \
+ "UniformOutput", false){:});
+ else
+ pnonscalar = false; # some less expensive interfaces later
+ prepidx = (1:nnames).';
+ epord = pord;
+ psubidx = ones (nnames, 1);
+ endif
+ else
+ pord = []; # spares checks for given but not needed
+ endif
+
+ if (pin_struct)
+ np = sum (pnel);
+ else
+ np = length (pin);
+ if (! isempty (pord) && np != sum (pnel))
+ error ("number of initial parameters not correct");
+ endif
+ endif
+
+ plabels = num2cell (num2cell ((1:np).'));
+ if (! isempty (pord))
+ plabels = cat (2, plabels, num2cell (epord), \
+ num2cell (num2cell (psubidx)));
+ endif
+
+ ## some useful vectors
+ zerosvec = zeros (np, 1);
+ NAvec = NA (np, 1);
+ Infvec = Inf (np, 1);
+ falsevec = false (np, 1);
+ sizevec = [np, 1];
+
+ ## collect parameter-related configuration
+ if (! isempty (pconf))
+ ## use supplied configuration structure
+
+ ## parameter-related configuration is either allowed by a structure
+ ## or by vectors
+ if (any_vector_conf)
+ error ("if param_config is given, its potential items must not \
+ be configured in another way");
+ endif
+
+ ## supplement parameter names lacking in param_config
+ nidx = ! isfield (pconf, pord);
+ pconf = cell2fields ({struct()}(ones (1, sum (nidx))), \
+ pord(nidx), 2, pconf);
+
+ pconf = structcat (1, fields2cell (pconf, pord){:});
+
+ ## in the following, use reshape with explicit dimensions (instead
+ ## of x(:)) so that errors are thrown if a configuration item has
+ ## incorrect number of elements
+
+ lbound = - Infvec;
+ if (isfield (pconf, "lbound"))
+ idx = ! fieldempty (pconf, "lbound");
+ if (pnonscalar)
+ lbound (idx(prepidx), 1) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).lbound}.', \
+ cpnel(idx), "UniformOutput", false){:});
+ else
+ lbound(idx, 1) = cat (1, pconf.lbound);
+ endif
+ endif
+
+ ubound = Infvec;
+ if (isfield (pconf, "ubound"))
+ idx = ! fieldempty (pconf, "ubound");
+ if (pnonscalar)
+ ubound (idx(prepidx), 1) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).ubound}.', \
+ cpnel(idx), "UniformOutput", false){:});
+ else
+ ubound(idx, 1) = cat (1, pconf.ubound);
+ endif
+ endif
+
+ max_fract_change = fract_prec = NAvec;
+
+ if (isfield (pconf, "max_fract_change"))
+ idx = ! fieldempty (pconf, "max_fract_change");
+ if (pnonscalar)
+ max_fract_change(idx(prepidx)) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).max_fract_change}.', \
+ cpnel(idx), \
+ "UniformOutput", false){:});
+ else
+ max_fract_change(idx) = [pconf.max_fract_change];
+ endif
+ endif
+
+ if (isfield (pconf, "fract_prec"))
+ idx = ! fieldempty (pconf, "fract_prec");
+ if (pnonscalar)
+ fract_prec(idx(prepidx)) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).fract_prec}.', cpnel(idx), \
+ "UniformOutput", false){:});
+ else
+ fract_prec(idx) = [pconf.fract_prec];
+ endif
+ endif
+
+ diffp = zerosvec;
+ diffp(:) = diffp_default;
+ if (isfield (pconf, "diffp"))
+ idx = ! fieldempty (pconf, "diffp");
+ if (pnonscalar)
+ diffp(idx(prepidx)) = \
+ cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).diffp}.', cpnel(idx), \
+ "UniformOutput", false){:});
+ else
+ diffp(idx) = [pconf.diffp];
+ endif
+ endif
+
+ diff_onesided = fixed = falsevec;
+
+ if (isfield (pconf, "diff_onesided"))
+ idx = ! fieldempty (pconf, "diff_onesided");
+ if (pnonscalar)
+ diff_onesided(idx(prepidx)) = \
+ logical \
+ (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).diff_onesided}.', cpnel(idx), \
+ "UniformOutput", false){:}));
+ else
+ diff_onesided(idx) = logical ([pconf.diff_onesided]);
+ endif
+ endif
+
+ if (isfield (pconf, "fixed"))
+ idx = ! fieldempty (pconf, "fixed");
+ if (pnonscalar)
+ fixed(idx(prepidx)) = \
+ logical \
+ (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).fixed}.', cpnel(idx), \
+ "UniformOutput", false){:}));
+ else
+ fixed(idx) = logical ([pconf.fixed]);
+ endif
+ endif
+
+ max_rand_step = NAvec;
+
+ if (isfield (pconf, "max_rand_step"))
+ idx = ! fieldempty (pconf, "max_rand_step");
+ if (pnonscalar)
+ max_rand_step(idx(prepidx)) = \
+ logical \
+ (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ {pconf(idx).max_rand_step}.',
+ cpnel(idx), \
+ "UniformOutput", false){:}));
+ else
+ max_rand_step(idx) = logical ([pconf.max_rand_step]);
+ endif
+ endif
+
+ else
+ ## use supplied configuration vectors
+
+ if (isempty (lbound))
+ lbound = - Infvec;
+ elseif (any (size (lbound) != sizevec))
+ error ("bounds: wrong dimensions");
+ endif
+
+ if (isempty (ubound))
+ ubound = Infvec;
+ elseif (any (size (ubound) != sizevec))
+ error ("bounds: wrong dimensions");
+ endif
+
+ if (isempty (max_fract_change))
+ max_fract_change = NAvec;
+ elseif (any (size (max_fract_change) != sizevec))
+ error ("max_fract_change: wrong dimensions");
+ endif
+
+ if (isempty (fract_prec))
+ fract_prec = NAvec;
+ elseif (any (size (fract_prec) != sizevec))
+ error ("fract_prec: wrong dimensions");
+ endif
+
+ if (isempty (diffp))
+ diffp = zerosvec;
+ diffp(:) = diffp_default;
+ else
+ if (any (size (diffp) != sizevec))
+ error ("diffp: wrong dimensions");
+ endif
+ diffp(isna (diffp)) = diffp_default;
+ endif
+
+ if (isempty (diff_onesided))
+ diff_onesided = falsevec;
+ else
+ if (any (size (diff_onesided) != sizevec))
+ error ("diff_onesided: wrong dimensions")
+ endif
+ diff_onesided(isna (diff_onesided)) = false;
+ diff_onesided = logical (diff_onesided);
+ endif
+
+ if (isempty (fixed))
+ fixed = falsevec;
+ else
+ if (any (size (fixed) != sizevec))
+ error ("fixed: wrong dimensions");
+ endif
+ fixed(isna (fixed)) = false;
+ fixed = logical (fixed);
+ endif
+
+ if (isempty (max_rand_step))
+ max_rand_step = NAvec;
+ elseif (any (size (max_rand_step) != sizevec))
+ error ("max_rand_step: wrong dimensions");
+ endif
+
+ endif
+
+ ## guaranty all (lbound <= ubound)
+ if (any (lbound > ubound))
+ error ("some lower bounds larger than upper bounds");
+ endif
+
+ #### consider whether initial parameters and functions are based on
+ #### parameter structures or parameter vectors; wrappers for call to
+ #### default function for jacobians
+
+ ## initial parameters
+ if (pin_struct)
+ if (pnonscalar)
+ pin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
+ fields2cell (pin, pord), cpnel, \
+ "UniformOutput", false){:});
+ else
+ pin = cat (1, fields2cell (pin, pord){:});
+ endif
+ endif
+
+ ## objective function
+ if (f_pstruct)
+ if (pnonscalar)
+ f = @ (p, varargin) \
+ f (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), varargin{:});
+ else
+ f = @ (p, varargin) \
+ f (cell2struct (num2cell (p), pord, 1), varargin{:});
+ endif
+ endif
+ f_pin = f (pin);
+
+ ## gradient of objective function
+ if (isempty (dfdp))
+ if (do_cstep)
+ dfdp = @ (p, hook) jacobs (p, f, hook);
+ else
+ dfdp = @ (p, hook) __dfdp__ (p, f, hook);
+ endif
+ endif
+ if (dfdp_pstruct)
+ if (pnonscalar)
+ dfdp = @ (p, hook) \
+ cat (2, \
+ fields2cell \
+ (dfdp (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), hook), \
+ pord){:});
+ else
+ dfdp = @ (p, hook) \
+ cat (2, \
+ fields2cell \
+ (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
+ pord){:});
+ endif
+ endif
+
+ ## hessian of objective function
+ if (hessian_pstruct)
+ if (pnonscalar)
+ hessian = @ (p) \
+ hessian_struct2mat \
+ (hessian (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1)), pord);
+ else
+ hessian = @ (p) \
+ hessian_struct2mat \
+ (hessian (cell2struct (num2cell (p), pord, 1)), pord);
+ endif
+ endif
+
+ ## function for general inequality constraints
+ if (f_inequc_pstruct)
+ if (pnonscalar)
+ f_genicstr = @ (p, varargin) \
+ f_genicstr (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), varargin{:});
+ else
+ f_genicstr = @ (p, varargin) \
+ f_genicstr \
+ (cell2struct (num2cell (p), pord, 1), varargin{:});
+ endif
+ endif
+
+ ## note this stage
+ possibly_pstruct_f_genicstr = f_genicstr;
+
+ ## jacobian of general inequality constraints
+ if (df_inequc_pstruct)
+ if (pnonscalar)
+ df_gencstr = @ (p, func, idx, hook) \
+ cat (2, \
+ fields2cell \
+ (df_gencstr \
+ (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1), \
+ func, idx, hook), \
+ pord){:});
+ else
+ df_gencstr = @ (p, func, idx, hook) \
+ cat (2, \
+ fields2cell \
+ (df_gencstr (cell2struct (num2cell (p), pord, 1), \
+ func, idx, hook), \
+ pord){:});
+ endif
+ endif
+
+ ## function for general equality constraints
+ if (f_equc_pstruct)
+ if (pnonscalar)
+ f_genecstr = @ (p, varargin) \
+ f_genecstr (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1), varargin{:});
+ else
+ f_genecstr = @ (p, varargin) \
+ f_genecstr \
+ (cell2struct (num2cell (p), pord, 1), varargin{:});
+ endif
+ endif
+
+ ## note this stage
+ possibly_pstruct_f_genecstr = f_genecstr;
+
+ ## jacobian of general equality constraints
+ if (df_equc_pstruct)
+ if (pnonscalar)
+ df_genecstr = @ (p, func, idx, hook) \
+ cat (2, \
+ fields2cell \
+ (df_genecstr \
+ (cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1), \
+ func, idx, hook), \
+ pord){:});
+ else
+ df_genecstr = @ (p, func, idx, hook) \
+ cat (2, \
+ fields2cell \
+ (df_genecstr (cell2struct (num2cell (p), pord, 1), \
+ func, idx, hook), \
+ pord){:});
+ endif
+ endif
+
+ ## linear inequality constraints
+ if (mc_struct)
+ idx = isfield (mc, pord);
+ if (rows (fieldnames (mc)) > sum (idx))
+ error ("unknown fields in structure of linear inequality constraints");
+ endif
+ smc = mc;
+ mc = zeros (np, rows (vc));
+ mc(idx(prepidx), :) = cat (1, fields2cell (smc, pord(idx)){:});
+ endif
+
+ ## linear equality constraints
+ if (emc_struct)
+ idx = isfield (emc, pord);
+ if (rows (fieldnames (emc)) > sum (idx))
+ error ("unknown fields in structure of linear equality constraints");
+ endif
+ semc = emc;
+ emc = zeros (np, rows (evc));
+ emc(idx(prepidx), :) = cat (1, fields2cell (semc, pord(idx)){:});
+ endif
+
+ ## parameter-related configuration for jacobian functions
+ if (dfdp_pstruct || df_inequc_pstruct || df_equc_pstruct)
+ if(pnonscalar)
+ s_diffp = cell2struct \
+ (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_diff_onesided = cell2struct \
+ (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_orig_lbound = cell2struct \
+ (cellfun (@ reshape, mat2cell (lbound, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_orig_ubound = cell2struct \
+ (cellfun (@ reshape, mat2cell (ubound, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ s_plabels = cell2struct \
+ (num2cell \
+ (cat (2, cellfun \
+ (@ (x) cellfun \
+ (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
+ pdims, "UniformOutput", false), \
+ num2cell (plabels, 1), "UniformOutput", false){:}), \
+ 2), \
+ pord, 1);
+ s_orig_fixed = cell2struct \
+ (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
+ pdims, "UniformOutput", false), pord, 1);
+ else
+ s_diffp = cell2struct (num2cell (diffp), pord, 1);
+ s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
+ s_orig_lbound = cell2struct (num2cell (lbound), pord, 1);
+ s_orig_ubound = cell2struct (num2cell (ubound), pord, 1);
+ s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
+ s_orig_fixed = cell2struct (num2cell (fixed), pord, 1);
+ endif
+ endif
+
+ #### some further values and checks
+
+ if (any (fixed & (pin < lbound | pin > ubound)))
+ warning ("some fixed parameters outside bounds");
+ endif
+
+ ## dimensions of linear constraints
+ if (isempty (mc))
+ mc = zeros (np, 0);
+ vc = zeros (0, 1);
+ endif
+ if (isempty (emc))
+ emc = zeros (np, 0);
+ evc = zeros (0, 1);
+ endif
+ [rm, cm] = size (mc);
+ [rv, cv] = size (vc);
+ if (rm != np || cm != rv || cv != 1)
+ error ("linear inequality constraints: wrong dimensions");
+ endif
+ [erm, ecm] = size (emc);
+ [erv, ecv] = size (evc);
+ if (erm != np || ecm != erv || ecv != 1)
+ error ("linear equality constraints: wrong dimensions");
+ endif
+
+ ## note initial values of linear constraits
+ pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
+ pin_cstr.equ.lin = emc.' * pin + evc;
+
+ ## note number and initial values of general constraints
+ if (isempty (f_genicstr))
+ pin_cstr.inequ.gen = [];
+ n_genicstr = 0;
+ else
+ n_genicstr = length (pin_cstr.inequ.gen = f_genicstr (pin));
+ endif
+ if (isempty (f_genecstr))
+ pin_cstr.equ.gen = [];
+ n_genecstr = 0;
+ else
+ n_genecstr = length (pin_cstr.equ.gen = f_genecstr (pin));
+ endif
+
+ #### collect remaining settings
+ hook.TolFun = optimget (settings, "TolFun", stol_default);
+ hook.MaxIter = optimget (settings, "MaxIter");
+ if (ischar (hook.cpiv = optimget (settings, "cpiv", @ cpiv_bard)))
+ hook.cpiv = str2func (hook.cpiv);
+ endif
+ hook.Display = optimget (settings, "Display", "off");
+ hook.testing = optimget (settings, "debug", false);
+ hook.new_s = optimget (settings, "lm_svd_feasible_alt_s", false);
+ hook.siman.T_init = optimget (settings, "T_init", .01);
+ hook.siman.T_min = optimget (settings, "T_min", 1.0e-5);
+ hook.siman.mu_T = optimget (settings, "mu_T", 1.005);
+ hook.siman.iters_fixed_T = optimget (settings, "iters_fixed_T", 10);
+ hook.stoch_regain_constr = \
+ optimget (settings, "stoch_regain_constr", false);
+ hook.trace_steps = \
+ optimget (settings, "trace_steps", false);
+ hook.siman_log = \
+ optimget (settings, "siman_log", false);
+ backend = optimget (settings, "Algorithm");
+ if (isempty (backend))
+ error ("At the moment, there is no default algorithm, it must be specified.");
+ endif
+ backend = map_matlab_algorithm_names (backend);
+ backend = map_backend (backend);
+
+ #### handle fixing of parameters
+ orig_lbound = lbound;
+ orig_ubound = ubound;
+ orig_fixed = fixed;
+ if (all (fixed))
+ error ("no free parameters");
+ endif
+
+ nonfixed = ! fixed;
+ if (any (fixed))
+ ## backend (returned values and initial parameters)
+ backend = @ (f, pin, hook) \
+ backend_wrapper (backend, fixed, f, pin, hook);
+
+ ## objective function
+ f = @ (p, varargin) f (assign (pin, nonfixed, p), varargin{:});
+
+ ## gradient of objective function
+ dfdp = @ (p, hook) \
+ dfdp (assign (pin, nonfixed, p), hook)(nonfixed);
+
+ ## hessian of objective function
+ if (! isempty (hessian))
+ hessian = @ (p) \
+ hessian (assign (pin, nonfixed, p))(nonfixed, nonfixed);
+ endif
+
+ ## function for general inequality constraints
+ f_genicstr = @ (p, varargin) \
+ f_genicstr (assign (pin, nonfixed, p), varargin{:});
+
+ ## jacobian of general inequality constraints
+ df_gencstr = @ (p, func, idx, hook) \
+ df_gencstr (assign (pin, nonfixed, p), func, idx, hook) \
+ (:, nonfixed);
+
+ ## function for general equality constraints
+ f_genecstr = @ (p, varargin) \
+ f_genecstr (assign (pin, nonfixed, p), varargin{:});
+
+ ## jacobian of general equality constraints
+ df_genecstr = @ (p, func, idx, hook) \
+ df_genecstr (assign (pin, nonfixed, p), func, idx, hook) \
+ (:, nonfixed);
+
+ ## linear inequality constraints
+ vc += mc(fixed, :).' * (tp = pin(fixed));
+ mc = mc(nonfixed, :);
+
+ ## linear equality constraints
+ evc += emc(fixed, :).' * tp;
+ emc = emc(nonfixed, :);
+
+ ## _last_ of all, vectors of parameter-related configuration,
+ ## including "fixed" itself
+ lbound = lbound(nonfixed, :);
+ ubound = ubound(nonfixed, :);
+ max_fract_change = max_fract_change(nonfixed);
+ fract_prec = fract_prec(nonfixed);
+ max_rand_step = max_rand_step(nonfixed);
+ fixed = fixed(nonfixed);
+ endif
+
+ #### supplement constants to jacobian functions
+
+ ## gradient of objective function
+ if (dfdp_pstruct)
+ dfdp = @ (p, hook) \
+ dfdp (p, cell2fields \
+ ({s_diffp, s_diff_onesided, s_orig_lbound, \
+ s_orig_ubound, s_plabels, \
+ cell2fields(num2cell(hook.fixed), pord(nonfixed), \
+ 1, s_orig_fixed), cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "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), \
+ cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "plabels", "fixed", "h"}, \
+ 2, hook));
+ endif
+
+ ## jacobian of general inequality constraints
+ if (df_inequc_pstruct)
+ df_gencstr = @ (p, func, idx, hook) \
+ df_gencstr (p, func, idx, cell2fields \
+ ({s_diffp, s_diff_onesided, s_orig_lbound, \
+ s_orig_ubound, s_plabels, \
+ cell2fields(num2cell(hook.fixed), pord(nonfixed), \
+ 1, s_orig_fixed), cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "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), cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "plabels", "fixed", "h"}, \
+ 2, hook));
+ endif
+
+ ## jacobian of general equality constraints
+ if (df_equc_pstruct)
+ df_genecstr = @ (p, func, idx, hook) \
+ df_genecstr (p, func, idx, cell2fields \
+ ({s_diffp, s_diff_onesided, s_orig_lbound, \
+ s_orig_ubound, s_plabels, \
+ cell2fields(num2cell(hook.fixed), pord(nonfixed), \
+ 1, s_orig_fixed), cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "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), cstep}, \
+ {"diffp", "diff_onesided", "lbound", "ubound", \
+ "plabels", "fixed", "h"}, \
+ 2, hook));
+ endif
+
+ #### interfaces to constraints
+
+ ## include bounds into linear inequality constraints
+ tp = eye (sum (nonfixed));
+ lidx = lbound != - Inf;
+ uidx = ubound != Inf;
+ mc = cat (2, tp(:, lidx), - tp(:, uidx), mc);
+ vc = cat (1, - lbound(lidx, 1), ubound(uidx, 1), vc);
+
+ ## concatenate linear inequality and equality constraints
+ mc = cat (2, mc, emc);
+ vc = cat (1, vc, evc);
+ n_lincstr = rows (vc);
+
+ ## concatenate general inequality and equality constraints
+ if (n_genecstr > 0)
+ if (n_genicstr > 0)
+ nidxi = 1 : n_genicstr;
+ nidxe = n_genicstr + 1 : n_genicstr + n_genecstr;
+ f_gencstr = @ (p, idx, varargin) \
+ cat (1, \
+ f_genicstr (p, idx(nidxi), varargin{:}), \
+ f_genecstr (p, idx(nidxe), varargin{:}));
+ df_gencstr = @ (p, idx, hook) \
+ cat (1, \
+ df_gencstr (p, @ (p, varargin) \
+ possibly_pstruct_f_genicstr \
+ (p, idx(nidxi), varargin{:}), \
+ idx(nidxi), \
+ setfield (hook, "f", \
+ hook.f(nidxi(idx(nidxi))))), \
+ df_genecstr (p, @ (p, varargin) \
+ possibly_pstruct_f_genecstr \
+ (p, idx(nidxe), varargin{:}), \
+ idx(nidxe), \
+ setfield (hook, "f", \
+ hook.f(nidxe(idx(nidxe))))));
+ else
+ f_gencstr = f_genecstr;
+ df_gencstr = @ (p, idx, hook) \
+ df_genecstr (p, \
+ @ (p, varargin) \
+ possibly_pstruct_f_genecstr \
+ (p, idx, varargin{:}), \
+ idx, \
+ setfield (hook, "f", hook.f(idx)));
+ endif
+ else
+ f_gencstr = f_genicstr;
+ df_gencstr = @ (p, idx, hook) \
+ df_gencstr (p, \
+ @ (p, varargin) \
+ possibly_pstruct_f_genicstr (p, idx, varargin{:}), \
+ idx, \
+ setfield (hook, "f", hook.f(idx)));
+ endif
+ n_gencstr = n_genicstr + n_genecstr;
+
+ ## concatenate linear and general constraints, defining the final
+ ## function interfaces
+ if (n_gencstr > 0)
+ nidxl = 1:n_lincstr;
+ nidxh = n_lincstr + 1 : n_lincstr + n_gencstr;
+ f_cstr = @ (p, idx, varargin) \
+ cat (1, \
+ mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), \
+ f_gencstr (p, idx(nidxh), varargin{:}));
+ df_cstr = @ (p, idx, hook) \
+ cat (1, \
+ mc(:, idx(nidxl)).', \
+ df_gencstr (p, idx(nidxh), \
+ setfield (hook, "f", \
+ hook.f(nidxh))));
+ else
+ f_cstr = @ (p, idx, varargin) mc(:, idx).' * p + vc(idx, 1);
+ df_cstr = @ (p, idx, hook) mc(:, idx).';
+ endif
+
+ ## define eq_idx (logical index of equality constraints within all
+ ## concatenated constraints
+ eq_idx = false (n_lincstr + n_gencstr, 1);
+ eq_idx(n_lincstr + 1 - rows (evc) : n_lincstr) = true;
+ n_cstr = n_lincstr + n_gencstr;
+ eq_idx(n_cstr + 1 - n_genecstr : n_cstr) = true;
+
+ #### prepare interface hook
+
+ ## passed constraints
+ hook.mc = mc;
+ hook.vc = vc;
+ hook.f_cstr = f_cstr;
+ hook.df_cstr = df_cstr;
+ hook.n_gencstr = n_gencstr;
+ hook.eq_idx = eq_idx;
+ hook.lbound = lbound;
+ hook.ubound = ubound;
+
+ ## passed values of constraints for initial parameters
+ hook.pin_cstr = pin_cstr;
+
+ ## passed function for gradient of objective function
+ hook.dfdp = dfdp;
+
+ ## passed function for hessian of objective function
+ hook.hessian = hessian;
+
+ ## passed function for complementary pivoting
+ ## hook.cpiv = cpiv; # set before
+
+ ## passed value of residual function for initial parameters
+ hook.f_pin = f_pin;
+
+ ## passed options
+ hook.max_fract_change = max_fract_change;
+ hook.fract_prec = fract_prec;
+ ## hook.TolFun = ; # set before
+ ## hook.MaxIter = ; # set before
+ hook.fixed = fixed;
+ hook.max_rand_step = max_rand_step;
+
+ #### call backend
+
+ [p, objf, cvg, outp] = backend (f, pin, hook);
+
+ if (pin_struct)
+ if (pnonscalar)
+ p = cell2struct \
+ (cellfun (@ reshape, mat2cell (p, ppartidx), \
+ pdims, "UniformOutput", false), \
+ pord, 1);
+ else
+ p = cell2struct (num2cell (p), pord, 1);
+ endif
+ endif
+
+endfunction
+
+function backend = map_matlab_algorithm_names (backend)
+
+ ## nothing done here at the moment
+
+endfunction
+
+function backend = map_backend (backend)
+
+ switch (backend)
+ case "bfgs_infeasible"
+ backend = "__sqp__";
+ case "sqp"
+ backend = "__sqp__";
+ case "bfgs_feasible"
+ backend = "__bfgs_feasible__";
+ case "siman"
+ backend = "__siman__";
+ otherwise
+ error ("no backend implemented for algorithm '%s'", backend);
+ endswitch
+
+ backend = str2func (backend);
+
+endfunction
+
+function [p, resid, cvg, outp] = backend_wrapper (backend, fixed, f, p, hook)
+
+ [tp, resid, cvg, outp] = backend (f, p(! fixed), hook);
+
+ p(! fixed) = tp;
+
+endfunction
+
+function lval = assign (lval, lidx, rval)
+
+ lval(lidx) = rval;
+
+endfunction
+
+function m = hessian_struct2mat (s, pord)
+
+ m = cell2mat (fields2cell \
+ (structcat (1, NA, fields2cell (s, pord){:}), pord));
+
+ idx = isna (m);
+
+ m(idx) = (m.')(idx);
+
+endfunction
+
+function ret = __optimget__ (s, name, default)
+
+ if (isfield (s, name))
+ ret = s.(name);
+ elseif (nargin > 2)
+ ret = default;
+ else
+ ret = [];
+ endif
+
+endfunction
+
+%!demo
+%! ## Example for simulated annealing, two parameters, "trace_steps"
+%! ## is true;
+%! t_init = .2;
+%! t_min = .002;
+%! mu_t = 1.002;
+%! iters_fixed_t = 10;
+%! init_p = [2; 2];
+%! max_rand_step = [.2; .2];
+%! [p, objf, cvg, outp] = nonlin_min (@ (p) (p(1)/10)^2 + (p(2)/10)^2 + .1 * (-cos(4*p(1)) - cos(4*p(2))), init_p, optimset ("algorithm", "siman", "max_rand_step", max_rand_step, "t_init", t_init, "T_min", t_min, "mu_t", mu_t, "iters_fixed_T", iters_fixed_t, "trace_steps", true));
+%! p
+%! objf
+%! x = (outp.trace(:, 1) - 1) * iters_fixed_t + outp.trace(:, 2);
+%! x(1) = 0;
+%! plot (x, cat (2, outp.trace(:, 3:end), t_init ./ (mu_t .^ outp.trace(:, 1))))
+%! legend ({"objective function value", "p(1)", "p(2)", "Temperature"})
+%! xlabel ("subiteration")
Added: trunk/octave-forge/main/optim/inst/private/__s2mat__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__s2mat__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__s2mat__.m 2012-05-24 06:03:47 UTC (rev 10512)
@@ -0,0 +1,50 @@
+## Copyright (C) 2010 Olaf Till <ola...@un...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## __s2mat__ (s, ord)
+##
+## Returns a matrix of second derivatives with respect to some
+## parameters from the structure-based representation of such a matrix
+## in s, using the order of parameter names ord. s has to contain all
+## fields named in ord. Each field has some subfields named in ord so
+## that each second derivative is represented at least in one of its two
+## possible orders. If it is represented differently in both orders, no
+## error is returned, but both entries might get into the final matrix
+## at symmetric positions.
+##
+## Should be included as a subfunction of a wrapper for optimization
+## functions possibly needing a Hessian.
+
+function ret = __s2mat__ (s, ord)
+
+ if (any (size (s) != [1, 1]))
+ error ("structure must be scalar");
+ endif
+
+ if (! (iscell (ord) && isvector (ord)))
+ error ("ord must be a one-dimensional cell-array");
+ endif
+
+ c = fields2cell (structcat (1, fields2cell (s, ord){:}), ord);
+
+ neidx = ! (eidx = cellfun ("isempty", c));
+
+ ret = zeros (length (ord));
+
+ ret(neidx) = [c{neidx}]; # faster than [c{:}] ?
+
+ ret(eidx) = ret.'(eidx);
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/optim/inst/private/__siman__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__siman__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__siman__.m 2012-05-24 06:03:47 UTC (rev 10512)
@@ -0,0 +1,249 @@
+## Copyright (C) 2012 Olaf Till <i7...@t-...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## The simulated annealing code is translated and adapted from siman.c,
+## written by Mark Galassi, of the GNU Scientific Library.
+
+function [p_res, objf, cvg, outp] = __siman__ (f, pin, hook)
+
+ ## needed for some anonymous functions
+ if (exist ("ifelse") != 5)
+ ifelse = @ scalar_ifelse;
+ endif
+
+ ## passed constraints
+ mc = hook.mc; # matrix of linear constraints
+ vc = hook.vc; # vector of linear constraints
+ f_cstr = hook.f_cstr; # function of all constraints
+ df_cstr = hook.df_cstr; # function of derivatives of all constraints
+ n_gencstr = hook.n_gencstr; # number of non-linear constraints
+ eq_idx = hook.eq_idx; # logical index of equality constraints in all
+ # constraints
+ lbound = hook.lbound; # bounds, subset of linear inequality
+ ubound = hook.ubound; # constraints in mc and vc
+
+ ## passed values of constraints for initial parameters
+ pin_cstr = hook.pin_cstr;
+
+ ## passed return value of f for initial parameters
+ f_pin = hook.f_pin;
+
+ ## passed function for complementary pivoting, currently sqp is used
+ ## instead
+ ##
+ ## cpiv = hook.cpiv;
+
+ ## passed simulated annealing parameters
+ T_init = hook.siman.T_init;
+ T_min = hook.siman.T_min;
+ mu_T = hook.siman.mu_T;
+ iters_fixed_T = hook.siman.iters_fixed_T;
+ max_rand_step = hook.max_rand_step;
+
+ ## passed options
+ fixed = hook.fixed;
+ verbose = strcmp (hook.Display, "iter");
+ regain_constraints = hook.stoch_regain_constr;
+ if ((siman_log = hook.siman_log))
+ log = zeros (0, 5);
+ endif
+ if ((trace_steps = hook.trace_steps))
+ trace = [0, 0, f_pin, pin.'];
+ endif
+
+ ## some useful variables derived from passed variables
+ n = length (pin);
+ sqp_hessian = 2 * eye (n);
+ n_lconstr = length (vc);
+ n_bounds = sum (lbound != -Inf) + sum (ubound != Inf);
+ bidx = false (n_lconstr + n_gencstr, 1);
+ bidx(1 : n_bounds) = true;
+ ac_idx = true (n_lconstr + n_gencstr, 1);
+ ineq_idx = ! eq_idx;
+ leq_idx = eq_idx(1:n_lconstr);
+ lineq_idx = ineq_idx(1:n_lconstr);
+ lfalse_idx = false(n_lconstr, 1);
+
+ nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
+
+ ## backend-specific checking of options and constraints
+ ##
+ ## equality constraints can not be met by chance
+ if ((any (eq_idx) || any (lbound == ubound)) && ! regain_constraints)
+ error ("If 'stoch_regain_constr' is not set, equality constraints or identical lower and upper bounds are not allowed by simulated annealing backend.");
+ endif
+ ##
+ if (any (pin < lbound | pin > ubound) ||
+ any (pin_cstr.inequ.lin_except_bounds < 0) ||
+ any (pin_cstr.inequ.gen < 0) ||
+ any (abs (pin_cstr.equ.lin)) >= nz ||
+ any (abs (pin_cstr.equ.gen)) >= nz)
+ error ("Initial parameters violate constraints.");
+ endif
+ ##
+ if (all (fixed))
+ error ("no free parameters");
+ endif
+ ##
+ idx = isna (max_rand_step);
+ max_rand_step(idx) = 0.005 * pin(idx);
+
+ ## fill constant fields of hook for derivative-functions; some fields
+ ## may be backend-specific
+ dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
+ # the backend still may add to it
+
+ ## set up for iterations
+ sizep = size (pin);
+ p = best_p = pin;
+ E = best_E = f_pin;
+ T = T_init;
+ n_evals = 1; # one has been done by frontend
+ n_iter = 0;
+ done = false;
+
+ cvg = 1;
+
+ ## simulated annealing
+ while (! done)
+
+ n_iter++;
+
+ n_accepts = n_rejects = n_eless = 0;
+
+ for id = 1 : iters_fixed_T
+
+ new_p = p + max_rand_step .* (2 * rand (sizep) - 1);
+
+ ## apply constraints
+ if (regain_constraints)
+ evidx = (abs ((ac = f_cstr (new_p, ac_idx))(eq_idx)) >= nz);
+ ividx = (ac(ineq_idx) < 0);
+ if (any (evidx) || any (ividx))
+ nv = sum (evidx) + sum (ividx);
+ if (sum (lbvidx = (new_p < lbound)) + \
+ sum (ubvidx = (new_p > ubound)) == \
+ nv)
+ ## special case only bounds violated, set back to bound
+ new_p(lbvidx) = lbound(lbvidx);
+ new_p(ubvidx) = ubound(ubvidx);
+ elseif (nv == 1 && \
+ sum (t_eq = (abs (ac(leq_idx)) >= nz)) + \
+ sum (t_inequ = (ac(lineq_idx) < 0)) == 1)
+ ## special case only one linear constraint violated, set
+ ## back perpendicularly to constraint
+ tidx = lfalse_idx;
+ tidx(leq_idx) = t_eq;
+ tidx(lineq_idx) = t_inequ;
+ c = mc(:, tidx);
+ d = ac(tidx);
+ new_p -= c * (d / (c.' * c));
+ else
+ ## other cases, set back keeping the distance to original
+ ## 'new_p' minimal, using quadratic programming, or
+ ## sequential quadratic programming for nonlinear
+ ## constraints
+ [new_p, discarded, sqp_info] = \
+ sqp (new_p, \
+ {@(x)sumsq(x-new_p), \
+ @(x)2*(x-new_p), \
+ @(x)sqp_hessian}, \
+ {@(x)f_cstr(x,eq_idx), \
+ @(x)df_cstr(x,eq_idx, \
+ setfield(hook,"f", \
+ f_cstr(x,ac_idx)))}, \
+ {@(x)f_cstr(x,ineq_idx), \
+ @(x)df_cstr(x,ineq_idx, \
+ setfield(hook,"f", \
+ f_cstr(x,ac_idx)))});
+ if (sqp_info != 101)
+ cvg = 0;
+ done = true;
+ break;
+ endif
+ endif
+ endif
+ else
+ n_retry_constr = 0;
+ while (any (abs ((ac = f_cstr (new_p, ac_idx))(eq_idx)) >= nz) \
+ || any (ac(ineq_idx) < 0))
+ new_p = p + max_rand_step .* (2 * rand (sizep) - 1);
+ n_retry_constr++;
+ endwhile
+ if (verbose && n_retry_constr)
+ printf ("%i additional tries of random step to meet constraints\n",
+ n_retry_constr);
+ endif
+ endif
+
+ new_E = f (new_p);
+ n_evals++;
+
+ if (new_E < best_E)
+ best_p = new_p;
+ best_E = new_E;
+ endif
+ if (new_E < E)
+ ## take a step
+ p = new_p;
+ E = new_E;
+ n_eless++;
+ if (trace_steps)
+ trace(end + 1, :) = [n_iter, id, E, p.'];
+ endif
+ elseif (rand (1) < exp (- (new_E - E) / T))
+ ## take a step
+ p = new_p;
+ E = new_E;
+ n_accepts++;
+ if (trace_steps)
+ trace(end + 1, :) = [n_iter, id, E, p.'];
+ endif
+ else
+ n_rejects++;
+ endif
+
+ endfor # iters_fixed_T
+
+ if (verbose)
+ printf ("temperature no. %i: %e, energy %e,\n", n_iter, T, E);
+ printf ("tries with energy less / not less but accepted / rejected:\n");
+ printf ("%i / %i / %i\n...
[truncated message content] |
|
From: <i7...@us...> - 2012-05-24 06:51:10
|
Revision: 10513
http://octave.svn.sourceforge.net/octave/?rev=10513&view=rev
Author: i7tiol
Date: 2012-05-24 06:50:59 +0000 (Thu, 24 May 2012)
Log Message:
-----------
Updated NEWS and nonlin_min help.
Modified Paths:
--------------
trunk/octave-forge/main/optim/INDEX
trunk/octave-forge/main/optim/NEWS
trunk/octave-forge/main/optim/inst/nonlin_min.m
Modified: trunk/octave-forge/main/optim/INDEX
===================================================================
--- trunk/octave-forge/main/optim/INDEX 2012-05-24 06:03:47 UTC (rev 10512)
+++ trunk/octave-forge/main/optim/INDEX 2012-05-24 06:50:59 UTC (rev 10513)
@@ -12,6 +12,7 @@
fminsearch
cg_min
de_min
+ nonlin_min
Data fitting
expfit
wpolyfit
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-05-24 06:03:47 UTC (rev 10512)
+++ trunk/octave-forge/main/optim/NEWS 2012-05-24 06:50:59 UTC (rev 10513)
@@ -3,7 +3,7 @@
** The following functions are new optim 1.1.0:
- powell cauchy
+ powell cauchy nonlin_min
** The following functions have been deprecated since they have been
implemented in Octave core:
@@ -23,6 +23,23 @@
deriv linprog
- ** The function `line_min' has a configurable setpesize and max evals
+ ** The function `line_min' has a configurable setpesize and max evals.
+ ** Added possibility to restrict a parameter to samin.
+
** Package is no longer automatically loaded.
+
+
+Some important changes of the last versions of optim before 1.1.0:
+------------------------------------------------------------------
+
+ ** New functions:
+
+ jacobs: complex step derivative approximation
+
+ nonlin_residmin, nonlin_curvefit: Frontends with a general
+ interface for constrained residual-based optimization. They
+ allow a.o. optimization of structure-based named parameters or
+ parameter-arrays. A backend is included, which is derived from
+ the algorithm of leasqr, but enables feasible-path optimization
+ with linear and general constraints.
Modified: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-05-24 06:03:47 UTC (rev 10512)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-05-24 06:50:59 UTC (rev 10513)
@@ -17,13 +17,14 @@
## @deftypefn {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin})
## @deftypefnx {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin}, @var{settings})
##
-## Frontend for nonlinear minimization of a scalar objective function.
-## The functions supplied by the user have a minimal interface; any
-## additionally needed constants can be supplied by wrapping the user
-## functions into anonymous functions.
+## Frontend for constrained nonlinear minimization of a scalar objective
+## function. The functions supplied by the user have a minimal
+## interface; any additionally needed constants can be supplied by
+## wrapping the user functions into anonymous functions.
##
## At the moment there is only a stochastic backend implemented, so that
-## some settings in the frontend have no use.
+## some settings in the frontend have no use. A backend for
+## feasible-path sequential quadratic programming is being worked on.
##
## The following description applies to usage with vector-based
## parameter handling. Differences in usage for structure-based
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-06-07 09:33:31
|
Revision: 10577
http://octave.svn.sourceforge.net/octave/?rev=10577&view=rev
Author: i7tiol
Date: 2012-06-07 09:33:21 +0000 (Thu, 07 Jun 2012)
Log Message:
-----------
Pass linker flags with better compatibility.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/inst/nonlin_min.m
trunk/octave-forge/main/optim/src/Makefile
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-06-07 08:27:45 UTC (rev 10576)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-06-07 09:33:21 UTC (rev 10577)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.1.0
-Date: 2012-05-24
+Version: 1.1.1
+Date: 2012-06-07
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Modified: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-07 08:27:45 UTC (rev 10576)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-07 09:33:21 UTC (rev 10577)
@@ -1283,12 +1283,12 @@
function backend = map_backend (backend)
switch (backend)
- case "bfgs_infeasible"
+ case "sqp_infeasible"
backend = "__sqp__";
case "sqp"
backend = "__sqp__";
- case "bfgs_feasible"
- backend = "__bfgs_feasible__";
+ case "sqp_feasible"
+ backend = "__sqp_feasible__";
case "siman"
backend = "__siman__";
otherwise
Modified: trunk/octave-forge/main/optim/src/Makefile
===================================================================
--- trunk/octave-forge/main/optim/src/Makefile 2012-06-07 08:27:45 UTC (rev 10576)
+++ trunk/octave-forge/main/optim/src/Makefile 2012-06-07 09:33:21 UTC (rev 10577)
@@ -5,10 +5,17 @@
# <car...@gm...>
BLAS_LIBS := $(shell $(MKOCTFILE) -p BLAS_LIBS)
+# Passing LFLAGS, supplemented with LAPACK_LIBS and BLAS_LIBS, in the
+# environment to mkoctfile is prefered over passing LAPACK_LIBS and
+# BLAS_LIBS in mkoctfiles commandline due to mkoctfiles difficulties
+# with non-standard flags on some systems (e.g. -framework ... on
+# Apple)
+LFLAGS := $(shell $(MKOCTFILE) -p LFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS)
+
all: __bfgsmin.oct numgradient.oct numhessian.oct samin.oct __disna_optim__.oct
%.oct: %.cc
- $(MKOCTFILE) -s $< $(LAPACK_LIBS) $(BLAS_LIBS)
+ LFLAGS="$(LFLAGS)" $(MKOCTFILE) -s $<
clean:
$(RM) *.o core octave-core *.oct *~
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-06-12 14:09:04
|
Revision: 10622
http://octave.svn.sourceforge.net/octave/?rev=10622&view=rev
Author: i7tiol
Date: 2012-06-12 14:08:57 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Completed optimization suite with new backend lm_feasible.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/NEWS
trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
trunk/octave-forge/main/optim/inst/nonlin_min.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.1.1
-Date: 2012-06-07
+Version: 1.2.0
+Date: 2012-06-12
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/NEWS 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,3 +1,26 @@
+Summary of important user-visible changes for optim 1.2.0:
+-------------------------------------------------------------------
+
+ ** Together with the new backend "lm_feasible" there is now a
+ complete suite of backends for optimization with linear and
+ general equality and inequality constraints, for scalar valued
+ objective functions and for array valued model function, which
+ features, a.o., honouring of constraints throughout optimization
+ and handling of structure-based parameters. The respective user
+ functions (frontends) are
+
+ nonlin_min nonlin_residmin nonlin_curvefit
+
+ together with a user function for statistics
+
+ residmin_stat
+
+ ** The requirement of nonlin_min, nonlin_residmin, and
+ nonlin_curvefit for the general constraint functions being able to
+ honour an index of constraints has been removed, the respective
+ feature is still available by setting some options.
+
+
Summary of important user-visible changes for optim 1.1.0:
-------------------------------------------------------------------
Modified: trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -81,3 +81,21 @@
endif
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_residmin')
+%!
+%! ## independents and observations
+%! indep = 1:5;
+%! obs = [1, 2, 4, 7, 14];
+%! ## model function:
+%! f = @ (p, x) p(1) * exp (p(2) * x);
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, model_values, cvg, outp] = nonlin_curvefit (f, init, indep, obs, settings)
Modified: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -22,10 +22,6 @@
## interface; any additionally needed constants can be supplied by
## wrapping the user functions into anonymous functions.
##
-## At the moment there is only a stochastic backend implemented, so that
-## some settings in the frontend have no use. A backend for
-## feasible-path sequential quadratic programming is being worked on.
-##
## The following description applies to usage with vector-based
## parameter handling. Differences in usage for structure-based
## parameter handling will be explained in a separate section below.
@@ -52,14 +48,14 @@
## currently can be @code{0} (maximum number of iterations exceeded),
## @code{1} (fixed number of iterations completed, e.g. in stochastic
## optimizers), @code{2} (parameter change less than specified precision
-## in two consecutive iterations), or @code{3} (improvement in objective
-## function less than specified).
+## in two consecutive iterations), @code{3} (improvement in objective
+## function less than specified), or @code{-4} (algorithm got stuck).
##
## @var{settings}:
##
-## @code{Algorithm}: String specifying the backend. Currently there is
-## no default, since only a stochastic backend is implemented
-## (@code{"siman"}). It is described in a separate section below.
+## @code{Algorithm}: String specifying the backend. Currently available
+## are @{"lm_feasible"} (default) and @code{"siman"}. They are described
+## in separate sections below.
##
## @code{objf_grad}: Function computing the gradient of the objective
## function with respect to the parameters, assuming residuals are
@@ -138,24 +134,29 @@
## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
## paraters and the optional argument @code{idx} is a logical index.
## @code{h} has to return the values of all constraints if @code{idx} is
-## not given, and has to return only the indexed constraints if
+## not given. It may choose to return only the indexed constraints if
## @code{idx} is given (so computation of the other constraints can be
-## spared). In gradient determination, this function may be called with
-## an informational third argument, whose content depends on the
-## function for gradient determination. If a second entry for general
-## inequality constraints is given, it must be a function computing the
-## jacobian of the constraints with respect to the parameters. For this
-## function, the description of @code{dfdp} above applies, except that
-## it is called with 3 arguments since it has an additional argument
-## @code{idx} --- a logical index --- at second position, indicating
-## which rows of the jacobian must be returned, and except that the
-## default function calls @code{h} with 3 arguments, since the argument
-## @code{idx} is also supplied. Note that specifying linear constraints
-## as general constraints will generally waste performance, even if
-## further, non-linear, general constraints are also specified.
+## spared); in this case, the additional setting @code{inequc_f_idx} has
+## to be set to @code{true}. In gradient determination, this function
+## may be called with an informational third argument, whose content
+## depends on the function for gradient determination. If a second entry
+## for general inequality constraints is given, it must be a function
+## computing the jacobian of the constraints with respect to the
+## parameters. For this function, the description of @code{dfdp} above
+## applies, with 2 exceptions: 1) it is called with 3 arguments since it
+## has an additional argument @code{idx} --- a logical index --- at
+## second position, indicating which rows of the jacobian must be
+## returned (if the function chooses to return only indexed rows, the
+## additional setting @code{inequc_df_idx} has to be set to
+## @code{true}). 2) the default jacobian function calls @code{h} with 3
+## arguments, since the argument @code{idx} is also supplied. Note that
+## specifying linear constraints as general constraints will generally
+## waste performance, even if further, non-linear, general constraints
+## are also specified.
##
## @code{equc}: Equality constraints. Specified the same way as
-## inequality constraints (see @code{inequc}).
+## inequality constraints (see @code{inequc}). The respective additional
+## settings are named @code{equc_f_idx} and @code{equc_df_idx}.
##
## @code{cpiv}: Function for complementary pivoting, usable in
## algorithms for constraints. Default: @ cpiv_bard. Only the default
@@ -242,8 +243,30 @@
## parameters must contain arrays with dimensions reshapable to those of
## the respective parameters.
##
-## Description of backends (currently only one)
+## Description of backends
##
+## "lm_feasible"
+##
+## A Levenberg/Marquardt-like optimizer, attempting to honour
+## constraints throughout the course of optimization. This means that
+## the initial parameters must not violate constraints (to find an
+## initial feasible set of parameters, e.g. Octaves @code{sqp} can be
+## used, by specifying an objective function which is constant or which
+## returns the quadratic distance to the initial values). If the
+## constraints need only be honoured in the result of the optimization,
+## Octaves @code{sqp} may be preferable. The Hessian is either supplied
+## by the user or is approximated by the BFGS algorithm.
+##
+## Returned value @var{cvg} will be @code{2} or @code{3} for success and
+## @code{0} or @code{-4} for failure (see above for meaning).
+##
+## Backend-specific defaults are: @code{MaxIter}: 20, @code{fract_prec}:
+## @code{zeros (size (parameters))}, @code{max_fract_change}: @code{Inf}
+## for all parameters.
+##
+## Interpretation of @code{Display}: if set to @code{"iter"}, currently
+## only information on applying @code{max_fract_change} is printed.
+##
## "siman"
##
## A simulated annealing (stochastic) optimizer, changing all parameters
@@ -335,10 +358,14 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"TolFun", stol_default, \
"MaxIter", [], \
"Display", "off", \
- "Algorithm", [], \ # set reasonable default, once available
+ "Algorithm", "lm_feasible", \
"T_init", .01, \
"T_min", 1.0e-5, \
"mu_T", 1.005, \
@@ -351,10 +378,12 @@
return;
endif
- assign = @ assign; # Is this faster in repeated calls?
+ if (nargin < 2 || nargin > 3)
+ print_usage ();
+ endif
- if (nargin != 3)
- error ("incorrect number of arguments");
+ if (nargin == 2)
+ settings = struct ();
endif
if (ischar (f))
@@ -418,13 +447,36 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
if (isempty (hessian)) hessian_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -988,7 +1040,6 @@
endif
hook.Display = optimget (settings, "Display", "off");
hook.testing = optimget (settings, "debug", false);
- hook.new_s = optimget (settings, "lm_svd_feasible_alt_s", false);
hook.siman.T_init = optimget (settings, "T_init", .01);
hook.siman.T_min = optimget (settings, "T_min", 1.0e-5);
hook.siman.mu_T = optimget (settings, "mu_T", 1.005);
@@ -999,10 +1050,7 @@
optimget (settings, "trace_steps", false);
hook.siman_log = \
optimget (settings, "siman_log", false);
- backend = optimget (settings, "Algorithm");
- if (isempty (backend))
- error ("At the moment, there is no default algorithm, it must be specified.");
- endif
+ backend = optimget (settings, "Algorithm", "lm_feasible");
backend = map_matlab_algorithm_names (backend);
backend = map_backend (backend);
@@ -1246,7 +1294,7 @@
## passed function for complementary pivoting
## hook.cpiv = cpiv; # set before
- ## passed value of residual function for initial parameters
+ ## passed value of objective function for initial parameters
hook.f_pin = f_pin;
## passed options
@@ -1283,12 +1331,12 @@
function backend = map_backend (backend)
switch (backend)
- case "sqp_infeasible"
- backend = "__sqp__";
- case "sqp"
- backend = "__sqp__";
- case "sqp_feasible"
- backend = "__sqp_feasible__";
+ ## case "sqp_infeasible"
+ ## backend = "__sqp__";
+ ## case "sqp"
+ ## backend = "__sqp__";
+ case "lm_feasible"
+ backend = "__lm_feasible__";
case "siman"
backend = "__siman__";
otherwise
@@ -1336,7 +1384,24 @@
endfunction
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
+
%!demo
+%! ## Example for default optimization (Levenberg/Marquardt with
+%! ## BFGS), one non-linear equality constraint. Constrained optimum is
+%! ## at p = [0; 1].
+%! objective_function = @ (p) p(1)^2 + p(2)^2;
+%! pin = [-2; 5];
+%! constraint_function = @ (p) p(1)^2 + 1 - p(2);
+%! [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset ("equc", {constraint_function}))
+
+%!demo
%! ## Example for simulated annealing, two parameters, "trace_steps"
%! ## is true;
%! t_init = .2;
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -129,21 +129,25 @@
## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
## paraters and the optional argument @code{idx} is a logical index.
## @code{h} has to return the values of all constraints if @code{idx} is
-## not given, and has to return only the indexed constraints if
+## not given. It may choose to return only the indexed constraints if
## @code{idx} is given (so computation of the other constraints can be
-## spared). In gradient determination, this function may be called with
-## an informational third argument, whose content depends on the
-## function for gradient determination. If a second entry for general
-## inequality constraints is given, it must be a function computing the
-## jacobian of the constraints with respect to the parameters. For this
-## function, the description of @code{dfdp} above applies, except that
-## it is called with 3 arguments since it has an additional argument
-## @code{idx} --- a logical index --- at second position, indicating
-## which rows of the jacobian must be returned, and except that the
-## default function calls @code{h} with 3 arguments, since the argument
-## @code{idx} is also supplied. Note that specifying linear constraints
-## as general constraints will generally waste performance, even if
-## further, non-linear, general constraints are also specified.
+## spared); in this case, the additional setting @code{inequc_f_idx} has
+## to be set to @code{true}. In gradient determination, this function
+## may be called with an informational third argument, whose content
+## depends on the function for gradient determination. If a second entry
+## for general inequality constraints is given, it must be a function
+## computing the jacobian of the constraints with respect to the
+## parameters. For this function, the description of @code{dfdp} above
+## applies, with 2 exceptions: 1) it is called with 3 arguments since it
+## has an additional argument @code{idx} --- a logical index --- at
+## second position, indicating which rows of the jacobian must be
+## returned (if the function chooses to return only indexed rows, the
+## additional setting @code{inequc_df_idx} has to be set to
+## @code{true}). 2) the default jacobian function calls @code{h} with 3
+## arguments, since the argument @code{idx} is also supplied. Note that
+## specifying linear constraints as general constraints will generally
+## waste performance, even if further, non-linear, general constraints
+## are also specified.
##
## @code{equc}: Equality constraints. Specified the same way as
## inequality constraints (see @code{inequc}).
@@ -281,3 +285,20 @@
[p, resid, cvg, outp] = __nonlin_residmin__ (varargin{:});
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_curvefit')
+%!
+%! ## independents
+%! indep = 1:5;
+%! ## residual function:
+%! f = @ (p) p(1) * exp (p(2) * indep) - [1, 2, 4, 7, 14];
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, residuals, cvg, outp] = nonlin_residmin (f, init, settings)
Added: trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -0,0 +1,505 @@
+## Copyright (C) 2012 Olaf Till <i7...@t-...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+function [p_res, objf, cvg, outp] = __lm_feasible__ (f, pin, hook)
+
+ ## some backend specific defaults
+ fract_prec_default = 0;
+ max_fract_step_default = Inf;
+
+ ## needed for some anonymous functions
+ if (exist ("ifelse") != 5)
+ ifelse = @ scalar_ifelse;
+ endif
+
+ n = length (pin);
+
+ ## passed constraints
+ mc = hook.mc; # matrix of linear constraints
+ vc = hook.vc; # vector of linear constraints
+ f_cstr = hook.f_cstr; # function of all constraints
+ df_cstr = hook.df_cstr; # function of derivatives of all constraints
+ n_gencstr = hook.n_gencstr; # number of non-linear constraints
+ eq_idx = hook.eq_idx; # logical index of equality constraints in all
+ # constraints
+ lbound = hook.lbound; # bounds, subset of linear inequality
+ ubound = hook.ubound; # constraints in mc and vc
+
+ ## passed values of constraints for initial parameters
+ pin_cstr = hook.pin_cstr;
+
+ ## passed return value of f for initial parameters
+ f_pin = hook.f_pin;
+
+ ## passed function for gradient of objective function
+ grad_f = hook.dfdp;
+
+ ## passed function for hessian of objective function
+ if (isempty (hessian = hook.hessian))
+ user_hessian = false;
+ A = eye (n);
+ else
+ user_hessian = true;
+ endif
+
+ ## passed function for complementary pivoting
+ cpiv = hook.cpiv;
+
+ ## passed options
+ ftol = hook.TolFun;
+ if (isempty (niter = hook.MaxIter)) niter = 20; endif
+ fixed = hook.fixed;
+ maxstep = hook.max_fract_change;
+ maxstep(isna (maxstep)) = max_fract_step_default;
+ pprec = hook.fract_prec;
+ pprec(isna (pprec)) = fract_prec_default;
+ verbose = strcmp (hook.Display, "iter");
+
+ ## some useful variables derived from passed variables
+ n_lcstr = size (vc, 1);
+ have_constraints_except_bounds = \
+ n_lcstr + n_gencstr > \
+ sum (lbound != -Inf) + sum (ubound != Inf);
+ ac_idx = true (n_lcstr + n_gencstr, 1); # index of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); # none of all constraints
+ gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); # gen. constr.
+
+ nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
+
+ ## backend-specific checking of options and constraints
+ ##
+ if (any (pin < lbound | pin > ubound) ||
+ any (pin_cstr.inequ.lin_except_bounds < 0) ||
+ any (pin_cstr.inequ.gen < 0) ||
+ any (abs (pin_cstr.equ.lin)) >= nz ||
+ any (abs (pin_cstr.equ.gen)) >= nz)
+ error ("Initial parameters violate constraints.");
+ endif
+ ##
+ idx = lbound == ubound;
+ if (any (idx))
+ warning ("lower and upper bounds identical for some parameters, fixing the respective parameters");
+ fixed(idx) = true;
+ endif
+ if (all (fixed))
+ error ("no free parameters");
+ endif
+ if (n_gencstr > 0 && any (! isinf (maxstep)))
+ warning ("setting both a maximum fractional step change of parameters and general constraints may result in inefficiency and failure");
+ endif
+
+ ## fill constant fields of hook for derivative-functions; some fields
+ ## may be backend-specific
+ dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
+ # the backend still may add to it
+
+ ## set up for iterations
+ p = pbest = pin;
+ vf = fbest = f_pin;
+ iter = 0;
+ done = false;
+ ll = 1;
+ ltab = [.1, 1, 1e2, 1e4, 1e6];
+ chgprev = Inf (n, 1);
+ df = [];
+ c_act = false (n, 1);
+ dca = zeros (n, 0);
+
+ while (! done)
+
+ iter++;
+
+ ## gradient of objective function
+ old_df = df;
+ df = grad_f (p, setfield (dfdp_hook, "f", f))(:);
+
+ ## constraints, preparation of some constants
+ v_cstr = f_cstr (p, ac_idx);
+ old_c_act = c_act;
+ old_dca = dca;
+ c_act = v_cstr < nz | eq_idx; # index of active constraints
+ if (any (c_act))
+
+ if (n_gencstr)
+ ## full gradient is needed later
+ dct = df_cstr (p, ac_idx, setfield (dfdp_hook, "f", v_cstr));
+ dct(:, fixed) = 0; # for user supplied dfdp; necessary?
+ dcat = dct(c_act, :);
+ else
+ dcat = df_cstr (p, c_act, setfield (dfdp_hook, "f", v_cstr));
+ dcat(:, fixed) = 0; # for user supplied dfdp; necessary?
+ endif
+
+ dca = dcat.';
+
+ a_eq_idx = eq_idx(c_act);
+
+ else
+
+ dca = zeros (n, 0);
+
+ endif
+
+ ## hessian of objectiv function
+ if (user_hessian)
+
+ A = hessian (p);
+ idx = isnan (A);
+ A(idx) = A.'(idx);
+ if (any (isnan (A(:))))
+ error ("some second derivatives undefined by user function");
+ endif
+ if (! isreal (A))
+ error ("second derivatives given by user function not real");
+ endif
+ if (! issymmetric (A))
+ error ("Hessian returned by user function not symmetric");
+ endif
+
+ elseif (iter > 1)
+
+ if (any (chg))
+
+ ## approximate Hessian of Lagrangian
+
+ ## I wonder if this hassle here and above with accounting for
+ ## changing active sets is indeed better than just approximating
+ ## the Hessian only of the objective function.
+ ##
+ ## index, over all constraints, of constraints active both
+ ## previously and currently
+ s_c_act = old_c_act & c_act;
+ ## index, over currently active constraints, of constraints
+ ## active both previously and currently
+ id_new = s_c_act(c_act);
+ ## index, over previously active constraints, of constraints
+ ## active both previously and currently
+ id_old = s_c_act(old_c_act);
+ ## gradients of currently active constraints which were also
+ ## active previously
+ dca_new_id = dca(:, id_new);
+ ## gradients of previously active constraints which are also
+ ## active currently
+ dca_old_id = old_dca(:, id_old);
+ ## index, over constraints active both previously and currently,
+ ## of (old) non-zero multipliers (bidx set below previously)
+ bidx_old_id = bidx(id_old);
+ ## index, over (old) non-zero multipliers, of constraints active
+ ## both previously and currently (bidx set below previously)
+ old_l_idx = id_old(bidx);
+
+ ## difference of derivatives of new and old active constraints,
+ ## multiplied by multipliers, as used for BFGS update (lb set
+ ## below previously)
+ dch = (dca_new_id(:, bidx_old_id) - \
+ dca_old_id(:, bidx_old_id)) * \
+ lb(old_l_idx);
+
+ y = df - old_df - dch;
+
+ ## Damped BFGS according to Nocedal & Wright, 2nd edition,
+ ## procedure 18.2.
+ chgt = chg.';
+ sAs = chgt * A * chg;
+ cy = chgt * y;
+ if (cy >= .2 * sAs)
+ th = 1;
+ else
+ if ((den1 = sAs - cy) == 0)
+ cvg = -4;
+ break;
+ endif
+ th = .8 * sAs / den1;
+ endif
+ Ac = A * chg;
+ r = th * y + (1 - th) * Ac;
+
+ if ((den2 = chgt * r) == 0 || sAs == 0)
+ cvg = -4;
+ break;
+ endif
+ A += r * r.' / den2 - Ac * Ac.' / sAs;
+
+ endif
+
+ endif
+
+ ## Inverse scaled decomposition A = G * (1 ./ L) * G.'
+ ##
+ ## make matrix Binv for scaling
+ Binv = diag (A);
+ nidx = ! (idx = Binv == 0);
+ Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
+ Binv(idx) = 1;
+ Binv = diag (Binv);
+ ## eigendecomposition of scaled A
+ [V, L] = eig (Binv * A * Binv);
+ L = diag (L);
+ ## A is symmetric, so V and L are real, delete any imaginary parts,
+ ## which might occur due to inaccuracy
+ V = real (V);
+ L = real (L);
+ ##
+ nminL = - min (L) * 1.1 / ltab(1);
+ G = Binv * V;
+
+ ## Levenberg/Marquardt
+ fgoal = (1 - ftol) * vf;
+ for l = ltab
+
+ ll = max (ll, nminL);
+ l = max (1e-7, ll * l);
+
+ R = G * diag (1 ./ (L + l)) * G.';
+
+ ## step computation
+ if (any (c_act))
+
+ ## some constraints are active, quadratic programming
+
+ tp = dcat * R;
+ [lb, bidx, ridx, tbl] = cpiv (- tp * df, tp * dca, a_eq_idx);
+ chg = R * (dca(:, bidx) * lb - df); # step direction
+
+ ## indices for different types of constraints
+ c_inact = ! c_act; # inactive constraints
+ c_binding = c_unbinding = nc_idx;
+ c_binding(c_act) = bidx; # constraints selected binding
+ c_unbinding(c_act) = ridx; # constraints unselected binding
+ c_nonbinding = c_act & ! (c_binding | c_unbinding); #
+ #constraints selected non-binding
+
+ else
+
+ ## no constraints are active, chg is the Levenberg/Marquardt step
+
+ chg = - R * df; # step direction
+
+ lb = zeros (0, 1);
+ bidx = false (0, 1);
+
+ ## indices for different types of constraints (meaning see above)
+ c_inact = ac_idx;
+ c_binding = nc_idx;
+ c_unbinding = nc_idx;
+ c_nonbinding = nc_idx;
+
+ endif
+
+ ## apply inactive and non-binding constraints to step width
+ ##
+ ## linear constraints
+ k = 1;
+ c_tp = c_inact(1:n_lcstr);
+ mcit = mc(:, c_tp).';
+ vci = vc(c_tp);
+ hstep = mcit * chg;
+ idx = hstep < 0;
+ if (any (idx))
+ k = min (1, min (- (vci(idx) + mcit(idx, :) * p) ./ \
+ hstep(idx)));
+ endif
+ ##
+ ## general constraints
+ if (n_gencstr)
+ c_tp = gc_idx & (c_nonbinding | c_inact);
+ if (any (c_tp) && any (f_cstr (p + k * chg, c_tp) < 0))
+ [k, fval, info] = \
+ fzero (@ (x) min (cat (1, \
+ f_cstr (p + x * chg, c_tp), \
+ k - x, \
+ ifelse (x < 0, -Inf, Inf))), \
+ 0);
+ if (info != 1 || abs (fval) >= nz)
+ error ("could not find stepwidth to satisfy inactive and non-binding general inequality constraints");
+ endif
+ endif
+ endif
+ ##
+ chg = k * chg;
+
+ ## if necessary, regain binding constraints and one of the
+ ## possibly active previously inactive or non-binding constraints
+ if (any (gc_idx & c_binding)) # none selected binding => none
+ # unselected binding
+ ptp1 = p + chg;
+
+ tp = true;
+ nt_nosuc = true;
+ lim = 20;
+ while (nt_nosuc && lim >= 0)
+ ## we keep d_p.' * inv (R) * d_p minimal in each step of the
+ ## inner loop
+ c_tp0 = c_inact | c_nonbinding;
+ c_tp1 = c_inact | (gc_idx & c_nonbinding);
+ btbl = tbl(bidx, bidx);
+ c_tp2 = c_binding;
+ ## once (any(tp)==false), it would not get true again even
+ ## with the following assignment
+ if (any (tp) && \
+ any (tp = f_cstr (ptp1, c_tp1) < nz))
+ ## keep only the first true entry in tp
+ tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
+ ## supplement binding index with one (the first) getting
+ ## binding in c_tp1
+ c_tp2(c_tp1) = tp;
+ ## gradient of this added constraint
+ caddt = dct(c_tp2 & ! c_binding, :);
+ cadd = caddt.';
+ C = dct(c_binding, :) * R * cadd;
+ Ct = C.';
+ T = [btbl, btbl * C; \
+ -Ct * btbl, caddt * R * cadd - Ct * btbl * C];
+ btbl = gjp (T, size (T, 1));
+ endif
+ dcbt = dct(c_tp2, :);
+ mfc = - R * dcbt.' * btbl;
+
+ ptp2 = ptp1;
+ nt_niter = nt_niter_start = 100;
+ while (nt_nosuc && nt_niter >= 0)
+ hv = f_cstr (ptp2, c_tp2);
+ if (all (abs (hv) < nz))
+ nt_nosuc = false;
+ chg = ptp2 - p;
+ else
+ ptp2 = ptp2 + mfc * hv; # step should be zero for each
+ # component for which the parameter is
+ # "fixed"
+ endif
+ nt_niter--;
+ endwhile
+
+ if (nt_nosuc || \
+ any (abs (chg) > abs (p .* maxstep)) || \
+ any (f_cstr (ptp2, c_tp0) < -nz))
+ ## if (nt_nosuc), regaining did not converge, else,
+ ## regaining violated type 3 and 4.
+ nt_nosuc = true;
+ ptp1 = (p + ptp1) / 2;
+ endif
+ if (! nt_nosuc && \
+ any ((tp = f_cstr (ptp2, c_unbinding)) < 0))
+ [discarded, id] = min(tp);
+ tid = find (ridx);
+ id = tid(id); # index within active constraints
+ unsuccessful_exchange = false;
+ if (abs (tbl(id, id)) < nz) # Bard: not absolute value
+ ## exchange this unselected binding constraint against a
+ ## binding constraint, but not against an equality
+ ## constraint
+ tbidx = bidx & ! a_eq_idx;
+ if (! any (tbidx))
+ unsuccessful_exchange = true;
+ else
+ [discarded, idm] = max (abs (tbl(tbidx, id)));
+ tid = find (tbidx);
+ idm = tid(idm); # -> index within active constraints
+ tbl = gjp (tbl, idm);
+ bidx(idm) = false;
+ ridx(idm) = true;
+ endif
+ endif
+ if (unsuccessful_exchange)
+ ## It probably doesn't look good now; this desperate last
+ ## attempt is not in the original algortithm, since that
+ ## didn't account for equality constraints.
+ ptp1 = (p + ptp1) / 2;
+ else
+ tbl = gjp (tbl, id);
+ bidx(id) = true;
+ ridx(id) = false;
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx;
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx;
+ endif
+ ## regaining violated type 2 constraints
+ nt_nosuc = true;
+ endif
+ lim--;
+ endwhile
+ if (nt_nosuc)
+ error ("could not regain binding constraints");
+ endif
+ else
+ ## check the maximal stepwidth and apply as necessary
+ ochg = chg;
+ idx = ! isinf (maxstep);
+ limit = abs (maxstep(idx) .* p(idx));
+ chg(idx) = min (max (chg(idx), - limit), limit);
+ if (verbose && any (ochg != chg))
+ printf ("Change in parameter(s): %s:maximal fractional stepwidth enforced", \
+ sprintf ("%d ", find (ochg != chg)));
+ endif
+ endif # regaining
+
+ aprec = abs (pprec .* pbest);
+ if (any (abs (chg) > 0.1 * aprec)) # only worth evaluating
+ # function if there is some
+ # non-miniscule change
+ p_chg = p + chg;
+ ## since the projection method may have slightly violated
+ ## constraints due to inaccuracy, correct parameters to bounds
+ ## --- but only if no further constraints are given, otherwise
+ ## the inaccuracy in honoring them might increase by this
+ if (! have_constraints_except_bounds)
+ lidx = p_chg < lbound;
+ uidx = p_chg > ubound;
+ p_chg(lidx, 1) = lbound(lidx, 1);
+ p_chg(uidx, 1) = ubound(uidx, 1);
+ chg(lidx, 1) = p_chg(lidx, 1) - p(lidx, 1);
+ chg(uidx, 1) = p_chg(uidx, 1) - p(uidx, 1);
+ endif
+ ##
+ if (! isreal (vf_chg = f (p_chg)))
+ error ("objective function not real");
+ endif
+ if (vf_chg < fbest)
+ pbest = p_chg;
+ fbest = vf_chg;
+ endif
+ if (vf_chg <= fgoal)
+ p = p_chg;
+ vf = vf_chg;
+ break;
+ endif
+ endif
+ endfor
+
+ ll = l;
+
+ aprec = abs (pprec .* pbest);
+ if (vf_chg < eps || vf_chg > fgoal)
+ cvg = 3;
+ done = true;
+ elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
+ cvg = 2;
+ done = true;
+ elseif (iter == niter)
+ cvg = 0;
+ done = true;
+ else
+ chgprev = chg;
+ endif
+
+ endwhile
+
+ ## return result
+ p_res = pbest;
+ objf = fbest;
+ outp.niter = iter;
+
+endfunction
Modified: trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -145,7 +145,7 @@
epsLlast=1;
epstab=[.1, 1, 1e2, 1e4, 1e6];
ac_idx = true (n_lcstr + n_gencstr, 1); % all constraints
- nc_idx = false (n_lcstr + n_gencstr, 1); % non of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); % none of all constraints
gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); % gen. constr.
lc_idx = ~gc_idx;
@@ -159,10 +159,10 @@
# violated at start
if (any (c_act))
if (n_gencstr > 0)
+ %% full gradient is needed later
dct = df_cstr (p, ac_idx, ...
setfield (dfdp_hook, 'f', v_cstr));
dct(:, fixed) = 0; % for user supplied dfdp; necessary?
- dc = dct.';
dcat = dct(c_act, :);
else
dcat = df_cstr (p, c_act, ...
@@ -420,6 +420,12 @@
if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating
% function if there is some non-miniscule
% change
+ %% In the code of the outer loop before the inner loop pbest is
+ %% actually identical to p, since once they deviate, the outer
+ %% loop will not be repeated. Though the inner loop can still be
+ %% repeated in this case, pbest is not used in it. Since pprev
+ %% is set from pbest in the outer loop before the inner loop, it
+ %% is also identical to p up to here.
p=chg+pprev;
%% since the projection method may have slightly violated
%% constraints due to inaccuracy, correct parameters to bounds
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -69,6 +69,10 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"weights", [], \
"TolFun", stol_default, \
"MaxIter", [], \
@@ -142,12 +146,35 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -995,3 +1022,11 @@
endif
endfunction
+
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-09-14 21:34:00
|
Revision: 11024
http://octave.svn.sourceforge.net/octave/?rev=11024&view=rev
Author: i7tiol
Date: 2012-09-14 21:33:53 +0000 (Fri, 14 Sep 2012)
Log Message:
-----------
Fixed wrong test for existence of 'options' in leasqr. New release.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/NEWS
trunk/octave-forge/main/optim/inst/leasqr.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-09-14 19:36:44 UTC (rev 11023)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-09-14 21:33:53 UTC (rev 11024)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.2.0
-Date: 2012-06-12
+Version: 1.2.1
+Date: 2012-09-14
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-09-14 19:36:44 UTC (rev 11023)
+++ trunk/octave-forge/main/optim/NEWS 2012-09-14 21:33:53 UTC (rev 11024)
@@ -1,3 +1,9 @@
+optim 1.2.1:
+------------
+
+ ** Bugfix in leasqr.m: errors when a function 'options' is in
+ namespace.
+
Summary of important user-visible changes for optim 1.2.0:
-------------------------------------------------------------------
Modified: trunk/octave-forge/main/optim/inst/leasqr.m
===================================================================
--- trunk/octave-forge/main/optim/inst/leasqr.m 2012-09-14 19:36:44 UTC (rev 11023)
+++ trunk/octave-forge/main/optim/inst/leasqr.m 2012-09-14 21:33:53 UTC (rev 11024)
@@ -570,7 +570,7 @@
%% only preliminary, for testing
hook.testing = false;
hook.new_s = false;
- if (exist ('options'))
+ if (nargin > 9)
if (isfield (options, 'testing'))
hook.testing = options.testing;
end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X