[Octave-cvsupdate] SF.net SVN: octave:[8374] trunk/octave-forge/main/optim

[<i7...@us...>]

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.