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