Deleted vectorized computation feature. Deleted neldermead_contour. Fixed the demos.
[scilab.git] / scilab / modules / optimization / help / en_US / optim.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - INRIA
5  * 
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at    
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  -->
13 <refentry version="5.0-subset Scilab" xml:id="optim" xml:lang="en"
14           xmlns="http://docbook.org/ns/docbook"
15           xmlns:xlink="http://www.w3.org/1999/xlink"
16           xmlns:svg="http://www.w3.org/2000/svg"
17           xmlns:ns3="http://www.w3.org/1999/xhtml"
18           xmlns:mml="http://www.w3.org/1998/Math/MathML"
19           xmlns:db="http://docbook.org/ns/docbook">
20   <info>
21     <pubdate>$LastChangedDate: 2008-04-28 09:36:26 +0200 (lun, 28 avr 2008)
22     $</pubdate>
23   </info>
24
25   <refnamediv>
26     <refname>optim</refname>
27
28     <refpurpose>non-linear optimization routine</refpurpose>
29   </refnamediv>
30
31   <refsynopsisdiv>
32     <title>Calling Sequence</title>
33
34     <synopsis>[f,xopt]=optim(costf,x0)
35 [f [,xopt [,gradopt [,work]]]]=optim(costf [,&lt;contr&gt;],x0 [,algo] [,df0 [,mem]] [,work] [,&lt;stop&gt;] [,&lt;params&gt;] [,imp=iflag])</synopsis>
36   </refsynopsisdiv>
37
38   <refsection>
39     <title>Parameters</title>
40
41     <variablelist>
42       <varlistentry>
43         <term>costf</term>
44
45         <listitem>
46           <para>external, i.e Scilab function list or string
47           (<literal>costf</literal> is the cost function, that is, a Scilab
48           script, a Fortran 77 routine or a C function.</para>
49         </listitem>
50       </varlistentry>
51
52       <varlistentry>
53         <term>x0</term>
54
55         <listitem>
56           <para>real vector (initial value of variable to be
57           minimized).</para>
58         </listitem>
59       </varlistentry>
60
61       <varlistentry>
62         <term>f</term>
63
64         <listitem>
65           <para>value of optimal cost
66           (<literal>f=costf(xopt)</literal>)</para>
67         </listitem>
68       </varlistentry>
69
70       <varlistentry>
71         <term>xopt</term>
72
73         <listitem>
74           <para>best value of <literal>x</literal> found.</para>
75         </listitem>
76       </varlistentry>
77
78       <varlistentry>
79         <term>&lt;contr&gt;</term>
80
81         <listitem>
82           <para>keyword representing the following sequence of arguments:
83           <literal>'b',binf,bsup</literal> with <literal>binf</literal> and
84           <literal>bsup</literal> are real vectors with same dimension as
85           <literal>x0</literal>. <literal>binf</literal> and
86           <literal>bsup</literal> are lower and upper bounds on
87           <literal>x</literal>.</para>
88         </listitem>
89       </varlistentry>
90
91       <varlistentry>
92         <term>algo</term>
93
94         <listitem>
95           <itemizedlist>
96             <listitem>
97               <para><literal>'qn'</literal> : quasi-Newton (this is the
98               default solver)</para>
99             </listitem>
100
101             <listitem>
102               <para><literal>'gc'</literal> : conjugate gradient</para>
103             </listitem>
104
105             <listitem>
106               <para><literal>'nd'</literal> : non-differentiable.</para>
107
108               <para>Note that the conjugate gradient solver does not accept
109               bounds on <literal>x</literal>.</para>
110             </listitem>
111           </itemizedlist>
112         </listitem>
113       </varlistentry>
114
115       <varlistentry>
116         <term>df0</term>
117
118         <listitem>
119           <para>real scalar. Guessed decreasing of <literal>f</literal> at
120           first iteration. (<literal>df0=1</literal> is the default
121           value).</para>
122         </listitem>
123       </varlistentry>
124
125       <varlistentry>
126         <term>mem :</term>
127
128         <listitem>
129           <para>integer, number of variables used to approximate the Hessian.
130           Default value is 10. This feature is available for the
131           Gradient-Conjugate algorithm "gc" without constraints and the
132           non-smooth algorithm "nd" without constraints.</para>
133         </listitem>
134       </varlistentry>
135
136       <varlistentry>
137         <term>&lt;stop&gt;</term>
138
139         <listitem>
140           <para>keyword representing the sequence of optional parameters
141           controlling the convergence of the algorithm. <literal>'ar',nap
142           [,iter [,epsg [,epsf [,epsx]]]]</literal></para>
143
144           <variablelist>
145             <varlistentry>
146               <term>"ar"</term>
147
148               <listitem>
149                 <para>reserved keyword for stopping rule selection defined as
150                 follows:</para>
151               </listitem>
152             </varlistentry>
153
154             <varlistentry>
155               <term>nap</term>
156
157               <listitem>
158                 <para>maximum number of calls to <literal>costf</literal>
159                 allowed (default is 100).</para>
160               </listitem>
161             </varlistentry>
162
163             <varlistentry>
164               <term>iter</term>
165
166               <listitem>
167                 <para>maximum number of iterations allowed (default is
168                 100).</para>
169               </listitem>
170             </varlistentry>
171
172             <varlistentry>
173               <term>epsg</term>
174
175               <listitem>
176                 <para>threshold on gradient norm.</para>
177               </listitem>
178             </varlistentry>
179
180             <varlistentry>
181               <term>epsf</term>
182
183               <listitem>
184                 <para>threshold controlling decreasing of
185                 <literal>f</literal></para>
186               </listitem>
187             </varlistentry>
188
189             <varlistentry>
190               <term>epsx</term>
191
192               <listitem>
193                 <para>threshold controlling variation of <literal>x</literal>.
194                 This vector (possibly matrix) of same size as
195                 <literal>x0</literal> can be used to scale
196                 <literal>x</literal>.</para>
197               </listitem>
198             </varlistentry>
199           </variablelist>
200         </listitem>
201       </varlistentry>
202
203       <varlistentry>
204         <term>&lt;params&gt;</term>
205
206         <listitem>
207           <para>keyword representing the method to initialize the arguments
208           <literal>ti, td</literal> passed to the objective function, provided
209           as a C or Fortran routine. This option has no meaning when the cost
210           function is a Scilab script. &lt;params&gt; can be set to only one
211           of the following values.</para>
212
213           <itemizedlist>
214             <listitem>
215               <para>"in"</para>
216
217               <para>That mode allows to allocate memory in the internal Scilab
218               workspace so that the objective function can get arrays with the
219               required size, but without directly allocating the memory. "in"
220               stands for "initialization". In that mode, before the value and
221               derivative of the objective function is to be computed, there is
222               a dialog between the optim Scilab primitive and the objective
223               function. In this dialog, the objective function is called two
224               times, with particular values of the "ind" parameter. The first
225               time, ind is set to 10 and the objective function is expected to
226               set the nizs, nrzs and ndzs integer parameters of the "nird"
227               common.</para>
228
229               <programlisting role = ""><![CDATA[ 
230 common /nird/ nizs,nrzs,ndzs
231  ]]></programlisting>
232
233               <para>This allows Scilab to allocate memory inside its internal
234               workspace. The second time the objective function is called, ind
235               is set to 11 and the objective function is expected to set the
236               ti, tr and tz arrays. After this initialization phase, each time
237               it is called, the objective function is ensured that the ti, tr
238               and tz arrays which are passed to it have the values that have
239               been previously initialized.</para>
240             </listitem>
241
242             <listitem>
243               <para>"ti",valti</para>
244
245               <para>In this mode, valti is expected to be a Scilab vector
246               variable containing integers. Whenever the objective function is
247               called, the ti array it receives contains the values of the
248               Scilab variable.</para>
249             </listitem>
250
251             <listitem>
252               <para>"td", valtd</para>
253
254               <para>In this mode, valtd is expected to be a Scilab vector
255               variable containing double values. Whenever the objective
256               function is called, the td array it receives contains the values
257               of the Scilab variable.</para>
258             </listitem>
259
260             <listitem>
261               <para>"ti",valti,"td",valtd</para>
262
263               <para>This mode combines the two previous.</para>
264             </listitem>
265           </itemizedlist>
266
267           <para>The <literal>ti, td</literal> arrays may be used so that the
268           objective function can be computed. For example, if the objective
269           function is a polynomial, the ti array may may be used to store the
270           coefficients of that polynomial.</para>
271
272           <para>Users should choose carefully between the "in" mode and the
273           "ti" and "td" mode, depending on the fact that the arrays are Scilab
274           variables or not. If the data is available as Scilab variables, then
275           the "ti", valti, "td", valtd mode should be chosen. If the data is
276           available directly from the objective function, the "in" mode should
277           be chosen. Notice that there is no "tr" mode, since, in Scilab, all
278           real values are of "double" type.</para>
279
280           <para>If neither the "in" mode, nor the "ti", "td" mode is chosen,
281           that is, if &lt;params&gt; is not present as an option of the optim
282           primitive, the user may should not assume that the ti,tr and td
283           arrays can be used : reading or writing the arrays may generate
284           unpredictable results.</para>
285         </listitem>
286       </varlistentry>
287
288       <varlistentry>
289         <term>"imp=iflag"</term>
290
291         <listitem>
292           <para>named argument used to set the trace mode. The possible values
293           for iflag are 0,1,2 and &gt;2. Use this option with caution : most
294           of these reports are written on the Scilab standard output.</para>
295
296           <itemizedlist>
297             <listitem>
298               <para>iflag=0: nothing (except errors) is reported (this is the
299               default),</para>
300             </listitem>
301
302             <listitem>
303               <para>iflag=1: initial and final reports,</para>
304             </listitem>
305
306             <listitem>
307               <para>iflag=2: adds a report per iteration,</para>
308             </listitem>
309
310             <listitem>
311               <para>iflag&gt;2: add reports on linear search.</para>
312             </listitem>
313
314             <listitem>
315               <para>iflag&lt;0: calls the cost function with ind=1 every -imp iterations.</para>
316             </listitem>
317           </itemizedlist>
318         </listitem>
319       </varlistentry>
320
321       <varlistentry>
322         <term>gradopt</term>
323
324         <listitem>
325           <para>gradient of <literal>costf</literal> at
326           <literal>xopt</literal></para>
327         </listitem>
328       </varlistentry>
329
330       <varlistentry>
331         <term>work</term>
332
333         <listitem>
334           <para>working array for hot restart for quasi-Newton method. This
335           array is automatically initialized by <literal>optim</literal> when
336           <literal>optim</literal> is invoked. It can be used as input
337           parameter to speed-up the calculations.</para>
338         </listitem>
339       </varlistentry>
340     </variablelist>
341   </refsection>
342
343   <refsection>
344     <title>Description</title>
345
346     <para>Non-linear optimization routine for programs without constraints or
347     with bound constraints:</para>
348
349     <programlisting role = ""><![CDATA[ 
350 min costf(x) w.r.t x.
351  ]]></programlisting>
352
353     <para><literal>costf</literal> is an "external" i.e a Scilab function, a
354     list or a string giving the name of a C or Fortran routine (see
355     "external"). This external must return the value <literal>f</literal> of
356     the cost function at the point <literal>x</literal> and the gradient
357     <literal>g</literal> of the cost function at the point
358     <literal>x</literal>.</para>
359
360     <variablelist>
361       <varlistentry>
362         <term>- Scilab function case</term>
363
364         <listitem>
365           <para>If <literal>costf</literal> is a Scilab function, the calling
366           sequence for <literal>costf</literal> must be:</para>
367
368           <programlisting role = ""><![CDATA[ 
369 [f,g,ind]=costf(x,ind)
370  ]]></programlisting>
371
372           <para>Here, <literal>costf</literal> is a function which returns
373           <literal>f</literal>, value (real number) of cost function at
374           <literal>x</literal>, and <literal>g</literal>, gradient vector of
375           cost function at <literal>x</literal>. The variable
376           <literal>ind</literal> is described below.</para>
377         </listitem>
378       </varlistentry>
379
380       <varlistentry>
381         <term>- List case</term>
382
383         <listitem>
384           <para>If <literal>costf</literal> is a list, it should be of the
385           form: <literal>list(real_costf, arg1,...,argn)</literal> with
386           <literal>real_costf</literal> a Scilab function with calling
387           sequence : <literal>[f,g,ind]=costf(x,ind,arg1,... argn)</literal>.
388           The <literal>x</literal>, <literal>f</literal>,
389           <literal>g</literal>, <literal>ind</literal> arguments have the same
390           meaning that above. <literal>argi</literal> arguments can be used to
391           pass function parameters.</para>
392         </listitem>
393       </varlistentry>
394
395       <varlistentry>
396         <term>- String case</term>
397
398         <listitem>
399           <para>If <literal>costf</literal> is a character string, it refers
400           to the name of a C or Fortran routine which must be linked to
401           Scilab</para>
402
403           <variablelist>
404             <varlistentry>
405               <term>* Fortran case</term>
406
407               <listitem>
408                 <para>The interface of the Fortran subroutine computing the
409                 objective must be :</para>
410
411                 <programlisting role = ""><![CDATA[ 
412 subroutine costf(ind,n,x,f,g,ti,tr,td)
413  ]]></programlisting>
414
415                 <para>with the following declarations:</para>
416
417                 <programlisting role = ""><![CDATA[ 
418 integer ind,n ti(*)
419 double precision x(n),f,g(n),td(*)
420 real tr(*)
421  ]]></programlisting>
422
423                 <para>The argument <literal>ind</literal> is described
424                 below.</para>
425
426                 <para>If ind = 2, 3 or 4, the inputs of the routine are :
427                 <literal>x, ind, n, ti, tr,td</literal>.</para>
428
429                 <para>If ind = 2, 3 or 4, the outputs of the routine are :
430                 <literal>f</literal> and <literal>g</literal>.</para>
431               </listitem>
432             </varlistentry>
433
434             <varlistentry>
435               <term>* C case</term>
436
437               <listitem>
438                 <para>The interface of the C function computing the objective
439                 must be :</para>
440
441                 <programlisting role = ""><![CDATA[ 
442 void costf(int *ind, int *n, double *x, double *f, double *g, int *ti, float *tr, double *td)
443  ]]></programlisting>
444
445                 <para>The argument <literal>ind</literal> is described
446                 below.</para>
447
448                 <para>The inputs and outputs of the function are the same as
449                 in the fortran case.</para>
450               </listitem>
451             </varlistentry>
452           </variablelist>
453         </listitem>
454       </varlistentry>
455     </variablelist>
456
457     <para>If <literal>ind=2</literal> (resp. <literal>3, 4</literal>),
458     <literal>costf</literal> must provide <literal>f</literal> (resp.
459     <literal>g, f</literal> and <literal>g</literal>).</para>
460
461     <para>If <literal>ind=1</literal> nothing is computed (used for display
462     purposes only).</para>
463
464     <para>On output, <literal>ind&lt;0</literal> means that
465     <literal>f</literal> cannot be evaluated at <literal>x</literal> and
466     <literal>ind=0</literal> interrupts the optimization.</para>
467   </refsection>
468
469   <refsection>
470     <title>Example #1 : Scilab function</title>
471
472     <para>The following is an example with a Scilab function. Notice, for
473     simplifications reasons, the Scilab function "cost" of the following
474     example computes the objective function f and its derivative no matter of
475     the value of ind. This allows to keep the example simple. In practical
476     situations though, the computation of "f" and "g" may raise performances
477     issues so that a direct optimization may be to use the value of "ind" to
478     compute "f" and "g" only when needed.</para>
479
480     <programlisting role="example"><![CDATA[ 
481 // External function written in Scilab
482 xref=[1;2;3];
483 x0=[1;-1;1];
484 function [f,g,ind] = cost(x,ind)
485   f=0.5*norm(x-xref)^2;
486   g=x-xref;
487 endfunction
488
489 // Simplest call
490 [f,xopt]=optim(cost,x0)
491
492 // By conjugate gradient - you can use 'qn', 'gc' or 'nd'
493 [f,xopt,gopt]=optim(cost,x0,'gc')
494
495 //Seen as non differentiable
496 [f,xopt,gopt]=optim(cost,x0,'nd')
497
498 // Upper and lower bounds on x
499 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0)
500
501 // Upper and lower bounds on x and setting up the algorithm to 'gc'
502 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc')
503
504 // Bound on the number of call to the objective function
505 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc','ar',3)
506
507 // Set max number of call to the objective function (3)
508 // Set max number of iterations (100)
509 // Set stopping threshold on the value of f (1e-6), 
510 // on the value of the norm of the gradient of the objective function (1e-6)
511 // on the improvement on the parameters x_opt (1e-6;1e-6;1e-6)
512 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc','ar',3,100,1e-6,1e-6,[1e-3;1e-3;1e-3])
513
514 // Print information messages while optimizing
515 // Be careful, some messages are printed in a terminal. You must
516 // Scilab from the command line to see these messages.
517 [f,xopt]=optim(cost,x0,imp=3)
518
519 // Use the 'derivative' function to compute the partial
520 // derivatives of the previous problem
521 deff('y=my_f(x)','y=0.5*norm(x-xref)^2');
522 deff('y=my_df(x)','y=derivative(my_f,x)');
523 deff('[f,g,ind]=cost(x,ind)','f=my_f(x); ...
524                               g=my_df(x)');
525
526 // Simplest call
527 xref=[1;2;3];x0=[1;-1;1]
528 [f,xopt]=optim(cost,x0)
529  ]]></programlisting>
530   </refsection>
531
532   <refsection>
533     <title>Example #2 : C function</title>
534
535     <para>The following is an example with a C function, where a C source code
536     is written into a file, dynamically compiled and loaded into Scilab, and
537     then used by the "optim" solver. The interface of the "rosenc" function is
538     fixed, even if the arguments are not really used in the cost function.
539     This is because the underlying optimization solvers must assume that the
540     objective function has a known, constant interface. In the following
541     example, the arrays ti and tr are not used, only the array "td" is used,
542     as a parameter of the Rosenbrock function. Notice that the content of the
543     arrays ti and td are the same that the content of the Scilab variable, as
544     expected.</para>
545
546     <programlisting role="example"><![CDATA[ 
547 // External function written in C (C compiler required)
548 // write down the C code (Rosenbrock problem)
549 C=['#include <math.h>'
550    'double sq(double x)'
551    '{ return x*x;}'
552    'void rosenc(int *ind, int *n, double *x, double *f, double *g, '
553    '                                int *ti, float *tr, double *td)'
554    '{'
555    '  double p;'
556    '  int i;'
557    '  p=td[0];'
558    '  if (*ind==2||*ind==4) {'
559    '    *f=1.0;'
560    '    for (i=1;i<*n;i++)'
561    '      *f+=p*sq(x[i]-sq(x[i-1]))+sq(1.0-x[i]);'
562    '  }'
563    '  if (*ind==3||*ind==4) {'
564    '    g[0]=-4.0*p*(x[1]-sq(x[0]))*x[0];'
565    '    for (i=1;i<*n-1;i++)'
566    '      g[i]=2.0*p*(x[i]-sq(x[i-1]))-4.0*p*(x[i+1]-sq(x[i]))*x[i]-2.0*(1.0-x[i]);'
567    '    g[*n-1]=2.0*p*(x[*n-1]-sq(x[*n-2]))-2.0*(1.0-x[*n-1]);'
568    '  }'
569    '}'];
570 mputl(C,TMPDIR+'/rosenc.c')
571
572 // compile the C code
573 l=ilib_for_link('rosenc','rosenc.o',[],'c',TMPDIR+'/Makefile');
574
575 // incremental linking
576 link(l,'rosenc','c')
577
578 //solve the problem
579 x0=[40;10;50];
580 p=100;
581 [f,xo,go]=optim('rosenc',x0,'td',p)
582  ]]></programlisting>
583   </refsection>
584
585   <refsection>
586     <title>Example #3 : Fortran function</title>
587
588     <para>The following is an example with a Fortran function.</para>
589
590     <programlisting role="example"><![CDATA[ 
591 // External function written in Fortran (Fortran compiler required)
592 // write down the Fortran  code (Rosenbrock problem)
593 F=[ '      subroutine rosenf(ind, n, x, f, g, ti, tr, td)'
594     '      integer ind,n,ti(*)'
595     '      double precision x(n),f,g(n),td(*)'
596     '      real tr(*)'
597     'c'
598     '      double precision y,p'
599     '      p=td(1)'
600     '      if (ind.eq.2.or.ind.eq.4) then'
601     '        f=1.0d0'
602     '        do i=2,n'
603     '          f=f+p*(x(i)-x(i-1)**2)**2+(1.0d0-x(i))**2'
604     '        enddo'
605     '      endif'
606     '      if (ind.eq.3.or.ind.eq.4) then'
607     '        g(1)=-4.0d0*p*(x(2)-x(1)**2)*x(1)'
608     '        if(n.gt.2) then'
609     '          do i=2,n-1'
610     '            g(i)=2.0d0*p*(x(i)-x(i-1)**2)-4.0d0*p*(x(i+1)-x(i)**2)*x(i)'
611     '     &           -2.0d0*(1.0d0-x(i))'
612     '          enddo'
613     '        endif'
614     '        g(n)=2.0d0*p*(x(n)-x(n-1)**2)-2.0d0*(1.0d0-x(n))'
615     '      endif'
616     '      return'
617     '      end'];
618
619 mputl(F,TMPDIR+'/rosenf.f')
620
621 // compile the Fortran code
622 l=ilib_for_link('rosenf','rosenf.o',[],'f',TMPDIR+'/Makefile');
623
624 // incremental linking
625 link(l,'rosenf','f')
626
627 //solve the problem
628 x0=[40;10;50];
629 p=100;
630 [f,xo,go]=optim('rosenf',x0,'td',p)
631  ]]></programlisting>
632   </refsection>
633
634   <refsection>
635     <title>Example #4 : Fortran function with initialization</title>
636
637     <para>The following is an example with a Fortran function in which the
638     "in" option is used to allocate memory inside the Scilab environment. In
639     this mode, there is a dialog between Scilab and the objective function.
640     The goal of this dialog is to initialize the parameters of the objective
641     function. Each part of this dialog is based on a specific value of the
642     "ind" parameter.</para>
643
644     <para>At the beginning, Scilab calls the objective function, with the ind
645     parameter equals to 10. This tells the objective function to initialize
646     the sizes of the arrays it needs by setting the nizs, nrzs and ndzs
647     integer parameters of the "nird" common. Then the objective function
648     returns. At this point, Scilab creates internal variables and allocate
649     memory for the variable izs, rzs and dzs. Scilab calls the objective
650     function back again, this time with ind equals to 11. This tells the
651     objective function to initialize the arrays izs, rzs and dzs. When the
652     objective function has done so, it returns. Then Scilab enters in the real
653     optimization mode and calls the optimization solver the user requested.
654     Whenever the objective function is called, the izs, rzs and dzs arrays
655     have the values that have been previously initialized.</para>
656
657     <programlisting role="example"><![CDATA[ 
658 // 
659 // Define a fortran source code and compile it (fortran compiler required)
660 //
661 fortransource=['      subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
662                'C     -------------------------------------------'
663                'c     Example of cost function given by a subroutine'
664                'c     if n<=2 returns ind=0'
665                'c     f.bonnans, oct 86'
666                '      implicit double precision (a-h,o-z)'
667                '      real rzs(1)'
668                '      double precision dzs(*)'
669                '      dimension x(n),g(n),izs(*)'
670                '      common/nird/nizs,nrzs,ndzs'
671                '      if (n.lt.3) then'
672                '        ind=0'
673                '        return'
674                '      endif'
675                '      if(ind.eq.10) then'
676                '         nizs=2'
677                '         nrzs=1'
678                '         ndzs=2'
679                '         return'
680                '      endif'
681                '      if(ind.eq.11) then'
682                '         izs(1)=5'
683                '         izs(2)=10'
684                '         dzs(2)=100.0d+0'
685                '         return'
686                '      endif'
687                '      if(ind.eq.2)go to 5'
688                '      if(ind.eq.3)go to 20'
689                '      if(ind.eq.4)go to 5'
690                '      ind=-1'
691                '      return'
692                '5     f=1.0d+0'
693                '      do 10 i=2,n'
694                '        im1=i-1'
695                '10      f=f + dzs(2)*(x(i)-x(im1)**2)**2 + (1.0d+0-x(i))**2'
696                '      if(ind.eq.2)return'
697                '20    g(1)=-4.0d+0*dzs(2)*(x(2)-x(1)**2)*x(1)'
698                '      nm1=n-1'
699                '      do 30 i=2,nm1'
700                '        im1=i-1'
701                '        ip1=i+1'
702                '        g(i)=2.0d+0*dzs(2)*(x(i)-x(im1)**2)'
703                '30      g(i)=g(i) -4.0d+0*dzs(2)*(x(ip1)-x(i)**2)*x(i) - '
704                '     &        2.0d+0*(1.0d+0-x(i))'
705                '      g(n)=2.0d+0*dzs(2)*(x(n)-x(nm1)**2) - 2.0d+0*(1.0d+0-x(n))'
706                '      return'
707                '      end'];
708 mputl(fortransource,TMPDIR+'/rosenf.f')
709
710 // compile the C code
711 libpath=ilib_for_link('rosenf','rosenf.o',[],'f',TMPDIR+'/Makefile');
712
713 // incremental linking
714 linkid=link(libpath,'rosenf','f');
715
716 x0=1.2*ones(1,5);
717 //
718 // Solve the problem
719 //
720 [f,x,g]=optim('rosenf',x0,'in');
721  ]]></programlisting>
722   </refsection>
723
724   <refsection>
725     <title>Example #5 : Fortran function with initialization on Windows with
726     Intel Fortran Compiler</title>
727
728     <para>Under the Windows operating system with Intel Fortran Compiler, one
729     must carefully design the fortran source code so that the dynamic link
730     works properly. On Scilab's side, the optimization component is
731     dynamically linked and the symbol "nird" is exported out of the
732     optimization dll. On the cost function's side, which is also dynamically
733     linked, the "nird" common must be imported in the cost function
734     dll.</para>
735
736     <para>The following example is a re-writing of the previous example, with
737     special attention for the Windows operating system with Intel Fortran
738     compiler as example. In that case, we introduce additionnal compiling
739     instructions, which allows the compiler to import the "nird"
740     symbol.</para>
741
742     <programlisting role="example"><![CDATA[ 
743 fortransource=['subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
744                'cDEC$ IF DEFINED (FORDLL)'
745                'cDEC$ ATTRIBUTES DLLIMPORT:: /nird/'
746                'cDEC$ ENDIF'
747                'C     -------------------------------------------'
748                'c     Example of cost function given by a subroutine'
749                'c     if n<=2 returns ind=0'
750                'c     f.bonnans, oct 86'
751                '      implicit double precision (a-h,o-z)'
752 [etc...]
753  ]]></programlisting>
754   </refsection>
755
756   <refsection>
757     <title>Example #6 : Logging features</title>
758
759     <para>The imp flag may take negative integer values, say k.
760     In that case, the cost function is called once every -k iterations.
761     This allows to draw the function value or write a log file.
762     </para>
763
764     <para>In the following example, we solve the Rosenbrock test case.
765     For each iteration of the algorithm, we print the value of x, f and g.
766     </para>
767
768 <programlisting role="example">
769 xref=[1;2;3];
770 x0=[1;-1;1];
771 function [f,g,ind] = cost(x,ind)
772   f=0.5*norm(x-xref)^2;
773   g=x-xref;
774 endfunction
775 function [f,g,ind] = cost(x,ind)
776   if ind == 2 | ind == 4 then
777     f=0.5*norm(x-xref)^2;
778   end
779   if ind == 3 | ind == 4 then
780     g=x-xref;
781   end
782   if ind == 1 then
783     mprintf("x = %s\n", strcat(string(x)," "))
784     mprintf("f = %e\n", f)
785     g=x-xref;
786     mprintf("g = %s\n", strcat(string(g)," "))
787   end
788 endfunction
789 [f,xopt]=optim(cost,x0,imp=-1)
790   </programlisting>
791
792     <para>In the following example, we solve the Rosenbrock test case.
793     For each iteration of the algorithm, we plot the current value of x
794     into a 2D graph containing the contours of Rosenbrock's function.
795     This allows to see the progress of the algorithm while the 
796     algorithm is performing. We could as well write the 
797     value of x, f and g into a log file if needed.
798     </para>
799
800 <programlisting role="example">
801 // 1. Define rosenbrock
802 function [ f , g , ind ] = rosenbrock ( x , ind )
803   if ((ind == 1) | (ind == 2) | (ind == 4)) then
804     f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
805   end
806   if ((ind == 1) | (ind == 3) | (ind == 4)) then
807     g(1) = - 400. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
808     g(2) = 200. * ( x(2) - x(1)**2 )
809   end
810   if (ind == 1) then
811     plot ( x(1) , x(2) , "g." )
812   end
813 endfunction
814 x0 = [-1.2 1.0];
815 xopt = [1.0 1.0];
816 // 2. Draw the contour of Rosenbrock's function
817 xmin = -2.0
818 xmax = 2.0
819 stepx = 0.1
820 ymin = -1.0
821 ymax = 2.0
822 stepy = 0.1
823 nx = 100
824 ny = 100
825 stepx = (xmax - xmin)/nx;
826 xdata = xmin:stepx:xmax;
827 stepy = (ymax - ymin)/ny;
828 ydata = ymin:stepy:ymax;
829 for ix = 1:length(xdata)
830     for iy = 1:length(ydata)
831       x = [xdata(ix) ydata(iy)];
832       f = rosenbrock ( x , 2 );
833       zdata ( ix , iy ) = f;
834     end
835 end
836 contour ( xdata , ydata , zdata , [1 10 100 500 1000])
837 plot(x0(1) , x0(2) , "b.")
838 plot(xopt(1) , xopt(2) , "r*")
839 // 3. Plot the optimization process, during optimization
840 [ fopt , xopt ] = optim ( rosenbrock , x0 , imp = -1)
841   </programlisting>
842   </refsection>
843
844
845   <refsection>
846     <title>Example #7 : Optimizing without derivatives</title>
847
848     <para>It is possible to optimize a problem without an explicit 
849     knowledge of the derivative of the cost function.
850     For this purpose, we can use the numdiff or derivative function
851     to compute a numerical derivative of the cost function.
852     </para>
853
854     <para>In the following example, we use the numdiff function to 
855     solve Rosenbrock's problem.
856     </para>
857
858
859 <programlisting role="example">
860 function f = rosenbrock ( x )
861   f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
862 endfunction
863 function [ f , g , ind ] = rosenbrockCost ( x , ind )
864   if ((ind == 1) | (ind == 2) | (ind == 4)) then
865     f = rosenbrock ( x );
866   end
867   if ((ind == 1) | (ind == 3) | (ind == 4)) then
868     g= numdiff ( rosenbrock , x );
869   end
870 endfunction
871 x0 = [-1.2 1.0];
872 [ fopt , xopt ] = optim ( rosenbrockCost , x0 )
873   </programlisting>
874
875     <para>In the following example, we use the derivative function to 
876     solve Rosenbrock's problem. Given that the step computation strategy
877     is not the same in numdiff and derivative, this might lead to improved
878     results.
879     </para>
880
881
882 <programlisting role="example">
883 function f = rosenbrock ( x )
884   f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
885 endfunction
886 function [ f , g , ind ] = rosenbrockCost2 ( x , ind )
887   if ((ind == 1) | (ind == 2) | (ind == 4)) then
888     f = rosenbrock ( x );
889   end
890   if ((ind == 1) | (ind == 3) | (ind == 4)) then
891     g= derivative ( rosenbrock , x.' , order = 4 );
892   end
893 endfunction
894 x0 = [-1.2 1.0];
895 [ fopt , xopt ] = optim ( rosenbrockCost2 , x0 )
896   </programlisting>
897
898   </refsection>
899
900
901   <refsection>
902     <title>See Also</title>
903
904     <simplelist type="inline">
905       <member><link linkend="external">external</link></member>
906
907       <member><link linkend="qpsolve">qpsolve</link></member>
908
909       <member><link linkend="datafit">datafit</link></member>
910
911       <member><link linkend="leastsq">leastsq</link></member>
912
913       <member><link linkend="numdiff">numdiff</link></member>
914
915       <member><link linkend="derivative">derivative</link></member>
916
917       <member><link linkend="NDcost">NDcost</link></member>
918     </simplelist>
919   </refsection>
920
921   <refsection>
922     <title>References</title>
923
924     <para>The following is a map from the various options to the underlying
925     solvers, with some comments about the algorithm, when available.</para>
926
927     <variablelist>
928       <varlistentry>
929         <term>"qn" without constraints</term>
930
931         <listitem>
932           <para>n1qn1 : a quasi-Newton method with a Wolfe-type line
933           search</para>
934         </listitem>
935       </varlistentry>
936
937       <varlistentry>
938         <term>"qn" with bounds constraints</term>
939
940         <listitem>
941           <para>qnbd : a quasi-Newton method with projection</para>
942
943           <para>RR-0242 - A variant of a projected variable metric method for
944           bound constrained optimization problems, Bonnans Frederic, Rapport
945           de recherche de l'INRIA - Rocquencourt, Octobre 1983</para>
946         </listitem>
947       </varlistentry>
948
949       <varlistentry>
950         <term>"gc" without constraints</term>
951
952         <listitem>
953           <para>n1qn3 : a conjugate gradient method with BFGS.</para>
954         </listitem>
955       </varlistentry>
956
957       <varlistentry>
958         <term>"gc" with bounds constraints</term>
959
960         <listitem>
961           <para>gcbd : a BFGS-type method with limited memory and
962           projection</para>
963         </listitem>
964       </varlistentry>
965
966       <varlistentry>
967         <term>"nd" without constraints</term>
968
969         <listitem>
970           <para>n1fc1 : a bundle method</para>
971         </listitem>
972       </varlistentry>
973
974       <varlistentry>
975         <term>"nd" with bounds constraints</term>
976
977         <listitem>
978           <para>not available</para>
979         </listitem>
980       </varlistentry>
981     </variablelist>
982   </refsection>
983   <refsection>
984     <title>Author</title>
985     <para>The Modulopt library : J.Frederic Bonnans, Jean-Charles Gilbert, Claude Lemarechal</para>
986     <para>The interfaces to the Modulopt library : J.Frederic Bonnans</para>
987     <para>This help : Michael Baudin</para>
988   </refsection>
989 </refentry>