# HG changeset patch
# User address@hidden
# Date 1232655579 -3600
# Node ID c103a79020fb957c57a99b59f9dfd38b8936d83b
# Parent  4a864c4e682de52e9b70b1f9c0886787202326a2
scripts/general/gradient.m: Add support for computing the gradient of a function handle

diff -r 4a864c4e682d -r c103a79020fb scripts/ChangeLog
--- a/scripts/ChangeLog	Thu Jan 22 13:50:08 2009 +0100
+++ b/scripts/ChangeLog	Thu Jan 22 21:19:39 2009 +0100
@@ -1,3 +1,8 @@
+2009-01-22  Soren Hauberg  <address@hidden>
+
+	* general/gradient.m: Add support for computing the gradient of a function
+	  handle.
+
 2009-01-22  Jaroslav Hajek  <address@hidden>
 
 	* optimization/fsolve.m: Undo the last change.
diff -r 4a864c4e682d -r c103a79020fb scripts/general/gradient.m
--- a/scripts/general/gradient.m	Thu Jan 22 13:50:08 2009 +0100
+++ b/scripts/general/gradient.m	Thu Jan 22 21:19:39 2009 +0100
@@ -21,8 +21,12 @@
 ## @deftypefnx {Function File} address@hidden, @var{y}, @dots{}] =} gradient (@var{M})
 ## @deftypefnx {Function File} address@hidden =} gradient (@var{M}, @var{s})
 ## @deftypefnx {Function File} address@hidden =} gradient (@var{M}, @var{dx}, @var{dy}, @dots{})
+## @deftypefnx {Function File} address@hidden =} gradient (@var{F}, @var{x0})
+## @deftypefnx {Function File} address@hidden =} gradient (@var{F}, @var{x0}, @var{s})
+## @deftypefnx {Function File} address@hidden =} gradient (@var{F}, @var{x0}, @var{dx}, @var{dy}, @dots{})
 ##
-## Calculates the gradient. @address@hidden = gradient (@var{M})}
+## Calculates the gradient of sampled data, or of a function.
+## @address@hidden = gradient (@var{M})}
 ## calculates the one dimensional gradient if @var{M} is a vector. If
 ## @var{M} is a matrix the gradient is calculated for each row.
 ##
@@ -46,6 +50,13 @@
 ## y'(i) = 1/(x(i+1)-x(i-1)) *(y(i-1)-y(i+1)).
 ## @end example
 ## 
+## If the first argument @var{F} is a function handle, the gradient of the
+## function at the points in @var{x0} is approximated using central difference.
+## For example, @code{gradient (@@cos, 0)} approximates the gradient of the cosine
+## function in the point @math{x0 = 0}. As with sampled data, the spacing values
+## between the points from which the gradient is estimated can be set via the
+## @var{s} or @var{dx}, @var{dy}, @dots{} arguments. By default a spacing of 1
+## is used.
 ## @end deftypefn
 
 ## Author:  Kai Habel <address@hidden>
@@ -57,6 +68,20 @@
     print_usage ()
   endif
 
+  nargout_with_ans = max(1,nargout);
+  if (ismatrix(M))
+    [varargout{1:nargout_with_ans}] = matrix_gradient(M, varargin{:});
+  elseif (isa(M, "function_handle"))
+    [varargout{1:nargout_with_ans}] = handle_gradient(M, varargin{:});
+  elseif (ischar(M))
+    [varargout{1:nargout_with_ans}] = handle_gradient(str2func(M), varargin{:});
+  else
+    error("gradient: first input must be an array or a function");
+  endif
+
+endfunction
+
+function [varargout] = matrix_gradient(M, varargin)
   transposed = false;
   if (isvector (M))
     ## Make a column vector.
@@ -138,3 +163,75 @@
     varargout{1} = varargout{1}.';
   endif
 endfunction
+
+function [varargout] = handle_gradient (f, p0, varargin)
+  ## Input checking
+  p0_size = size (p0);
+
+  if (numel (p0_size) != 2)
+    error ("gradient: the second input argument should either be a vector or a matrix");
+  endif
+
+  if (any (p0_size == 1))
+    p0 = p0 (:);
+    dim = 1;
+    num_points = numel (p0);
+  else
+    num_points = p0_size (1);
+    dim = p0_size (2);
+  endif
+  
+  if (length (varargin) == 0)
+    delta = 1;
+  elseif (length (varargin) == 1 || length (varargin) == dim)
+    try
+      delta = [varargin{:}];
+    catch
+      error ("gradient: spacing parameters must be scalars or a vector");
+    end_try_catch
+  else
+    error ("gradient: incorrect number of spacing parameters");
+  endif
+  
+  if (isscalar (delta))
+    delta = repmat (delta, 1, dim);
+  elseif (!isvector (delta))
+    error ("gradient: spacing values must be scalars or a vector");
+  endif
+  
+  ## Calculate the gradient
+  p0 = mat2cell (p0, num_points, ones (1, dim));
+  varargout = cell (1,dim);
+  for d = 1:dim
+    s = delta (d);
+    df_dx = (f (p0{1:d-1}, p0{d}+s, p0{d+1:end})
+           - f (p0{1:d-1}, p0{d}-s, p0{d+1:end})) ./ (2*s);
+    if (dim == 1)
+      varargout{d} = reshape (df_dx, p0_size);
+    else
+      varargout{d} = df_dx;
+    endif
+  endfor
+endfunction
+
+%!test
+%! data = [1, 2, 4, 2];
+%! dx = gradient (data);
+%! assert (dx, [1, 3/2, 0, -2]);
+
+%!test
+%! x = 0:10;
+%! f = @cos;
+%! df_dx = @(x) -sin (x);
+%! assert (gradient (f, x), df_dx (x), 0.2);
+%! assert (gradient (f, x, 0.5), df_dx (x), 0.1);
+
+%!test
+%! xy = reshape (1:10, 5, 2);
+%! f = @(x,y) sin (x) .* cos (y);
+%! df_dx = @(x, y) cos (x) .* cos (y);
+%! df_dy = @(x, y) -sin (x) .* sin (y);
+%! [dx, dy] = gradient (f, xy);
+%! assert (dx, df_dx (xy (:, 1), xy (:, 2)), 0.1)
+%! assert (dy, df_dy (xy (:, 1), xy (:, 2)), 0.1)
+