Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / optim-1.2.0 / jacobs.m

diff --git a/octave_packages/optim-1.2.0/jacobs.m b/octave_packages/optim-1.2.0/jacobs.m
new file mode 100644 (file)
index 0000000..521b259
--- /dev/null
+++ b/octave_packages/optim-1.2.0/jacobs.m
@@ -0,0 +1,142 @@
+## Copyright (C) 2011 Fotios Kasolis <fotios.kasolis@gmail.com>
+##
+## 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} {Df =} jacobs (@var{x}, @var{f})
+## @deftypefnx {Function File} {Df =} jacobs (@var{x}, @var{f}, @var{hook})
+## Calculate the jacobian of a function using the complex step method.
+##
+## Let @var{f} be a user-supplied function. Given a point @var{x} at
+## which we seek for the Jacobian, the function @command{jacobs} returns
+## the Jacobian matrix @code{d(f(1), @dots{}, df(end))/d(x(1), @dots{},
+## x(n))}. The function uses the complex step method and thus can be
+## applied to real analytic functions.
+##
+## The optional argument @var{hook} is a structure with additional options. @var{hook}
+## can have the following fields:
+## @itemize @bullet
+## @item
+## @code{h} - can be used to define the magnitude of the complex step and defaults
+## to 1e-20; steps larger than 1e-3 are not allowed.
+## @item
+## @code{fixed} - is a logical vector internally usable by some optimization
+## functions; it indicates for which elements of @var{x} no gradient should be
+## computed, but zero should be returned.
+## @end itemize
+##
+## For example:
+## 
+## @example
+## @group
+## f = @@(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
+## Df = jacobs ([1, 2], f)
+## @end group
+## @end example
+## @end deftypefn
+
+function Df = jacobs (x, f, hook)
+
+  if ( (nargin < 2) || (nargin > 3) )
+    print_usage ();
+  endif
+
+  if (ischar (f))
+    f = str2func (f, "global");
+  endif
+
+  n  = numel (x);
+
+  default_h = 1e-20;
+  max_h = 1e-3; # must be positive
+
+  if (nargin > 2)
+
+    if (isfield (hook, "h"))
+      if (! (isscalar (hook.h)))
+        error ("complex step magnitude must be a scalar");
+      endif
+      if (abs (hook.h) > max_h)
+        warning ("complex step magnitude larger than allowed, set to %e", ...
+                 max_h)
+        h = max_h;
+      else
+        h = hook.h;
+      endif
+    else
+      h = default_h;
+    endif
+
+    if (isfield (hook, "fixed"))
+      if (numel (hook.fixed) != n)
+        error ("index of fixed parameters has wrong dimensions");
+      endif
+      fixed = hook.fixed;
+    else
+      fixed = false (n, 1);
+    endif
+
+  else
+
+    h = default_h;
+    fixed = false (n, 1);
+
+  endif
+
+  if (all (fixed))
+    error ("all elements of 'x' are fixed");
+  endif
+
+  x = repmat (x(:), 1, n) + h * 1i * eye (n);
+
+  idx = find (! fixed);
+
+  ## after first evaluation, dimensionness of 'f' is known
+  t_Df = imag (f (x(:, idx(1)))(:));
+  dim = numel (t_Df);
+
+  Df = zeros (dim, n);
+
+  Df(:, idx(1)) = t_Df;
+
+  for count = idx(2:end)
+    Df(:, count) = imag (f (x(:, count))(:));
+  endfor
+
+  Df /=  h;
+
+endfunction
+
+%!assert (jacobs (1, @(x) x), 1)
+%!assert (jacobs (6, @(x) x^2), 12)
+%!assert (jacobs ([1; 1], @(x) [x(1)^2; x(1)*x(2)]), [2, 0; 1, 1])
+%!assert (jacobs ([1; 2], @(x) [x(1)^2 + x(2); x(2)*exp(x(1))]), [2, 1; 2*exp(1), exp(1)])
+
+%% Test input validation
+%!error jacobs ()
+%!error jacobs (1)
+%!error jacobs (1, 2, 3, 4)
+%!error jacobs (@sin, 1, [1, 1])
+%!error jacobs (@sin, 1, ones(2, 2))
+
+%!demo
+%! # Relative error against several h-values
+%! k = 3:20; h = 10 .^ (-k); x = 0.3*pi;
+%! err = zeros (1, numel (k));
+%! for count = 1 : numel (k)
+%!   err(count) = abs (jacobs (x, @sin,        struct ("h", h(count))) - cos (x)) / abs (cos (x)) + eps;
+%! endfor
+%! loglog (h, err); grid minor;
+%! xlabel ("h"); ylabel ("|Df(x) - cos(x)| / |cos(x)|")
+%! title ("f(x)=sin(x), f'(x)=cos(x) at x = 0.3pi")