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