14485a8f40dee9e2baf11ae60d9a05dd76a916ff
[scilab.git] / scilab / modules / optimization / macros / numderivative.sci
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-2011 - 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
12
13 function [J, H] = numderivative(varargin)
14     //
15     // Check input arguments
16     [lhs, rhs] = argn();
17     if (rhs < 2 | rhs > 6) then
18         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"), "numderivative", 2, 6));
19     end
20     if (lhs < 1 | lhs > 2) then
21         error(msprintf(gettext("%s: Wrong number of output arguments: %d to %d expected.\n"), "numderivative", 1, 2));
22     end
23     //
24     // Get input arguments
25     __numderivative_f__ = varargin(1)
26     //
27     // Manage x, to get the size n.
28     x = varargin(2);
29     if type(x) ~= 1 then
30         error(msprintf(gettext("%s: Wrong type for argument #%d: Matrix expected.\n"), "numderivative", 2));
31     end
32     [n, p] = size(x);
33     if (n <> 1 & p <> 1) then
34         error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector expected.\n"), "numderivative", 2));
35     end
36     // Make x a column vector, if required
37     if p <> 1 then
38         x = x(:);
39         [n, p] = size(x);
40     end
41     //
42     // Manage h: make it a column vector, if required.
43     h = [];
44     if rhs >= 3 then
45         h = varargin(3);
46     end
47     if type(h) ~= 1 then
48         error(msprintf(gettext("%s: Wrong type for argument #%d: Matrix expected.\n"), "numderivative", 3));
49     end
50     if h <> [] then
51         if size(h, "*") <> 1 then
52             [nrows, ncols] = size(h);
53             if (nrows <> 1 & ncols <> 1) then
54                 error(msprintf(gettext("%s: Wrong size for input argument #%d: Vector expected.\n"), "numderivative", 3));
55             end
56             if ncols <> 1 then
57                 h = h(:);
58             end
59             if or(size(h) <> [n 1]) then
60                 error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"), "numderivative", 3, 1));
61             end
62         end
63     end
64     order = 2;
65     if (rhs >= 4 & varargin(4) <> []) then
66         order = varargin(4);
67     end
68     H_form = "default";
69     if (rhs >= 5 & varargin(5) <> []) then
70         H_form = varargin(5);
71     end
72     Q = eye(n, n);
73     Q_not_given = %t;
74     if (rhs >= 6 & varargin(6) <> []) then
75         Q = varargin(6);
76         Q_not_given = %f;
77     end
78     //
79     // Check types
80     if and(type(__numderivative_f__) <> [11 13 15 130]) then
81         // Must be a function (uncompiled or compiled) or a list
82         error(msprintf(gettext("%s: Wrong type for argument #%d: Function or list expected.\n"), "numderivative", 1));
83     end
84     if type(__numderivative_f__) == 15 then
85         // List case
86         // Check that the first element in the list is a function
87         if and(type(__numderivative_f__(1)) <> [11 13]) then
88             error(msprintf(gettext("%s: Wrong type for argument #%d: Function expected in first element of list.\n"), "numderivative", 1));
89         end
90     end
91     if type(order) ~= 1 then
92         error(msprintf(gettext("%s: Wrong type for argument #%d: Matrix expected.\n"), "numderivative", 4));
93     end
94     if type(H_form) ~= 10 then
95         error(msprintf(gettext("%s: Wrong type for input argument #%d: String array expected.\n"), "numderivative", 5));
96     end
97     if type(Q) ~= 1 then
98         error(msprintf(gettext("%s: Wrong type for argument #%d: Matrix expected.\n"), "numderivative", 6));
99     end
100     //
101     // Check sizes
102     if type(__numderivative_f__) == 15 then
103         // List case
104         if length(__numderivative_f__) < 2 then
105             error(msprintf(gettext("%s: Wrong number of elements in input argument #%d: At least %d elements expected, but current number is %d.\n"), "numderivative", 1, 2, length(__numderivative_f__)));
106         end
107     end
108     if or(size(order) <> [1 1]) then
109         error(msprintf(gettext("%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n"), "numderivative", 4, 1, 1));
110     end
111     if or(size(H_form) <> [1 1]) then
112         error(msprintf(gettext("%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n"), "numderivative", 5, 1, 1));
113     end
114     if or(size(Q) <> [n n]) then
115         error(msprintf(gettext("%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n"), "numderivative", 6, n ,n));
116     end
117     //
118     // Check value
119     if h <> [] then
120         if or(h < 0) then
121             error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be > %d.\n"), "numderivative", 3, 0));
122         end
123     end
124     if and(order <> [1 2 4]) then
125         error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n"), "numderivative", 4, "1, 2, 4"));
126     end
127     if and(H_form <> ["default" "blockmat" "hypermat"]) then
128         error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n"), "numderivative", 5, "default, blockmat, hypermat"));
129     end
130     if norm(clean(Q*Q'-eye(n, n))) > 0 then
131         error(msprintf(gettext("%s: Q must be orthogonal.\n"), "numderivative"));
132     end
133     //
134     // Proceed...
135     if h == [] then
136         h_not_given = %t;
137     else
138         h_not_given = %f;
139         // If h is scalar, expand to the same size as x.
140         if size(h) == [1 1] then
141             h = h * ones(x);
142         end
143     end
144     //
145     // Compute Jacobian
146     if ( h_not_given ) then
147         h = numderivative_step(x, order, 1);
148     end
149     J = numderivative_deriv1(__numderivative_f__, x, h, order, Q);
150     m = size(J, 1);
151     //
152     // Quick return if possible
153     if lhs == 1 then
154         return
155     end
156     //
157     // Compute Hessian matrix
158     if ( h_not_given ) then
159         h = numderivative_step(x, order, 2);
160     end
161     funForHList = list(numderivative_funForH, __numderivative_f__, h, order, Q);
162     if ~Q_not_given then
163         H = numderivative_deriv1(funForHList, x, h, order, Q);
164     else
165         H = numderivative_deriv2(funForHList, x, h, order, Q);
166     end
167     //
168     // At this point, H is a m*n-by-n block matrix.
169     // Update the format of the Hessian
170     if H_form == "default" then
171         // H has the old scilab form
172         H = matrix(H', n*n, m)'
173     end
174     if H_form == "hypermat" then
175         if m > 1 then
176             // H is a hypermatrix if m > 1
177             H = H';
178             H = hypermat([n n m], H(:));
179         end
180     end
181 endfunction
182
183 //
184 // numderivative_step --
185 //   Returns the step for given x, given order and given derivative:
186 //   d = 1 is for Jacobian
187 //   d = 2 is for Hessian
188 // Uses the optimal step.
189 // Then scale the step depending on abs(x).
190 function h = numderivative_step(x, order, d)
191     n = size(x, "*");
192     select d
193     case 1
194         // For Jacobian
195         select order
196         case 1
197             hdefault = sqrt(%eps);
198         case 2
199             hdefault = %eps^(1/3);
200         case 4
201             hdefault = %eps^(1/5);
202         else
203             lclmsg = gettext("%s: Unknown value %s for option %s.\n");
204             error(msprintf(lclmsg,"numderivative_step", string(d), "d"));
205         end
206     case 2
207         // For Hessian
208         select order
209         case 1
210             hdefault = %eps^(1/3);
211         case 2
212             hdefault = %eps^(1/4);
213         case 4
214             hdefault = %eps^(1/6);
215         else
216             lclmsg = gettext("%s: Unknown value %s for option %s.\n");
217             error(msprintf(lclmsg, "numderivative_step", string(d), "d"));
218         end
219     else
220         lclmsg = gettext("%s: Unknown value %s for option %s.\n");
221         error(msprintf(lclmsg, "numderivative_step", string(order), "order"));
222     end
223     // Convert this scalar into a vector, with same size as x
224     // For zero entries in x, use the default.
225     // For nonzero entries in x, scales by abs(x).
226     h = hdefault * abs(x);
227     h(x==0) = hdefault;
228 endfunction
229
230 //
231 // numderivative_funForH --
232 //   Returns the numerical derivative of __numderivative_f__.
233 //   This function is called to compute the numerical Hessian.
234 function J = numderivative_funForH(x, __numderivative_f__, h, order, Q)
235     // Transpose !
236     J = numderivative_deriv1(__numderivative_f__, x, h, order, Q)';
237     J = J(:);
238 endfunction
239
240 // numderivative_deriv1 --
241 //   Computes the numerical gradient of __numderivative_f__, using the given step h.
242 //   This function is used for the computation of the jacobian matrix.
243 function g = numderivative_deriv1(__numderivative_f__, x, h, order, Q)
244     n = size(x, "*");
245     Dy = []; // At this point, we do not know 'm' yet, so we cannot allocate Dy.
246     select order
247     case 1
248         D = Q * diag(h);
249         y = numderivative_evalf(__numderivative_f__, x);
250         for i=1:n
251             d = D(:, i);
252             yplus = numderivative_evalf(__numderivative_f__, x+d);
253             Dyi = (yplus-y)/h(i);
254             Dy = [Dy Dyi];
255         end
256         g = Dy*Q';
257     case 2
258         D = Q * diag(h);
259         for i=1:n
260             d = D(:, i);
261             yplus = numderivative_evalf(__numderivative_f__, x+d);
262             yminus = numderivative_evalf(__numderivative_f__, x-d);
263             Dyi = (yplus-yminus)/(2*h(i));
264             Dy = [Dy Dyi];
265         end
266         g = Dy*Q';
267     case 4
268         D = Q * diag(h);
269         for i=1:n
270             d = D(:, i);
271             yplus = numderivative_evalf(__numderivative_f__, x+d);
272             yminus = numderivative_evalf(__numderivative_f__, x-d);
273             yplus2 = numderivative_evalf(__numderivative_f__, x+2*d);
274             yminus2 = numderivative_evalf(__numderivative_f__, x-2*d);
275             dFh =  (yplus-yminus)/(2*h(i));
276             dF2h = (yplus2-yminus2)/(4*h(i));
277             Dyi = (4*dFh - dF2h)/3;
278             Dy = [Dy Dyi];
279         end
280         g = Dy*Q';
281     end
282 endfunction
283
284 // numderivative_deriv2 --
285 //   Computes the numerical gradient of the argument __numderivative_f__, using the given step h.
286 //   This function is used for the computation of the hessian matrix, to take advantage of its symmetry
287 function g = numderivative_deriv2(__numderivative_f__, x, h, order, Q)
288     n = size(x, "*");
289     Dy = zeros(m*n, n); // 'm' is known at this point, so we can allocate Dy to reduce memory operations
290     select order
291     case 1
292         D = Q * diag(h);
293         y = numderivative_evalf(__numderivative_f__, x);
294         for i=1:n
295             d = D(:, i);
296             yplus = numderivative_evalf(__numderivative_f__, x+d);
297             for j=0:m-1
298                 Dyi(1+j*n:i-1+j*n) = Dy(i+j*n, 1:i-1)'; // Retrieving symmetric elements (will not be done for the first vector)
299                 Dyi(i+j*n:(j+1)*n) = (yplus(i+j*n:(j+1)*n)-y(i+j*n:(j+1)*n))/h(i); // Computing the new ones
300             end
301             Dy(:, i) = Dyi;
302         end
303         g = Dy*Q';
304     case 2
305         D = Q * diag(h);
306         for i=1:n
307             d = D(:, i);
308             yplus = numderivative_evalf(__numderivative_f__, x+d);
309             yminus = numderivative_evalf(__numderivative_f__, x-d);
310             for j=0:m-1
311                 Dyi(1+j*n:i-1+j*n) = Dy(i+j*n, 1:i-1)'; // Retrieving symmetric elements (will not be done for the first vector)
312                 Dyi(i+j*n:(j+1)*n) = (yplus(i+j*n:(j+1)*n)-yminus(i+j*n:(j+1)*n))/(2*h(i)); // Computing the new ones
313             end
314             Dy(:, i) = Dyi;
315         end
316         g = Dy*Q';
317     case 4
318         D = Q * diag(h);
319         for i=1:n
320             d = D(:, i);
321             yplus = numderivative_evalf(__numderivative_f__, x+d);
322             yminus = numderivative_evalf(__numderivative_f__, x-d);
323             yplus2 = numderivative_evalf(__numderivative_f__, x+2*d);
324             yminus2 = numderivative_evalf(__numderivative_f__, x-2*d);
325             for j=0:m-1
326                 dFh(1+j*n:i-1+j*n) = Dy(i+j*n, 1:i-1)'; // Retrieving symmetric elements (will not be done for the first vector)
327                 dFh(i+j*n:(j+1)*n) = (yplus(i+j*n:(j+1)*n)-yminus(i+j*n:(j+1)*n))/(2*h(i)); // Computing the new ones
328                 dF2h(1+j*n:i-1+j*n) = Dy(i+j*n, 1:i-1)'; // Retrieving symmetric elements (will not be done for the first vector)
329                 dF2h(i+j*n:(j+1)*n) = (yplus2(i+j*n:(j+1)*n)-yminus2(i+j*n:(j+1)*n))/(4*h(i)); // Computing the new ones
330             end
331             Dyi = (4*dFh - dF2h)/3;
332             Dy(:, i) = Dyi;
333         end
334         g = Dy*Q';
335     end
336 endfunction
337
338 // numderivative_evalf --
339 // Computes the value of __numderivative_f__ at the point x.
340 // The argument __numderivative_f__ can be a function (macro or linked code) or a list.
341 function y = numderivative_evalf(__numderivative_f__, x)
342     if type(__numderivative_f__) == 15 then
343         // List case
344         __numderivative_fun__ = __numderivative_f__(1);
345         instr = "y = __numderivative_fun__(x, __numderivative_f__(2:$))";
346     elseif or(type(__numderivative_f__) == [11 13 130]) then
347         // Function case
348         instr = "y = __numderivative_f__(x)";
349     else
350         error(msprintf(gettext("%s: Wrong type for input argument #%d: A function expected.\n"), "numderivative", 1));
351     end
352     ierr = execstr(instr, "errcatch")
353     if ierr <> 0 then
354         lamsg = lasterror();
355         lclmsg = "%s: Error while evaluating the function: ""%s""\n";
356         error(msprintf(gettext(lclmsg), "numderivative", lamsg));
357     end
358 endfunction