Mercurial > hg > octave-kai
changeset 8601:b297b86f4ad9
gradient.m: Add support for computing the gradient of a function handle
author | sh@sh-t400 |
---|---|
date | Tue, 27 Jan 2009 12:28:22 -0500 |
parents | a6c1aa6f5915 |
children | 827f0285a201 |
files | scripts/ChangeLog scripts/general/gradient.m |
diffstat | 2 files changed, 126 insertions(+), 25 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2009-01-27 Søren Hauberg <hauberg@gmail.com> + + * general/gradient.m: Handle computing the gradient of a function + handle. + 2009-01-27 Jaroslav Hajek <highegg@gmail.com> * optimization/lsqnonneg.m: Reimplement using QR updating for @@ -39,7 +44,7 @@ * sparse/svds.m: svds.m: skip tests if ARPACK is missing. -2009-01-23 Søren Hauberg <hauberg@gmail.com> +2009-01-23 Søren Hauberg <hauberg@gmail.com> * help/type.m: Make 'type X' work, when X is the name of a variable.
--- a/scripts/general/gradient.m +++ b/scripts/general/gradient.m @@ -17,17 +17,20 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{x} =} gradient (@var{M}) -## @deftypefnx {Function File} {[@var{x}, @var{y}, @dots{}] =} gradient (@var{M}) -## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{M}, @var{s}) -## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{M}, @var{dx}, @var{dy}, @dots{}) +## @deftypefn {Function File} {@var{x} =} gradient (@var{m}) +## @deftypefnx {Function File} {[@var{x}, @var{y}, @dots{}] =} gradient (@var{m}) +## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{s}) +## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{dx}, @var{dy}, @dots{}) +## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}) +## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{s}) +## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{dx}, @var{dy}, @dots{}) ## -## Calculates the gradient. @code{@var{x} = 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. +## Calculate the gradient of sampled data, or of a function. If @var{m} +## is a vector, calculate the one dimensional gradient of @var{m}. If +## @var{m} is a matrix the gradient is calculated for each row. ## -## @code{[@var{x}, @var{y}] = gradient (@var{M})} calculates the one -## dimensional gradient for each direction if @var{M} if @var{M} is a +## @code{[@var{x}, @var{y}] = gradient (@var{m})} calculates the one +## dimensional gradient for each direction if @var{m} if @var{m} is a ## matrix. Additional return arguments can be use for multi-dimensional ## matrices. ## @@ -37,7 +40,7 @@ ## values of the spacing can be supplied by the @var{dx}, etc variables. ## A scalar value specifies an equidistant spacing, while a vector value ## can be used to specify a variable spacing. The length must match -## their respective dimension of @var{M}. +## their respective dimension of @var{m}. ## ## At boundary points a linear extrapolation is applied. Interior points ## are calculated with the first approximation of the numerical gradient @@ -46,26 +49,47 @@ ## 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 <kai.habel@gmx.de> ## Modified: David Bateman <dbateman@free.fr> Added NDArray support -function varargout = gradient (M, varargin) +function varargout = gradient (m, varargin) if (nargin < 1) print_usage () endif - transposed = false; - if (isvector (M)) - ## Make a column vector. - transposed = (size (M, 2) == 1); - M = M(:)'; + 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 - nd = ndims (M); - sz = size (M); +endfunction + +function varargout = matrix_gradient (m, varargin) + transposed = false; + if (isvector (m)) + ## make a column vector. + transposed = (size (m, 2) == 1); + m = m(:)'; + endif + + nd = ndims (m); + sz = size (m); if (nargin > 2 && nargin != nd + 1) print_usage () endif @@ -112,21 +136,21 @@ for i = 1:max (2, min (nd, nargout)) mr = sz(i); mc = prod ([sz(1:i-1), sz(i+1:nd)]); - Y = zeros (size (M), class (M)); + Y = zeros (size (m), class (m)); if (mr > 1) ## Top and bottom boundary. - Y(1,:) = diff (M(1:2,:)) / d{i}(1); - Y(mr,:) = diff (M(mr-1:mr,:)) / d{i}(mr-1); + Y(1,:) = diff (m(1:2,:)) / d{i}(1); + Y(mr,:) = diff (m(mr-1:mr,:)) / d{i}(mr-1); endif if (mr > 2) ## Interior points. - Y(2:mr-1,:) = (M(3:mr,:) .- M(1:mr-2,:)) ... - ./ kron (d{i}(1:mr-2) .+ d{i}(2:mr-1), ones (1, mc)); + Y(2:mr-1,:) = ((m(3:mr,:) - m(1:mr-2,:)) + ./ kron (d{i}(1:mr-2) + d{i}(2:mr-1), ones (1, mc))); endif varargout{i} = ipermute (Y, [i:nd,1:i-1]); - M = permute (M, [2:nd,1]); + m = permute (m, [2:nd,1]); endfor ## Why the hell did Matlab decide to swap these two values? @@ -138,3 +162,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) +