matrix high-order differentials
FORMAT [dfdx] = pr_spm_diff(f,x,...,n)
f - [inline] function f(x{1},...)
x - input argument[s]
n - arguments to differentiate w.r.t.
dfdx - df/dx{i} ; n = i
dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p) ; n = [i j ... k]
- a cunning recursive routine
__________________________________________________________________________
Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0001 function [J] = pr_spm_diff(varargin) 0002 % matrix high-order differentials 0003 % FORMAT [dfdx] = pr_spm_diff(f,x,...,n) 0004 % 0005 % f - [inline] function f(x{1},...) 0006 % x - input argument[s] 0007 % n - arguments to differentiate w.r.t. 0008 % 0009 % dfdx - df/dx{i} ; n = i 0010 % dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p) ; n = [i j ... k] 0011 % 0012 % - a cunning recursive routine 0013 %__________________________________________________________________________ 0014 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience 0015 0016 % Karl Friston 0017 % $Id: spm_diff.m 417 2006-02-01 13:50:14Z karl $ 0018 0019 % create inline object 0020 %-------------------------------------------------------------------------- 0021 f = fcnchk(varargin{1}); 0022 x = varargin(2:(end - 1)); 0023 n = varargin{end}; 0024 m = n(end); 0025 xm = pr_spm_vec(x{m}); 0026 dx = exp(-8); 0027 J = cell(1,length(xm)); 0028 0029 % proceed to derivatives 0030 %========================================================================== 0031 if length(n) == 1 0032 0033 % dfdx 0034 %---------------------------------------------------------------------- 0035 f0 = feval(f,x{:}); 0036 for i = 1:length(J) 0037 xi = x; 0038 xmi = xm; 0039 xmi(i) = xmi(i) + dx; 0040 xi{n} = pr_spm_unvec(xmi,xi{n}); 0041 fi = feval(f,xi{:}); 0042 J{i} = pr_spm_dfdx(fi,f0,dx); 0043 end 0044 0045 else 0046 0047 % dfdxdxdx.... 0048 %---------------------------------------------------------------------- 0049 f0 = pr_spm_diff(f,x{:},n(1:end - 1)); 0050 for i = 1:length(J) 0051 xi = x; 0052 xmi = xm; 0053 xmi(i) = xmi(i) + dx; 0054 xi{m} = pr_spm_unvec(xmi,xi{m}); 0055 fi = pr_spm_diff(f,xi{:},n(1:end - 1)); 0056 J{i} = pr_spm_dfdx(fi,f0,dx); 0057 end 0058 return 0059 0060 end 0061 0062 % return numeric array for first order derivatives 0063 %========================================================================== 0064 0065 % vectorise f 0066 %-------------------------------------------------------------------------- 0067 f = pr_spm_vec(f0); 0068 0069 % if there are no arguments to differentiate w.r.t. ... 0070 %-------------------------------------------------------------------------- 0071 if ~length(xm) 0072 J = sparse(length(f),0); 0073 return 0074 end 0075 0076 % if there are no arguments to differentiate 0077 %-------------------------------------------------------------------------- 0078 if ~length(f) 0079 J = sparse(0,length(xm)); 0080 return 0081 end 0082 0083 % if f is a scalar 0084 %-------------------------------------------------------------------------- 0085 if length(f) == 1 0086 J = pr_spm_cat(J); 0087 return 0088 end 0089 0090 % if x{n} is a scalar 0091 %-------------------------------------------------------------------------- 0092 if length(xm) == 1 0093 J = pr_spm_cat(J); 0094 return 0095 end 0096 0097 % else f and xm are vectors return numeric array 0098 %-------------------------------------------------------------------------- 0099 for i = 1:length(J) 0100 J{i} = pr_spm_vec(J{i}); 0101 end 0102 J = pr_spm_cat(J); 0103 return 0104 0105 function dfdx = pr_spm_dfdx(f,f0,dx) 0106 % cell subtraction 0107 %-------------------------------------------------------------------------- 0108 if iscell(f) 0109 dfdx = f; 0110 for i = 1:length(f(:)) 0111 dfdx{i} = pr_spm_dfdx(f{i},f0{i},dx); 0112 end 0113 else 0114 dfdx = (f - f0)/dx; 0115 end