1## Copyright (C) 2011 Fotios Kasolis <fotios.kasolis@gmail.com>
2## Copyright (C) 2013-2019 Olaf Till <i7tiol@t-online.de>
3##
4## This program is free software; you can redistribute it and/or modify it under
5## the terms of the GNU General Public License as published by the Free Software
6## Foundation; either version 3 of the License, or (at your option) any later
7## version.
8##
9## This program is distributed in the hope that it will be useful, but WITHOUT
10## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12## details.
13##
14## You should have received a copy of the GNU General Public License along with
15## this program; if not, see <http://www.gnu.org/licenses/>.
16
17## -*- texinfo -*-
18## @deftypefn {Function File} {Df =} jacobs (@var{x}, @var{f})
19## @deftypefnx {Function File} {Df =} jacobs (@var{x}, @var{f}, @var{hook})
20## Calculate the jacobian of a function using the complex step method.
21##
22## Let @var{f} be a user-supplied function. Given a point @var{x} at
23## which we seek for the Jacobian, the function @command{jacobs} returns
24## the Jacobian matrix @code{d(f(1), @dots{}, df(end))/d(x(1), @dots{},
25## x(n))}. The function uses the complex step method and thus can be
26## applied to real analytic functions.
27##
28## The optional argument @var{hook} is a structure with additional options. @var{hook}
29## can have the following fields:
30## @itemize @bullet
31## @item
32## @code{h} - can be used to define the magnitude of the complex step and defaults
33## to 1e-20; steps larger than 1e-3 are not allowed.
34## @item
35## @code{fixed} - is a logical vector internally usable by some optimization
36## functions; it indicates for which elements of @var{x} no gradient should be
37## computed, but zero should be returned.
38## @end itemize
39##
40## For example:
41##
42## @example
43## @group
44## f = @@(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
45## Df = jacobs ([1, 2], f)
46## @end group
47## @end example
48## @end deftypefn
49
50function Df = jacobs (x, f, hook)
51
52  if ( (nargin < 2) || (nargin > 3) )
53    print_usage ();
54  endif
55
56  if (ischar (f))
57    f = str2func (f, "global");
58  endif
59
60  n  = numel (x);
61
62  parclass = class (x);
63
64  default_h = 1e-20 * eps (parclass) / eps;
65  max_h = 1e-3; # must be positive
66
67  if (nargin > 2)
68
69    if (isfield (hook, "h"))
70      if (! (isscalar (hook.h)))
71        error ("complex step magnitude must be a scalar");
72      endif
73      if (abs (hook.h) > max_h)
74        warning ("complex step magnitude larger than allowed, set to %e", ...
75                 max_h)
76        h = max_h;
77      else
78        h = hook.h;
79      endif
80    else
81      h = default_h;
82    endif
83
84    if (isfield (hook, "fixed"))
85      if (numel (hook.fixed) != n)
86        error ("index of fixed parameters has wrong dimensions");
87      endif
88      fixed = hook.fixed(:);
89    else
90      fixed = false (n, 1);
91    endif
92
93  else
94    h = default_h;
95    fixed = false (n, 1);
96  endif
97
98  if (all (fixed))
99    error ("all elements of 'x' are fixed");
100  endif
101
102  x = repmat (x(:), 1, n) + h * 1i * eye (n);
103
104  idx = find (! fixed);
105
106  if (nargin > 2)
107
108    if (isfield (hook, 'parallel_local'))
109      parallel_local = hook.parallel_local;
110    else
111      parallel_local = false;
112    end
113
114    if (isfield (hook, "parallel_net"))
115      parallel_net = hook.parallel_net;
116    else
117      parallel_net = [];
118    endif
119
120    if (parallel_local || ! isempty (parallel_net))
121
122      parallel = true;
123      ## user function
124      func = @ (id) {imag(f(x(:, id))(:)) / h, false, []}{:};
125      ## error handler
126      errh = @ (s, id) {[], true, s}{:};
127
128      if (parallel_local && ! isempty (parallel_net))
129        error ("If option 'parallel_net' is not empty, option 'parallel_local' must not be true.");
130      endif
131
132      if (parallel_local)
133
134        if (parallel_local > 1)
135          npr = parallel_local;
136        else
137          npr = nproc ("current");
138        endif
139
140        parfun = @ () pararrayfun (npr, func, idx,
141                                   "UniformOutput", false,
142                                   "VerboseLevel", 0,
143                                   "ErrorHandler", errh);
144
145      else # ! isempty (parallel_net)
146        parfun = @ () netarrayfun (parallel_net, func, idx,
147                                   "UniformOutput", false,
148                                   "ErrorHandler", errh);
149      endif
150
151    else
152      parallel = false;
153    endif
154
155  else
156    parallel = false;
157  endif
158
159  if (parallel)
160
161    [t_Df, err, info] = parfun ();
162
163    ## check for errors
164    if (any ((err = [err{:}])))
165      id = find (err, 1);
166      error ("Some subprocesses, calling model function for complex step derivatives, returned and error. Message of first of these (with id %i): %s%s",
167             id, info{id}.message, print_stack (info{id}));
168    endif
169
170    ## process output
171    t_Df = horzcat (t_Df{:});
172    Df = zeros (rows (t_Df), n, parclass);
173    Df(:, idx) = t_Df;
174
175  else # not parallel
176
177    ## after first evaluation, dimensionness of 'f' is known
178    t_Df = imag (f (x(:, idx(1)))(:));
179    dim = numel (t_Df);
180
181    Df = zeros (dim, n, parclass);
182
183    Df(:, idx(1)) = t_Df;
184
185    for count = (idx.')(2:end)
186      Df(:, count) = imag (f (x(:, count))(:));
187    endfor
188
189    Df /=  h;
190
191  endif
192
193endfunction
194
195function ret = print_stack (info)
196
197  ret = "";
198
199  if (isfield (info, "stack"))
200    for id = 1 : numel (info.stack)
201      ret = cstrcat (ret, sprintf ("\n    %s at line %i comumn %i",
202                                   info.stack(id).name,
203                                   info.stack(id).line,
204                                   info.stack(id).column));
205    endfor
206  endif
207
208endfunction
209
210%!assert (jacobs (1, @(x) x), 1)
211%!assert (jacobs (6, @(x) x^2), 12)
212%!assert (jacobs ([1; 1], @(x) [x(1)^2; x(1)*x(2)]), [2, 0; 1, 1])
213%!assert (jacobs ([1; 2], @(x) [x(1)^2 + x(2); x(2)*exp(x(1))]), [2, 1; 2*exp(1), exp(1)])
214
215%% Test input validation
216%!error jacobs ()
217%!error jacobs (1)
218%!error jacobs (1, 2, 3, 4)
219%!error jacobs (@sin, 1, [1, 1])
220%!error jacobs (@sin, 1, ones(2, 2))
221
222%!demo
223%! # Relative error against several h-values
224%! k = 3:20; h = 10 .^ (-k); x = 0.3*pi;
225%! err = zeros (1, numel (k));
226%! for count = 1 : numel (k)
227%!   err(count) = abs (jacobs (x, @sin,	struct ("h", h(count))) - cos (x)) / abs (cos (x)) + eps;
228%! endfor
229%! loglog (h, err); grid minor;
230%! xlabel ("h"); ylabel ("|Df(x) - cos(x)| / |cos(x)|")
231%! title ("f(x)=sin(x), f'(x)=cos(x) at x = 0.3pi")
232