1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2008 - Rainer von Seggern
3 // Copyright (C) 2008 - Bruno Pincon
4 // Copyright (C) 2009 - INRIA - Michael Baudin
5 // Copyright (C) 2010 - DIGITEO - Michael Baudin
6 //
7 // This file must be used under the terms of the CeCILL.
8 // This source file is licensed as described in the file COPYING, which
9 // you should have received as part of this distribution.  The terms
10 // are also available at
11 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
13 //
14 //  PURPOSE
15 //     First and second order numerical derivatives of a function  F: R^n --> R^m
16 //     by finite differences.
17 //     J=J(x) is the m x n Jacobian (the gradient for m=1), and H=H(x) the Hessian
18 //     of the m components of F at x. The default form of H is a mxn^2 matrix;
19 //     in this form the Taylor series of F up to second order terms is given by:
20 //
21 //        F(x+dx) = F(x) + J(x)*dx + 1/2*H(x)*(dx.*.dx) +...
22 //
23 //  NOTES
24 //     1/ See derivative.cat for details of the parameters
25 //
26 //     2/ This function uses the 3 "internal" functions (following
27 //        this one in this file) :
28 //
29 //          %DF_      => used to compute the Hessian by "differentiating
30 //                       the derivative"
31 //          %deriv1_  => contains the various finite difference formulae
32 //          %R_       => to deal with F as this arg may be a scilab
33 //                       function or a list embedding a function with
34 //                       its parameters
35 //
36 function [J,H] = derivative(F, x, h, order, H_form, Q , verbose )
37     warnobsolete("numderivative","6.0")
38     [lhs,rhs]=argn();
39     if rhs<2 | rhs>6 then
40         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"derivative",2,6));
41     end
42     if type(x) ~= 1 then
43         error(msprintf(gettext("%s: Wrong type for input argument #%d: N-dimensional array expected.\n"),"derivative",2));
44     end
45     [n,p] = size(x)
46     if p ~= 1 then
47         error(msprintf(gettext("%s: Wrong size for input argument #%d: A column vector expected.\n"),"derivative",2));
48     end
49     if ~exists("order","local") then
50         order = 2
51     elseif (order ~= 1 & order ~= 2 & order ~= 4) then
52         error(msprintf(gettext("%s: Order must be 1, 2 or 4.\n"),"derivative"));
53     end
55     if ~exists("H_form","local") then
56         H_form = "default"
57     end
59     if ~exists("Q","local") then
60         Q = eye(n,n);
61     else
62         if norm(clean(Q*Q'-eye(n,n)))>0 then
63             error(msprintf(gettext("%s: Q must be orthogonal.\n"),"derivative"));
64         end
65     end
66     if ~exists("h","local") then
67         h_not_given = %t
68         // stepsizes for approximation of first derivatives
69         select order
70         case 1
71             h = sqrt(%eps)
72         case 2
73             h = %eps^(1/3)
74         case 4
75             // TODO : check this, i found 1/5 instead !
76             h = %eps^(1/4)
77         end
78     else
79         h_not_given = %f
80     end
82     if ~exists("verbose","local") then
83         verbose = 0
84     end
86     if verbose == 1 then
87         mprintf("h = %s\n", string(h))
88         mprintf("order = %d\n", order)
89         mprintf("H_form = %s\n", H_form)
90         mprintf("Q = \n")
91         for i = 1:n
92             mprintf("%s\n", strcat(string(Q(i,:)), " "))
93         end
94     end
96     J = %deriv1_(F, x, h, order, Q)
97     m = size(J,1);
99     if lhs == 1 then
100         return
101     end
103     if h_not_given then
104         // stepsizes for approximation of second derivatives
105         select order
106         case 1
107             h = %eps^(1/3)
108         case 2
109             h = %eps^(1/4)
110         case 4
111             h = %eps^(1/6)
112         end
113     end
114     // H is a mxn^2 block matrix
115     H = %deriv1_(%DF_, x, h, order, Q)
116     if     H_form == "default" then
117         // H has the old scilab form
118         H = matrix(H',n*n,m)'
119     end
120     if H_form == "hypermat" then
121         if m>1 then
122             // H is a hypermatrix if m>1
123             H=H';
124             H=hypermat([n n m],H(:));
125         end
126     end
127     if (H_form ~= "blockmat")&(H_form ~= "default")&(H_form ~= "hypermat") then
128         error(msprintf(gettext("%s: H_form must be ""default"",""blockmat"" or ""hypermat"", but current H_form=%s\n"),"derivative",H_form));
129     end
130 endfunction
132 function z=%DF_(x)
133     // Transpose !
134     z = %deriv1_(F, x, h, order, Q)';
135     z = z(:);
136 endfunction
138 function g=%deriv1_(F_, x, h, order, Q)
139     n=size(x,"*")
140     Dy=[];
141     select order
142     case 1
143         D = h*Q;
144         y=%R_(F_,x);
145         for d=D
146             Dy=[Dy %R_(F_,x+d)-y]
147         end
148         g=Dy*Q'/h
149     case 2
150         D = h*Q;
151         for d=D
152             Dy=[Dy %R_(F_,x+d)-%R_(F_,x-d)]
153         end
154         g=Dy*Q'/(2*h)
155     case 4
156         D = h*Q;
157         for d=D
158             dFh =  (%R_(F_,x+d)-%R_(F_,x-d))/(2*h)
159             dF2h = (%R_(F_,x+2*d)-%R_(F_,x-2*d))/(4*h)
160             Dy=[Dy (4*dFh - dF2h)/3]
161         end
162         g = Dy*Q'
163     end
164 endfunction
166 function y=%R_(F_,x)
167     if type(F_)==15 then
168         if ( length(F_) < 2 ) then
169             error(msprintf(gettext("%s: Wrong number of elements in input argument #%d: At least %d elements expected, but current number is %d.\n"),"derivative",1,2,length(F_)));
170         end
171         R=F_(1);
172         if ( and(typeof(R) <> ["function" "funptr"]) ) then
173             error(msprintf(gettext("%s: Wrong type for element #%d in input argument #%d: A function is expected, but current type is %s.\n"),"derivative",1,1,typeof(R)));
174         end
175         y=R(x,F_(2:\$));
176         // See extraction, list or tlist case: ...
177         // But if the extraction syntax is used within a function
178         // input calling sequence each returned list component is
179         // added to the function calling sequence.
180     elseif  type(F_)==13 | type(F_)==11 then
181         y=F_(x);
182     else
183         error(msprintf(gettext("%s: Wrong type for input argument #%d: A function expected.\n"),"derivative",1));
184     end
185 endfunction