Optimization help: Minor typo fix
[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  * Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
6  * Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
7  *
8  * This file must be used under the terms of the CeCILL.
9  * This source file is licensed as described in the file COPYING, which
10  * you should have received as part of this distribution.  The terms
11  * are also available at
12  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
13  *
14  -->
15 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"  xml:id="optim" xml:lang="en">
16     <refnamediv>
17         <refname>optim</refname>
18         <refpurpose>non-linear optimization routine</refpurpose>
19     </refnamediv>
20     <refsynopsisdiv>
21         <title>Calling Sequence</title>
22         <synopsis>
23             fopt = optim(costf, x0)
24             fopt = optim(costf [,&lt;contr&gt;],x0 [,algo] [,df0 [,mem]] [,work] [,&lt;stop&gt;] [,&lt;params&gt;] [,imp=iflag])
25             [fopt, xopt] = optim(...)
26             [fopt, xopt, gopt] = optim(...)
27             [fopt, xopt, gopt, work] = optim(...)
28             [fopt, xopt, gopt, work, iters] = optim(...)
29             [fopt, xopt, gopt, work, iters, evals] = optim(...)
30             [fopt, xopt, gopt, work, iters, evals, err] = optim(...)
31             [fopt, xopt, [gopt, work, iters, evals, err], ti, td] = optim(..., params="si","sd")
32         </synopsis>
33     </refsynopsisdiv>
34     <refsection>
35         <title>Arguments</title>
36         <variablelist>
37             <varlistentry>
38                 <term>costf</term>
39                 <listitem>
40                     <para>a function, a list or a string, the objective function.</para>
41                 </listitem>
42             </varlistentry>
43             <varlistentry>
44                 <term>x0</term>
45                 <listitem>
46                     <para>real vector, the initial guess for
47                         <literal>x</literal>.
48                     </para>
49                 </listitem>
50             </varlistentry>
51             <varlistentry>
52                 <term>&lt;contr&gt;</term>
53                 <listitem>
54                     <para>an optional sequence of arguments containing the lower and
55                         upper bounds on <literal>x</literal>. If bounds are required, this
56                         sequence of arguments must be <literal>"b",binf,bsup</literal> where
57                         <literal>binf</literal> and <literal>bsup</literal> are real vectors
58                         with same dimension as <literal>x0</literal>.
59                     </para>
60                 </listitem>
61             </varlistentry>
62             <varlistentry>
63                 <term>algo</term>
64                 <listitem>
65                     <para>a string, the algorithm to use (default
66                         <literal>algo="qn"</literal>).
67                     </para>
68                     <para>The available algorithms are:</para>
69                     <itemizedlist>
70                         <listitem>
71                             <para>
72                                 <literal>"qn"</literal>: Quasi-Newton with BFGS
73                             </para>
74                         </listitem>
75                         <listitem>
76                             <para>
77                                 <literal>"gc"</literal>: limited memory BFGS
78                             </para>
79                         </listitem>
80                         <listitem>
81                             <para>
82                                 <literal>"nd"</literal>: non-differentiable.
83                             </para>
84                             <para>
85                                 The <literal>"nd"</literal> algorithm does not accept
86                                 bounds on <literal>x</literal>.
87                             </para>
88                             <para>
89                                 The <literal>"qn"</literal> cannot be used if
90                                 <literal>size(x0)&gt;46333</literal>.
91                             </para>
92                         </listitem>
93                     </itemizedlist>
94                 </listitem>
95             </varlistentry>
96             <varlistentry>
97                 <term>df0</term>
98                 <listitem>
99                     <para>
100                         real scalar, a guess of the decreasing of <literal>f</literal>
101                         at first iteration. (default <literal>df0=1</literal>).
102                     </para>
103                 </listitem>
104             </varlistentry>
105             <varlistentry>
106                 <term>mem</term>
107                 <listitem>
108                     <para>integer, the number of variables used to approximate the
109                         Hessian (default <literal>mem=10</literal>). This feature is
110                         available for the <literal>"gc"</literal> algorithm without
111                         constraints and the non-smooth algorithm <literal>"nd"</literal>
112                         without constraints.
113                     </para>
114                 </listitem>
115             </varlistentry>
116             <varlistentry>
117                 <term>&lt;stop&gt;</term>
118                 <listitem>
119                     <para>a sequence of arguments containing the parameters controlling
120                         the convergence of the algorithm. The following sequences are
121                         available: <screen>              "ar",nap
122                             "ar",nap,iter
123                             "ar",nap,iter,epsg
124                             "ar",nap,iter,epsg,epsf
125                             "ar",nap,iter,epsg,epsf,epsx
126                         </screen>
127                     </para>
128                     <para>where:</para>
129                     <variablelist>
130                         <varlistentry>
131                             <term>nap</term>
132                             <listitem>
133                                 <para>
134                                     maximum number of calls to <literal>costf</literal>
135                                     allowed (default <literal>nap=100</literal>).
136                                 </para>
137                             </listitem>
138                         </varlistentry>
139                         <varlistentry>
140                             <term>iter</term>
141                             <listitem>
142                                 <para>maximum number of iterations allowed (default
143                                     <literal>iter=100</literal>).
144                                 </para>
145                             </listitem>
146                         </varlistentry>
147                         <varlistentry>
148                             <term>epsg</term>
149                             <listitem>
150                                 <para>threshold on gradient norm (default
151                                     <literal>epsg= %eps</literal>).
152                                 </para>
153                             </listitem>
154                         </varlistentry>
155                         <varlistentry>
156                             <term>epsf</term>
157                             <listitem>
158                                 <para>
159                                     threshold controlling decreasing of <literal>f</literal>
160                                     (default <literal>epsf=0</literal>).
161                                 </para>
162                             </listitem>
163                         </varlistentry>
164                         <varlistentry>
165                             <term>epsx</term>
166                             <listitem>
167                                 <para>
168                                     threshold controlling variation of <literal>x</literal>
169                                     (default <literal>epsx=0</literal>). This vector (possibly
170                                     matrix) with same size as <literal>x0</literal> can be used to
171                                     scale <literal>x</literal>.
172                                 </para>
173                             </listitem>
174                         </varlistentry>
175                     </variablelist>
176                 </listitem>
177             </varlistentry>
178             <varlistentry>
179                 <term>&lt;params&gt;</term>
180                 <listitem>
181                     <para>in the case where the objective function is a C or Fortran
182                         routine, a sequence of arguments containing the method to
183                         communicate with the objective function. This option has no meaning
184                         when the cost function is a Scilab script.
185                     </para>
186                     <para>The available values for &lt;params&gt; are the
187                         following.
188                     </para>
189                     <itemizedlist>
190                         <listitem>
191                             <para>
192                                 <literal>"in"</literal>
193                             </para>
194                             <para>That mode allows to allocate memory in the internal Scilab
195                                 workspace so that the objective function can get arrays with the
196                                 required size, but without directly allocating the memory. The
197                                 <literal>"in"</literal> value stands for "initialization". In
198                                 that mode, before the value and derivative of the objective
199                                 function is to be computed, there is a dialog between the
200                                 <literal>optim</literal> Scilab primitive and the objective
201                                 function <literal>costf</literal>. In this dialog, the objective
202                                 function is called two times, with particular values of the
203                                 <literal>ind</literal> parameter. The first time,
204                                 <literal>ind</literal> is set to 10 and the objective function
205                                 is expected to set the <literal>nizs</literal>,
206                                 <literal>nrzs</literal> and <literal>ndzs</literal> integer
207                                 parameters of the <literal>nird</literal> common, which is
208                                 defined as:
209                             </para>
210                             <screen>common /nird/ nizs,nrzs,ndzs    </screen>
211                             <para>This allows Scilab to allocate memory inside its internal
212                                 workspace. The second time the objective function is called,
213                                 <literal>ind</literal> is set to 11 and the objective function
214                                 is expected to set the <literal>ti</literal>,
215                                 <literal>tr</literal> and <literal>tz</literal> arrays. After
216                                 this initialization phase, each time it is called, the objective
217                                 function is ensured that the <literal>ti</literal>,
218                                 <literal>tr</literal> and <literal>tz</literal> arrays which are
219                                 passed to it have the values that have been previously
220                                 initialized.
221                             </para>
222                         </listitem>
223                         <listitem>
224                             <para>
225                                 <literal>"ti",valti</literal>
226                             </para>
227                             <para>
228                                 In this mode, <literal>valti</literal> is expected to be a
229                                 Scilab vector variable containing integers. Whenever the
230                                 objective function is called, the <literal>ti</literal> array it
231                                 receives contains the values of the Scilab variable.
232                             </para>
233                         </listitem>
234                         <listitem>
235                             <para>
236                                 <literal>"td", valtd</literal>
237                             </para>
238                             <para>
239                                 In this mode, <literal>valtd</literal> is expected to be a
240                                 Scilab vector variable containing double values. Whenever the
241                                 objective function is called, the <literal>td</literal> array it
242                                 receives contains the values of the Scilab variable.
243                             </para>
244                         </listitem>
245                         <listitem>
246                             <para>
247                                 <literal>"ti",valti,"td",valtd</literal>
248                             </para>
249                             <para>This mode combines the two previous modes.</para>
250                         </listitem>
251                         <listitem>
252                             <para>
253                                 <literal>"si"</literal>
254                             </para>
255                             <para>
256                                 In this mode, <literal>valti</literal> is saved in the output
257                                 variable <varname>ti</varname>.
258                             </para>
259                         </listitem>
260                         <listitem>
261                             <para>
262                                 <literal>"sd"</literal>
263                             </para>
264                             <para>
265                                 In this mode, <literal>valtd</literal> is saved in the output
266                                 variable <varname>td</varname>.
267                             </para>
268                         </listitem>
269                         <listitem>
270                             <para>
271                                 <literal>"si","sd"</literal>
272                             </para>
273                             <para>This mode combines the two previous modes.</para>
274                         </listitem>
275                     </itemizedlist>
276                     <para>
277                         The <literal>ti, td</literal> arrays may be used so that the
278                         objective function can be computed. For example, if the objective
279                         function is a polynomial, the ti array may may be used to store the
280                         coefficients of that polynomial.
281                     </para>
282                     <para>Users should choose carefully between the
283                         <literal>"in"</literal> mode and the <literal>"ti"</literal> and
284                         <literal>"td"</literal> mode, depending on the fact that the arrays
285                         are Scilab variables or not. If the data is available as Scilab
286                         variables, then the <literal>"ti", valti, "td", valtd</literal> mode
287                         should be chosen. If the data is available directly from the
288                         objective function, the <literal>"in"</literal> mode should be
289                         chosen. Notice that there is no <literal>"tr"</literal> mode, since,
290                         in Scilab, all real values are doubles.
291                     </para>
292                     <para>If neither the "in" mode, nor the "ti", "td" mode is chosen,
293                         that is, if &lt;params&gt; is not present as an option of the optim
294                         primitive, the user may should not assume that the ti,tr and td
295                         arrays can be used : reading or writing the arrays may generate
296                         unpredictable results.
297                     </para>
298                 </listitem>
299             </varlistentry>
300             <varlistentry>
301                 <term>"imp=iflag"</term>
302                 <listitem>
303                     <para>named argument used to set the trace mode (default
304                         <literal>imp=0</literal>, which prints no messages). If <varname>imp</varname>
305                         is greater or equal to 1, more information are printed, depending on the
306                         algorithm chosen. More precisely:
307                     </para>
308                     <itemizedlist>
309                         <listitem>
310                             <para>
311                                 <literal>"qn"</literal> without constraints: from <literal>iflag=1</literal>
312                                 to <literal>iflag=3</literal>.
313                             </para>
314                             <itemizedlist>
315                                 <listitem>
316                                     <para>
317                                         <literal>iflag>=1</literal>: initial and final print,
318                                     </para>
319                                 </listitem>
320                                 <listitem>
321                                     <para>
322                                         <literal>iflag>=2</literal>: one line per iteration (number of iterations,
323                                         number of calls to f, value of f),
324                                     </para>
325                                 </listitem>
326                                 <listitem>
327                                     <para>
328                                         <literal>iflag>=3</literal>: extra information on line searches.
329                                     </para>
330                                 </listitem>
331                             </itemizedlist>
332                         </listitem>
333                         <listitem>
334                             <para>
335                                 <literal>"qn"</literal> with bounds constraints: from <literal>iflag=1</literal>
336                                 to <literal>iflag=4</literal>.
337                             </para>
338                             <itemizedlist>
339                                 <listitem>
340                                     <para>
341                                         <literal>iflag>=1</literal>: initial and final print,
342                                     </para>
343                                 </listitem>
344                                 <listitem>
345                                     <para>
346                                         <literal>iflag>=2</literal>: one line per iteration (number of iterations,
347                                         number of calls to f, value of f),
348                                     </para>
349                                 </listitem>
350                                 <listitem>
351                                     <para>
352                                         <literal>iflag>=3</literal>: extra information on line searches.
353                                     </para>
354                                 </listitem>
355                             </itemizedlist>
356                         </listitem>
357                         <listitem>
358                             <para>
359                                 <literal>"gc"</literal> without constraints: from <literal>iflag=1</literal>
360                                 to <literal>iflag=5</literal>.
361                             </para>
362                             <itemizedlist>
363                                 <listitem>
364                                     <para>
365                                         <literal>iflag>=1</literal> and <literal>iflag>=2</literal>: initial and final print,
366                                     </para>
367                                 </listitem>
368                                 <listitem>
369                                     <para>
370                                         <literal>iflag=3</literal>: one line per iteration (number of iterations,
371                                         number of calls to f, value of f),
372                                     </para>
373                                 </listitem>
374                                 <listitem>
375                                     <para>
376                                         <literal>iflag>=4</literal>: extra information on lines searches.
377                                     </para>
378                                 </listitem>
379                             </itemizedlist>
380                         </listitem>
381                         <listitem>
382                             <para>
383                                 <literal>"gc"</literal> with bounds constraints: from <literal>iflag=1</literal>
384                                 to <literal>iflag=3</literal>.
385                             </para>
386                             <itemizedlist>
387                                 <listitem>
388                                     <para>
389                                         <literal>iflag>=1</literal>: initial and final print,
390                                     </para>
391                                 </listitem>
392                                 <listitem>
393                                     <para>
394                                         <literal>iflag>=2</literal>: one print per iteration,
395                                     </para>
396                                 </listitem>
397                                 <listitem>
398                                     <para>
399                                         <literal>iflag=3</literal>: extra information.
400                                     </para>
401                                 </listitem>
402                             </itemizedlist>
403                         </listitem>
404                         <listitem>
405                             <para>
406                                 <literal>"nd"</literal> with bounds constraints: from <literal>iflag=1</literal>
407                                 to <literal>iflag=8</literal>.
408                             </para>
409                             <itemizedlist>
410                                 <listitem>
411                                     <para>
412                                         <literal>iflag>=1</literal>: initial and final print,
413                                     </para>
414                                 </listitem>
415                                 <listitem>
416                                     <para>
417                                         <literal>iflag>=2</literal>: one print on each convergence,
418                                     </para>
419                                 </listitem>
420                                 <listitem>
421                                     <para>
422                                         <literal>iflag>=3</literal>: one print per iteration,
423                                     </para>
424                                 </listitem>
425                                 <listitem>
426                                     <para>
427                                         <literal>iflag>=4</literal>: line search,
428                                     </para>
429                                 </listitem>
430                                 <listitem>
431                                     <para>
432                                         <literal>iflag>=5</literal>: various tolerances,
433                                     </para>
434                                 </listitem>
435                                 <listitem>
436                                     <para>
437                                         <literal>iflag>=6</literal>: weight and information on the computation of direction.
438                                     </para>
439                                 </listitem>
440                             </itemizedlist>
441                         </listitem>
442                     </itemizedlist>
443                     <para>
444                         If <varname>imp</varname> is lower than 0, then the cost function is evaluated every
445                         <varname>-m</varname> iterations, with <literal>ind=1</literal>.
446                     </para>
447                 </listitem>
448             </varlistentry>
449             <varlistentry>
450                 <term>fopt</term>
451                 <listitem>
452                     <para>the value of the objective function at the point
453                         <literal>xopt</literal>
454                     </para>
455                 </listitem>
456             </varlistentry>
457             <varlistentry>
458                 <term>xopt</term>
459                 <listitem>
460                     <para>
461                         best value of <literal>x</literal> found.
462                     </para>
463                 </listitem>
464             </varlistentry>
465             <varlistentry>
466                 <term>gopt</term>
467                 <listitem>
468                     <para>the gradient of the objective function at the point
469                         <literal>xopt</literal>
470                     </para>
471                 </listitem>
472             </varlistentry>
473             <varlistentry>
474                 <term>work</term>
475                 <listitem>
476                     <para>working array for hot restart for quasi-Newton method. This
477                         array is automatically initialized by <literal>optim</literal> when
478                         <literal>optim</literal> is invoked. It can be used as input
479                         parameter to speed-up the calculations.
480                     </para>
481                 </listitem>
482             </varlistentry>
483             <varlistentry>
484                 <term>iters</term>
485                 <listitem>
486                     <para>
487                         scalar, the number of iterations that is displayed when <literal>imp=2</literal>.
488                     </para>
489                 </listitem>
490             </varlistentry>
491             <varlistentry>
492                 <term>evals</term>
493                 <listitem>
494                     <para>
495                         scalar, the number of <literal>cost</literal> function evaluations
496                         that is displayed when <literal>imp=2</literal>.
497                     </para>
498                 </listitem>
499             </varlistentry>
500             <varlistentry>
501                 <term>err</term>
502                 <listitem>
503                     <para>
504                         scalar, a termination indicator.
505                         The success flag is <literal>9</literal>.
506                         <literal>err=1</literal>: Norm of projected gradient lower than...
507                         <literal>err=2</literal>: At last iteration f decreases by less than...
508                         <literal>err=3</literal>: Optimization stops because of too small variations for x.
509                         <literal>err=4</literal>: Optim stops: maximum number of calls to f is reached.
510                         <literal>err=5</literal>: Optim stops: maximum number of iterations is reached.
511                         <literal>err=6</literal>: Optim stops: too small variations in gradient direction.
512                         <literal>err=7</literal>: Stop during calculation of descent direction.
513                         <literal>err=8</literal>: Stop during calculation of estimated hessian.
514                         <literal>err=9</literal>: End of optimization, successful completion.
515                         <literal>err=10</literal>: End of optimization (linear search fails).
516                     </para>
517                 </listitem>
518             </varlistentry>
519         </variablelist>
520     </refsection>
521     <refsection>
522         <title>Description</title>
523         <para>This function solves unconstrained nonlinear optimization
524             problems:
525         </para>
526         <screen>min f(x)      </screen>
527         <para>
528             where <literal>x</literal> is a vector and <literal>f(x)</literal>
529             is a function that returns a scalar. This function can also solve bound
530             constrained nonlinear optimization problems:
531         </para>
532         <screen>min f(x)
533             binf &lt;= x &lt;= bsup
534         </screen>
535         <para>
536             where <literal>binf</literal> is the lower bound and
537             <literal>bsup</literal> is the upper bound on <literal>x</literal>.
538         </para>
539         <para>
540             The <literal>costf</literal> argument can be a Scilab function, a
541             list or a string giving the name of a C or Fortran routine (see
542             "external"). This external must return the value <literal>f</literal> of
543             the cost function at the point <literal>x</literal> and the gradient
544             <literal>g</literal> of the cost function at the point
545             <literal>x</literal>.
546         </para>
547         <variablelist>
548             <varlistentry>
549                 <term>Scilab function case</term>
550                 <listitem>
551                     <para>
552                         If <literal>costf</literal> is a Scilab function, its calling
553                         sequence must be:
554                     </para>
555                     <screen>[f, g, ind] = costf(x, ind)      </screen>
556                     <para>
557                         where <literal>x</literal> is the current point,
558                         <literal>ind</literal> is an integer flag described below,
559                         <literal>f</literal> is the real value of the objective function at
560                         the point <literal>x</literal> and <literal>g</literal> is a vector
561                         containing the gradient of the objective function at
562                         <literal>x</literal>. The variable <literal>ind</literal> is
563                         described below.
564                     </para>
565                 </listitem>
566             </varlistentry>
567             <varlistentry>
568                 <term>List case</term>
569                 <listitem>
570                     <para>It may happen that objective function requires extra
571                         arguments. In this case, we can use the following feature. The
572                         <literal>costf</literal> argument can be the list
573                         <literal>(real_costf, arg1,...,argn)</literal>. In this case,
574                         <literal>real_costf</literal>, the first element in the list, must
575                         be a Scilab function with calling sequence: <screen>        [f,g,ind]=real_costf(x,ind,arg1,...,argn)      </screen>
576                         The <literal>x</literal>, <literal>f</literal>,
577                         <literal>g</literal>, <literal>ind</literal> arguments have the same
578                         meaning as before. In this case, each time the objective function is
579                         called back, the arguments <literal>arg1,...,argn</literal> are
580                         automatically appended at the end of the calling sequence of
581                         <literal>real_costf</literal>.
582                     </para>
583                 </listitem>
584             </varlistentry>
585             <varlistentry>
586                 <term>String case</term>
587                 <listitem>
588                     <para>
589                         If <literal>costf</literal> is a string, it refers to the name
590                         of a C or Fortran routine which must be linked to Scilab
591                     </para>
592                     <variablelist>
593                         <varlistentry>
594                             <term>Fortran case</term>
595                             <listitem>
596                                 <para>The calling sequence of the Fortran subroutine computing
597                                     the objective must be:
598                                 </para>
599                                 <screen>subroutine costf(ind,n,x,f,g,ti,tr,td)      </screen>
600                                 <para>with the following declarations:</para>
601                                 <screen>integer ind,n ti(*)
602                                     double precision x(n),f,g(n),td(*)
603                                     real tr(*)
604                                 </screen>
605                                 <para>
606                                     The argument <literal>ind</literal> is described
607                                     below.
608                                 </para>
609                                 <para>If ind = 2, 3 or 4, the inputs of the routine are :
610                                     <literal>x, ind, n, ti, tr,td</literal>.
611                                 </para>
612                                 <para>If ind = 2, 3 or 4, the outputs of the routine are :
613                                     <literal>f</literal> and <literal>g</literal>.
614                                 </para>
615                             </listitem>
616                         </varlistentry>
617                         <varlistentry>
618                             <term>C case</term>
619                             <listitem>
620                                 <para>The calling sequence of the C function computing the
621                                     objective must be:
622                                 </para>
623                                 <screen>void costf(int *ind, int *n, double *x, double *f, double *g, int *ti, float *tr, double *td)      </screen>
624                                 <para>
625                                     The argument <literal>ind</literal> is described
626                                     below.
627                                 </para>
628                                 <para>The inputs and outputs of the function are the same as
629                                     in the fortran case.
630                                 </para>
631                             </listitem>
632                         </varlistentry>
633                     </variablelist>
634                 </listitem>
635             </varlistentry>
636         </variablelist>
637         <para>
638             On output, <literal>ind&lt;0</literal> means that
639             <literal>f</literal> cannot be evaluated at <literal>x</literal> and
640             <literal>ind=0</literal> interrupts the optimization.
641         </para>
642     </refsection>
643     <refsection>
644         <title>Termination criteria</title>
645         <para>Each algorithm has its own termination criteria, which may use the
646             parameters given by the user, that is <literal>nap</literal>,
647             <literal>iter</literal>, <literal>epsg</literal>, <literal>epsf</literal>
648             and <literal>epsx</literal>. Not all the parameters are taken into
649             account. In the table below, we present the specific termination
650             parameters which are taken into account by each algorithm. The
651             unconstrained solver is identified by "UNC" while the bound constrained
652             solver is identified by "BND". An empty entry means that the parameter is
653             ignored by the algorithm.
654         </para>
655         <para>
656             <informaltable border="1">
657                 <tr>
658                     <td>Solver</td>
659                     <td>nap</td>
660                     <td>iter</td>
661                     <td>epsg</td>
662                     <td>epsf</td>
663                     <td>epsx</td>
664                 </tr>
665                 <tr>
666                     <td>optim/"qn" UNC</td>
667                     <td>X</td>
668                     <td>X</td>
669                     <td>X</td>
670                     <td/>
671                     <td/>
672                 </tr>
673                 <tr>
674                     <td>optim/"qn" BND</td>
675                     <td>X</td>
676                     <td>X</td>
677                     <td>X</td>
678                     <td>X</td>
679                     <td>X</td>
680                 </tr>
681                 <tr>
682                     <td>optim/"gc" UNC</td>
683                     <td>X</td>
684                     <td>X</td>
685                     <td>X</td>
686                     <td/>
687                     <td>X</td>
688                 </tr>
689                 <tr>
690                     <td>optim/"gc" BND</td>
691                     <td>X</td>
692                     <td>X</td>
693                     <td>X</td>
694                     <td>X</td>
695                     <td>X</td>
696                 </tr>
697                 <tr>
698                     <td>optim/"nd" UNC</td>
699                     <td>X</td>
700                     <td>X</td>
701                     <td/>
702                     <td>X</td>
703                     <td>X</td>
704                 </tr>
705             </informaltable>
706         </para>
707     </refsection>
708     <refsection>
709         <title>Example: Scilab function</title>
710         <para>The following is an example with a Scilab function. Notice, for
711             simplifications reasons, the Scilab function "cost" of the following
712             example computes the objective function f and its derivative no matter of
713             the value of ind. This allows to keep the example simple. In practical
714             situations though, the computation of "f" and "g" may raise performances
715             issues so that a direct optimization may be to use the value of "ind" to
716             compute "f" and "g" only when needed.
717         </para>
718         <programlisting role="example">function [f, g, ind] = cost(x, ind)
719             xref = [1; 2; 3];
720             f = 0.5 * norm(x - xref)^2;
721             g = x - xref;
722             endfunction
723             
724             // Simplest call
725             x0 = [1; -1; 1];
726             [fopt, xopt] = optim(cost, x0)
727             
728             // Use "gc" algorithm
729             [fopt, xopt, gopt] = optim(cost, x0, "gc")
730             
731             // Use "nd" algorithm
732             [fopt, xopt, gopt] = optim(cost, x0, "nd")
733             
734             // Upper and lower bounds on x
735             [fopt, xopt, gopt] = optim(cost, "b", [-1;0;2], [0.5;1;4], x0)
736             
737             // Upper and lower bounds on x and setting up the algorithm to "gc"
738             [fopt, xopt, gopt] = optim(cost, "b", [-1; 0; 2], [0.5; 1; 4], x0, "gc")
739             
740             // Bound on the number of calls to the objective function
741             [fopt, xopt, gopt] = optim(cost, "b", [-1; 0; 2], [0.5; 1; 4], x0, "gc", "ar", 3)
742             
743             // Set max number of calls to the objective function (3)
744             // Set max number of iterations (100)
745             // Set stopping threshold on the value of f (1e-6),
746             // on the value of the norm of the gradient of the objective function (1e-6)
747             // on the improvement on the parameters x_opt (1e-6;1e-6;1e-6)
748             [fopt, 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])
749             
750             // Additional messages are printed in the console.
751             [fopt, xopt] = optim(cost, x0, imp = 3)
752         </programlisting>
753     </refsection>
754     <refsection>
755         <title>Example: Print messages</title>
756         <para>
757             The <literal>imp</literal> flag may take negative integer values,
758             say k. In that case, the cost function is called once every -k iterations.
759             This allows to draw the function value or write a log file.
760         </para>
761         <para>
762             This feature is available only with the <literal>"qn"</literal>
763             algorithm without constraints.
764         </para>
765         <para>In the following example, we solve the Rosenbrock test case. For
766             each iteration of the algorithm, we print the value of x, f and g.
767         </para>
768         <programlisting role="example">function [f, g, ind] = cost(x, ind)
769             xref = [1; 2; 3];
770             f = 0.5 * norm(x - xref)^2;
771             g = x - xref;
772             if (ind == 1) then
773             mprintf("f(x) = %s, |g(x)|=%s\n", string(f), string(norm(g)))
774             end
775             endfunction
776             
777             x0 = [1; -1; 1];
778             [fopt, xopt] = optim(cost, x0, imp = -1)
779         </programlisting>
780         <para>The previous script produces the following output.</para>
781         <screen>--&gt;[fopt, xopt] = optim(cost, x0, imp = -1)
782             f(x) = 6.5, |g(x)|=3.6055513
783             f(x) = 2.8888889, |g(x)|=2.4037009
784             f(x) = 9.861D-31, |g(x)|=1.404D-15
785             f(x) = 0, |g(x)|=0
786             Norm of projected gradient lower than   0.0000000D+00.
787             xopt  =
788             1.
789             2.
790             3.
791             fopt  =
792             0.
793         </screen>
794         <para>In the following example, we solve the Rosenbrock test case. For
795             each iteration of the algorithm, we plot the current value of x into a 2D
796             graph containing the contours of Rosenbrock's function. This allows to see
797             the progress of the algorithm while the algorithm is performing. We could
798             as well write the value of x, f and g into a log file if needed.
799         </para>
800         <programlisting role="example">// 1. Define Rosenbrock for optimization
801             function [f , g , ind] = rosenbrock (x , ind)
802             f = 100.0 *(x(2) - x(1)^2)^2 + (1 - x(1))^2;
803             g(1) = - 400. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
804             g(2) = 200. * ( x(2) - x(1)**2 )
805             endfunction
806             
807             // 2. Define rosenbrock for contouring
808             function f = rosenbrockC ( x1 , x2 )
809             x = [x1 x2]
810             ind = 4
811             [ f , g , ind ] = rosenbrock ( x , ind )
812             endfunction
813             
814             // 3. Define Rosenbrock for plotting
815             function [ f , g , ind ] = rosenbrockPlot ( x , ind )
816             [ f , g , ind ] = rosenbrock ( x , ind )
817             if (ind == 1) then
818             plot ( x(1) , x(2) , "g." )
819             end
820             endfunction
821             
822             // 4. Draw the contour of Rosenbrock's function
823             x0 = [-1.2 1.0];
824             xopt = [1.0 1.0];
825             xdata = linspace(-2,2,100);
826             ydata = linspace(-2,2,100);
827             contour ( xdata , ydata , rosenbrockC , [1 10 100 500 1000])
828             plot(x0(1) , x0(2) , "b.")
829             plot(xopt(1) , xopt(2) , "r*")
830             
831             // 5. Plot the optimization process, during optimization
832             [fopt, xopt] = optim ( rosenbrockPlot , x0 , imp = -1)
833         </programlisting>
834         <scilab:image>
835             function [f, g, ind]=rosenbrock(x, ind)
836             f = 100.0 *(x(2) - x(1)^2)^2 + (1 - x(1))^2;
837             g(1) = - 400. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
838             g(2) = 200. * ( x(2) - x(1)**2 )
839             endfunction
840             
841             function f=rosenbrockC(x1, x2)
842             x = [x1 x2]
843             ind = 4
844             [ f , g , ind ] = rosenbrock ( x , ind )
845             endfunction
846             
847             function [f, g, ind]=rosenbrockPlot(x, ind)
848             [ f , g , ind ] = rosenbrock ( x , ind )
849             if (ind == 1) then
850             plot ( x(1) , x(2) , "g." )
851             end
852             endfunction
853             
854             x0 = [-1.2 1.0];
855             xopt = [1.0 1.0];
856             xdata = linspace(-2,2,100);
857             ydata = linspace(-2,2,100);
858             contour ( xdata , ydata , rosenbrockC , [1 10 100 500 1000])
859             plot(x0(1) , x0(2) , "b.")
860             plot(xopt(1) , xopt(2) , "r*")
861             [fopt, xopt] = optim ( rosenbrockPlot , x0 , imp = -1)
862         </scilab:image>
863         
864     </refsection>
865     <refsection>
866         <title>Example: Optimizing with numerical derivatives</title>
867         <para>It is possible to optimize a problem without an explicit knowledge
868             of the derivative of the cost function. For this purpose, we can use the
869             numdiff or derivative function to compute a numerical derivative of the
870             cost function.
871         </para>
872         <para>In the following example, we use the numdiff function to solve
873             Rosenbrock's problem.
874         </para>
875         <programlisting role="example">function f = rosenbrock ( x )
876             f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
877             endfunction
878             
879             function [ f , g , ind ] = rosenbrockCost ( x , ind )
880             f = rosenbrock ( x );
881             g= numdiff ( rosenbrock , x );
882             endfunction
883             
884             x0 = [-1.2 1.0];
885             
886             [ fopt , xopt ] = optim ( rosenbrockCost , x0 )
887         </programlisting>
888         <para>In the following example, we use the derivative function to solve
889             Rosenbrock's problem. Given that the step computation strategy is not the
890             same in numdiff and derivative, this might lead to improved
891             results.
892         </para>
893         <programlisting role="example">function f = rosenbrock ( x )
894             f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
895             endfunction
896             
897             function [ f , g , ind ] = rosenbrockCost2 ( x , ind )
898             f = rosenbrock ( x );
899             g = derivative ( rosenbrock , x.' , order = 4 );
900             endfunction
901             
902             x0 = [-1.2 1.0];
903             [fopt , xopt] = optim ( rosenbrockCost2 , x0 )
904         </programlisting>
905     </refsection>
906     <refsection>
907         <title>Example: Counting function evaluations and number of
908             iterations
909         </title>
910         <para>
911             The <literal>imp</literal> option can take negative values. If the
912             <literal>imp</literal> is equal to <literal>m</literal> where
913             <literal>m</literal> is a negative integer, then the cost function is
914             evaluated every -<literal>m</literal> iterations, with the
915             <literal>ind</literal> input argument equal to 1. The following example
916             uses this feature to compute the number of iterations. The global variable
917             <literal>mydata</literal> is used to store the number of function
918             evaluations as well as the number of iterations.
919         </para>
920         <programlisting role="example"><![CDATA[
921 function [f, g, ind] = cost(x, ind)
922     global _MYDATA_
923     if ( ind == 1 )
924         _MYDATA_.niter = _MYDATA_.niter + 1;
925     end
926     _MYDATA_.nfevals = _MYDATA_.nfevals + 1;
927     xref = [1; 2; 3];
928     if ( ind == 2 | ind == 4 ) then
929         f = 0.5*norm(x-xref)^2;
930     else
931         f = 0;
932     end
933     if ( ind == 3 | ind == 4 ) then
934         g = x-xref;
935     else
936         g = zeros(3, 1);
937     end
938 endfunction
939
940 x0 = [1; -1; 1];
941 global _MYDATA_
942 _MYDATA_ = tlist ( ["MYDATA", "niter", "nfevals"]);
943 _MYDATA_.niter = 0;
944 _MYDATA_.nfevals = 0;
945
946 [f, xopt] = optim(cost, x0, imp=-1);
947 mprintf ( "Number of function evaluations: %d\n", _MYDATA_.nfevals );
948 mprintf ( "Number of iterations: %d\n", _MYDATA_.niter );
949         ]]></programlisting>
950         <para>While the previous example perfectly works, there is a risk that the
951             same variable <literal>_MYDATA_</literal> is used by some internal
952             function used by <literal>optim</literal>. In this case, the value may be
953             wrong. This is why a sufficiently weird variable name has been
954             used.
955         </para>
956     </refsection>
957     <refsection>
958         <title>Example : Passing extra parameters</title>
959         <para>In most practical situations, the cost function depends on extra
960             parameters which are required to evaluate the cost function. There are
961             several methods to achieve this goal.
962         </para>
963         <para>In the following example, the cost function uses 4 parameters
964             <literal>a, b, c</literal> and <literal>d</literal>. We define the cost
965             function with additional input arguments, which are declared after the
966             index argument. Then we pass a list as the first input argument of the
967             <literal>optim</literal> solver. The first element of the list is the cost
968             function. The additional variables are directly passed to the cost
969             function.
970         </para>
971         <programlisting role="example">function [ f , g , ind ] = costfunction ( x , ind , a , b , c , d )
972             f = a * ( x(1) - c ) ^2 + b * ( x(2) - d )^2
973             g(1) = 2 * a * ( x(1) - c )
974             g(2) = 2 * b * ( x(2) - d )
975             endfunction
976             
977             x0 = [1 1];
978             a = 1.0;
979             b = 2.0;
980             c = 3.0;
981             d = 4.0;
982             costf = list ( costfunction , a , b , c, d );
983             [fopt , xopt] = optim ( costf , x0 , imp = 2)
984         </programlisting>
985         <para>In complex cases, the cost function may have so many parameters,
986             that having a function which takes all arguments as inputs is not
987             convenient. For example, consider the situation where the cost function
988             needs 12 parameters. Then, designing a function with 14 input arguments
989             (x, index and the 12 parameters) is difficult to manage. Instead, we can
990             use a more complex data structure to store our data. In the following
991             example, we use a tlist to store the 4 input arguments. This method can
992             easily be expanded to an arbitrary number of parameters.
993         </para>
994         <programlisting role="example">function [f , g , ind] = costfunction ( x , ind , parameters)
995             // Get the parameters
996             a = parameters.a
997             b = parameters.b
998             c = parameters.c
999             d = parameters.d
1000             f = a * ( x(1) - c ) ^2 + b * ( x(2) - d )^2
1001             g(1) = 2 * a * ( x(1) - c )
1002             g(2) = 2 * b * ( x(2) - d )
1003             endfunction
1004             
1005             x0 = [1 1];
1006             a = 1.0;
1007             b = 2.0;
1008             c = 3.0;
1009             d = 4.0;
1010             // Store the parameters
1011             parameters = tlist ( [
1012             "T_MYPARAMS"
1013             "a"
1014             "b"
1015             "c"
1016             "d"
1017             ]);
1018             
1019             parameters.a = a;
1020             parameters.b = b;
1021             parameters.c = c;
1022             parameters.d = d;
1023             costf = list ( costfunction , parameters );
1024             [fopt , xopt] = optim ( costf , x0 , imp = 2)
1025         </programlisting>
1026         <para>In the following example, the parameters are defined before the
1027             optimizer is called. They are directly used in the cost function.
1028         </para>
1029         <programlisting role="example">// The example NOT to follow
1030             function [ f , g , ind ] = costfunction ( x , ind )
1031             f = a * ( x(1) - c ) ^2 + b * ( x(2) - d )^2
1032             g(1) = 2 * a * ( x(1) - c )
1033             g(2) = 2 * b * ( x(2) - d )
1034             endfunction
1035             x0 = [1 1];
1036             a = 1.0;
1037             b = 2.0;
1038             c = 3.0;
1039             d = 4.0;
1040             [ fopt , xopt ] = optim ( costfunction , x0 , imp = 2 )
1041         </programlisting>
1042         <para>While the previous example perfectly works, there is a risk that the
1043             same variables are used by some internal function used by
1044             <literal>optim</literal>. In this case, the value of the parameters are
1045             not what is expected and the optimization can fail or, worse, give a wrong
1046             result. It is also difficult to manage such a function, which requires
1047             that all the parameters are defined in the calling context.
1048         </para>
1049         <para>In the following example, we define the cost function with the
1050             classical header. Inside the function definition, we declare that the
1051             parameters <literal>a, b, c</literal> and <literal>d</literal> are global
1052             variables. Then we declare and set the global variables.
1053         </para>
1054         <programlisting role="example">// Another example NOT to follow
1055             function [ f , g , ind ] = costfunction ( x , ind )
1056             global a b c d
1057             f = a * ( x(1) - c ) ^2 + b * ( x(2) - d )^2
1058             g(1) = 2 * a * ( x(1) - c )
1059             g(2) = 2 * b * ( x(2) - d )
1060             endfunction
1061             x0 = [1 1];
1062             global a b c d
1063             a = 1.0;
1064             b = 2.0;
1065             c = 3.0;
1066             d = 4.0;
1067             [ fopt , xopt ] = optim ( costfunction , x0 , imp = 2 )
1068         </programlisting>
1069         <para>While the previous example perfectly works, there is a risk that the
1070             same variables are used by some internal function used by
1071             <literal>optim</literal>. In this case, the value of the parameters are
1072             not what is expected and the optimization can fail or, worse, give a wrong
1073             result.
1074         </para>
1075     </refsection>
1076     <refsection>
1077         <title>Example : Checking that derivatives are correct</title>
1078         <para>Many optimization problem can be avoided if the derivatives are
1079             computed correctly. One common reason for failure in the step-length
1080             procedure is an error in the calculation of the cost function and its
1081             gradient. Incorrect calculation of derivatives is by far the most common
1082             user error.
1083         </para>
1084         <para>In the following example, we give a false implementation of
1085             Rosenbrock's gradient. In order to check the computation of the
1086             derivatives, we use the <literal>derivative</literal> function. We define
1087             the <literal>simplified</literal> function, which delegates the
1088             computation of <literal>f</literal> to the rosenbrock function. The
1089             <literal>simplified</literal> function is passed as an input argument of
1090             the <literal>derivative</literal> function.
1091         </para>
1092         <programlisting role="example">function [ f , g , index ] = rosenbrock ( x , index )
1093             f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
1094             // Exact :
1095             g(1) = - 400. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
1096             // Wrong :
1097             g(1) = - 1200. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
1098             g(2) = 200. * ( x(2) - x(1)**2 )
1099             endfunction
1100             
1101             function f = simplified ( x )
1102             index = 1;
1103             [ f , g , index ] = rosenbrock ( x , index )
1104             endfunction
1105             
1106             x0 = [-1.2 1];
1107             index = 1;
1108             [ f , g , index ] = rosenbrock ( x0 , index );
1109             gnd = derivative ( simplified , x0.' );
1110             mprintf("Exact derivative:[%s]\n" , strcat ( string(g) , " " ));
1111             mprintf("Numerical derivative:[%s]\n" , strcat ( string(gnd) , " " ));
1112         </programlisting>
1113         <para>The previous script produces the following output. Obviously, the
1114             difference between the two gradient is enormous, which shows that the
1115             wrong formula has been used in the gradient.
1116         </para>
1117         <programlisting role="example">      Exact derivative:[-638 -88]
1118             Numerical derivative:[-215.6 -88]
1119         </programlisting>
1120     </refsection>
1121     <refsection>
1122         <title>Example: C function</title>
1123         <para>The following is an example with a C function, where a C source code
1124             is written into a file, dynamically compiled and loaded into Scilab, and
1125             then used by the "optim" solver. The interface of the "rosenc" function is
1126             fixed, even if the arguments are not really used in the cost function.
1127             This is because the underlying optimization solvers must assume that the
1128             objective function has a known, constant interface. In the following
1129             example, the arrays ti and tr are not used, only the array "td" is used,
1130             as a parameter of the Rosenbrock function. Notice that the content of the
1131             arrays ti and td are the same that the content of the Scilab variable, as
1132             expected.
1133         </para>
1134         <programlisting role="example">// External function written in C (C compiler required)
1135             // write down the C code (Rosenbrock problem)
1136             C=['#include &lt;math.h&gt;'
1137             'double sq(double x)'
1138             '{ return x*x;}'
1139             'void rosenc(int *ind, int *n, double *x, double *f, double *g, '
1140             '                                int *ti, float *tr, double *td)'
1141             '{'
1142             '  double p;'
1143             '  int i;'
1144             '  p=td[0];'
1145             '  if (*ind==2||*ind==4) {'
1146             '    *f=1.0;'
1147             '    for (i=1;i&lt;*n;i++)'
1148             '      *f+=p*sq(x[i]-sq(x[i-1]))+sq(1.0-x[i]);'
1149             '  }'
1150             '  if (*ind==3||*ind==4) {'
1151             '    g[0]=-4.0*p*(x[1]-sq(x[0]))*x[0];'
1152             '    for (i=1;i&lt;*n-1;i++)'
1153             '      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]);'
1154             '    g[*n-1]=2.0*p*(x[*n-1]-sq(x[*n-2]))-2.0*(1.0-x[*n-1]);'
1155             '  }'
1156             '}'];
1157             cd TMPDIR;
1158             mputl(C, TMPDIR+'/rosenc.c')
1159             
1160             // compile the C code
1161             l = ilib_for_link('rosenc', 'rosenc.c', [], 'c');
1162             
1163             // incremental linking
1164             link(l, 'rosenc', 'c')
1165             
1166             //solve the problem
1167             x0 = [40; 10; 50];
1168             p = 100;
1169             [f, xo, go] = optim('rosenc', x0, 'td', p)
1170         </programlisting>
1171     </refsection>
1172     <refsection>
1173         <title>Example: Fortran function</title>
1174         <para>The following is an example with a Fortran function.</para>
1175         <programlisting role="example">// External function written in Fortran (Fortran compiler required)
1176             // write down the Fortran  code (Rosenbrock problem)
1177             F = [ '      subroutine rosenf(ind, n, x, f, g, ti, tr, td)'
1178             '      integer ind,n,ti(*)'
1179             '      double precision x(n),f,g(n),td(*)'
1180             '      real tr(*)'
1181             'c'
1182             '      double precision y,p'
1183             '      p=td(1)'
1184             '      if (ind.eq.2.or.ind.eq.4) then'
1185             '        f=1.0d0'
1186             '        do i=2,n'
1187             '          f=f+p*(x(i)-x(i-1)**2)**2+(1.0d0-x(i))**2'
1188             '        enddo'
1189             '      endif'
1190             '      if (ind.eq.3.or.ind.eq.4) then'
1191             '        g(1)=-4.0d0*p*(x(2)-x(1)**2)*x(1)'
1192             '        if(n.gt.2) then'
1193             '          do i=2,n-1'
1194             '            g(i)=2.0d0*p*(x(i)-x(i-1)**2)-4.0d0*p*(x(i+1)-x(i)**2)*x(i)'
1195             '     &amp;           -2.0d0*(1.0d0-x(i))'
1196             '          enddo'
1197             '        endif'
1198             '        g(n)=2.0d0*p*(x(n)-x(n-1)**2)-2.0d0*(1.0d0-x(n))'
1199             '      endif'
1200             '      return'
1201             '      end'];
1202             cd TMPDIR;
1203             mputl(F, TMPDIR+'/rosenf.f')
1204             
1205             // compile the Fortran code
1206             l = ilib_for_link('rosenf', 'rosenf.f', [], 'f');
1207             
1208             // incremental linking
1209             link(l, 'rosenf', 'f')
1210             
1211             //solve the problem
1212             x0 = [40; 10; 50];
1213             p = 100;
1214             [f, xo, go] = optim('rosenf', x0, 'td', p)
1215         </programlisting>
1216     </refsection>
1217     <refsection>
1218         <title>Example: Fortran function with initialization</title>
1219         <para>The following is an example with a Fortran function in which the
1220             "in" option is used to allocate memory inside the Scilab environment. In
1221             this mode, there is a dialog between Scilab and the objective function.
1222             The goal of this dialog is to initialize the parameters of the objective
1223             function. Each part of this dialog is based on a specific value of the
1224             "ind" parameter.
1225         </para>
1226         <para>At the beginning, Scilab calls the objective function, with the ind
1227             parameter equal to 10. This tells the objective function to initialize
1228             the sizes of the arrays it needs by setting the nizs, nrzs and ndzs
1229             integer parameters of the "nird" common. Then the objective function
1230             returns. At this point, Scilab creates internal variables and allocate
1231             memory for the variable izs, rzs and dzs. Scilab calls the objective
1232             function back again, this time with ind equal to 11. This tells the
1233             objective function to initialize the arrays izs, rzs and dzs. When the
1234             objective function has done so, it returns. Then Scilab enters in the real
1235             optimization mode and calls the optimization solver the user requested.
1236             Whenever the objective function is called, the izs, rzs and dzs arrays
1237             have the values that have been previously initialized.
1238         </para>
1239         <programlisting role="example">//
1240             // Define a fortran source code and compile it (fortran compiler required)
1241             //
1242             fortransource = ['      subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
1243             'C     -------------------------------------------'
1244             'c     Example of cost function given by a subroutine'
1245             'c     if n&lt;=2 returns ind=0'
1246             'c     f.bonnans, oct 86'
1247             '      implicit double precision (a-h,o-z)'
1248             '      real rzs(1)'
1249             '      double precision dzs(*)'
1250             '      dimension x(n),g(n),izs(*)'
1251             '      common/nird/nizs,nrzs,ndzs'
1252             '      if (n.lt.3) then'
1253             '        ind=0'
1254             '        return'
1255             '      endif'
1256             '      if(ind.eq.10) then'
1257             '         nizs=2'
1258             '         nrzs=1'
1259             '         ndzs=1'
1260             '         return'
1261             '      endif'
1262             '      if(ind.eq.11) then'
1263             '         izs(1)=5'
1264             '         izs(2)=10'
1265             '         dzs(1)=100.0d+0'
1266             '         return'
1267             '      endif'
1268             '      if(ind.eq.2)go to 5'
1269             '      if(ind.eq.3)go to 20'
1270             '      if(ind.eq.4)go to 5'
1271             '      ind=-1'
1272             '      return'
1273             '5     f=1.0d+0'
1274             '      do 10 i=2,n'
1275             '        im1=i-1'
1276             '10      f=f + dzs(1)*(x(i)-x(im1)**2)**2 + (1.0d+0-x(i))**2'
1277             '      if(ind.eq.2)return'
1278             '20    g(1)=-4.0d+0*dzs(1)*(x(2)-x(1)**2)*x(1)'
1279             '      nm1=n-1'
1280             '      do 30 i=2,nm1'
1281             '        im1=i-1'
1282             '        ip1=i+1'
1283             '        g(i)=2.0d+0*dzs(1)*(x(i)-x(im1)**2)'
1284             '30      g(i)=g(i) -4.0d+0*dzs(1)*(x(ip1)-x(i)**2)*x(i) - '
1285             '     &amp;        2.0d+0*(1.0d+0-x(i))'
1286             '      g(n)=2.0d+0*dzs(1)*(x(n)-x(nm1)**2) - 2.0d+0*(1.0d+0-x(n))'
1287             '      return'
1288             '      end'];
1289             cd TMPDIR;
1290             mputl(fortransource, TMPDIR + '/rosenf.f')
1291             
1292             // compile the C code
1293             libpath = ilib_for_link('rosenf', 'rosenf.f', [], 'f');
1294             
1295             // incremental linking
1296             linkid = link(libpath, 'rosenf', 'f');
1297             
1298             x0 = 1.2 * ones(1, 5);
1299             //
1300             // Solve the problem
1301             //
1302             [f, x, g] = optim('rosenf', x0, 'in');
1303         </programlisting>
1304     </refsection>
1305     <refsection>
1306         <title>Example: Fortran function with initialization on Windows with Intel
1307             Fortran Compiler
1308         </title>
1309         <para>Under the Windows operating system with Intel Fortran Compiler, one
1310             must carefully design the fortran source code so that the dynamic link
1311             works properly. On Scilab's side, the optimization component is
1312             dynamically linked and the symbol "nird" is exported out of the
1313             optimization dll. On the cost function's side, which is also dynamically
1314             linked, the "nird" common must be imported in the cost function
1315             dll.
1316         </para>
1317         <para>The following example is a re-writing of the previous example, with
1318             special attention for the Windows operating system with Intel Fortran
1319             compiler as example. In that case, we introduce additional compiling
1320             instructions, which allows the compiler to import the "nird"
1321             symbol.
1322         </para>
1323         <programlisting role="example">fortransource = ['subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
1324             'cDEC$ IF DEFINED (FORDLL)'
1325             'cDEC$ ATTRIBUTES DLLIMPORT:: /nird/'
1326             'cDEC$ ENDIF'
1327             'C     -------------------------------------------'
1328             'c     Example of cost function given by a subroutine'
1329             'c     if n&lt;=2 returns ind=0'
1330             'c     f.bonnans, oct 86'
1331             '      implicit double precision (a-h,o-z)'
1332             [etc...]
1333         </programlisting>
1334     </refsection>
1335     <refsection>
1336         <title>Example: Fortran function with tab saving
1337         </title>
1338         <para>
1339             Adding the "si" and/or "sd" options to <literal>optim</literal> implies
1340             the addition of the output arguments "ti" and/or "td", which will represent
1341             the working tables of the objective function.
1342         </para>
1343         <para>
1344             Those output arguments must be placed at the end of the output list.
1345         </para>
1346         <programlisting role="example">
1347             // Move into the temporary directory to create the temporary files there
1348             cur_dir = pwd();
1349             chdir(TMPDIR);
1350             fortransource = ['      subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
1351             'C     -------------------------------------------'
1352             'c (DLL Digital Visual Fortran)'
1353             'c On Windows , we need to import common nird from scilab'
1354             'cDEC$ IF DEFINED (FORDLL)'
1355             'cDEC$ ATTRIBUTES DLLIMPORT:: /nird/'
1356             'cDEC$ ENDIF'
1357             'C     -------------------------------------------'
1358             'c     Example of cost function given by a subroutine'
1359             'c     if n.le.2 returns ind=0'
1360             'c     f.bonnans, oct 86'
1361             '      implicit double precision (a-h,o-z)'
1362             '      real rzs(1)'
1363             '      double precision dzs(*)'
1364             '      dimension x(n),g(n),izs(*)'
1365             '      common/nird/nizs,nrzs,ndzs'
1366             '      if (n.lt.3) then'
1367             '        ind=0'
1368             '        return'
1369             '      endif'
1370             '      if(ind.eq.10) then'
1371             '         nizs=2'
1372             '         nrzs=1'
1373             '         ndzs=1'
1374             '         return'
1375             '      endif'
1376             '      if(ind.eq.11) then'
1377             '         izs(1)=5'
1378             '         izs(2)=10'
1379             '         dzs(1)=100.0d+0'
1380             '         return'
1381             '      endif'
1382             '      if(ind.eq.2)go to 5'
1383             '      if(ind.eq.3)go to 20'
1384             '      if(ind.eq.4)go to 5'
1385             '      ind=-1'
1386             '      return'
1387             '5     f=1.0d+0'
1388             '      do 10 i=2,n'
1389             '        im1=i-1'
1390             '10      f=f + dzs(1)*(x(i)-x(im1)**2)**2 + (1.0d+0-x(i))**2'
1391             '      if(ind.eq.2)return'
1392             '20    g(1)=-4.0d+0*dzs(1)*(x(2)-x(1)**2)*x(1)'
1393             '      nm1=n-1'
1394             '      do 30 i=2,nm1'
1395             '        im1=i-1'
1396             '        ip1=i+1'
1397             '        g(i)=2.0d+0*dzs(1)*(x(i)-x(im1)**2)'
1398             '30      g(i)=g(i) -4.0d+0*dzs(1)*(x(ip1)-x(i)**2)*x(i) - '
1399             '     &amp;        2.0d+0*(1.0d+0-x(i))'
1400             '      g(n)=2.0d+0*dzs(1)*(x(n)-x(nm1)**2) - 2.0d+0*(1.0d+0-x(n))'
1401             '      return'
1402             '      end'];
1403             mputl(fortransource, TMPDIR + '/rosenf.f');
1404             ilib_for_link('rosenf', 'rosenf.f', [], 'f');
1405             exec loader.sce;
1406             chdir(cur_dir);
1407             //
1408             // Define some constants
1409             //
1410             Leps = 10e3 * 8.e-5;
1411             bs = 10.*ones(1, 5);
1412             bi = -bs;
1413             x0 = 0.12 * bs;
1414             epsx = 1.e-15 * x0;
1415             xopt = .1*bs;
1416             
1417             // 'ti' and 'td' always at the end of the output sequence
1418             [f, x, ti, td]                           = optim('rosenf', x0, 'in', 'si', 'sd')
1419             [f, x, g, ti, td]                        = optim('rosenf', x0, 'in', 'si', 'sd')
1420             [f, x, g, w, ti, td]                     = optim('rosenf', x0, 'in', 'si', 'sd')
1421             [f, x, g, w, niter, nevals, ti, td]      = optim('rosenf', x0, 'in', 'si', 'sd')
1422             [f, x, g, w, niter, nevals, err, ti, td] = optim('rosenf', x0, 'in', 'si', 'sd')
1423             
1424             // With input argument 'in', ti and td will be initialized by rosenf function.
1425             [f, x, ti, td] = optim('rosenf', x0, 'in', 'si', 'sd')
1426             // Reuses the last ti and td for the next call and return it again.
1427             [f, x, ti, td] = optim('rosenf', x0, 'ti', ti, 'td', td, 'si', 'sd')
1428             
1429             // Initializes ti and td but return only ti
1430             [f, x, ti] = optim('rosenf', x0, 'in', 'si')
1431             [f, x, ti] = optim('rosenf', x0, 'ti', ti, 'si')
1432             
1433             // Initializes ti and td but return only td
1434             [f, x, td] = optim('rosenf', x0, 'in', 'sd')
1435             [f, x, td] = optim('rosenf', x0, 'td', td, 'sd')
1436         </programlisting>
1437     </refsection>
1438     <refsection role="see also">
1439         <title>See Also</title>
1440         <simplelist type="inline">
1441             <member>
1442                 <link linkend="external">external</link>
1443             </member>
1444             <member>
1445                 <link linkend="qpsolve">qpsolve</link>
1446             </member>
1447             <member>
1448                 <link linkend="datafit">datafit</link>
1449             </member>
1450             <member>
1451                 <link linkend="leastsq">leastsq</link>
1452             </member>
1453             <member>
1454                 <link linkend="numdiff">numdiff</link>
1455             </member>
1456             <member>
1457                 <link linkend="derivative">derivative</link>
1458             </member>
1459             <member>
1460                 <link linkend="NDcost">NDcost</link>
1461             </member>
1462         </simplelist>
1463     </refsection>
1464     <refsection>
1465         <title>References</title>
1466         <para>The following is a map from the various options to the underlying
1467             solvers.
1468         </para>
1469         <variablelist>
1470             <varlistentry>
1471                 <term>"qn" without constraints</term>
1472                 <listitem>
1473                     <para>n1qn1 : a quasi-Newton method with a Wolfe-type line
1474                         search
1475                     </para>
1476                 </listitem>
1477             </varlistentry>
1478             <varlistentry>
1479                 <term>"qn" with bounds constraints</term>
1480                 <listitem>
1481                     <para>qnbd : a quasi-Newton method with projection</para>
1482                     <para>RR-0242 - A variant of a projected variable metric method for
1483                         bound constrained optimization problems, Bonnans Frederic, Rapport
1484                         de recherche de l'INRIA - Rocquencourt, Octobre 1983
1485                     </para>
1486                 </listitem>
1487             </varlistentry>
1488             <varlistentry>
1489                 <term>"gc" without constraints</term>
1490                 <listitem>
1491                     <para>n1qn3 : a Quasi-Newton limited memory method with BFGS.</para>
1492                 </listitem>
1493             </varlistentry>
1494             <varlistentry>
1495                 <term>"gc" with bounds constraints</term>
1496                 <listitem>
1497                     <para>gcbd : a BFGS-type method with limited memory and
1498                         projection
1499                     </para>
1500                 </listitem>
1501             </varlistentry>
1502             <varlistentry>
1503                 <term>"nd" without constraints</term>
1504                 <listitem>
1505                     <para>n1fc1 : a bundle method</para>
1506                 </listitem>
1507             </varlistentry>
1508             <varlistentry>
1509                 <term>"nd" with bounds constraints</term>
1510                 <listitem>
1511                     <para>not available</para>
1512                 </listitem>
1513             </varlistentry>
1514         </variablelist>
1515     </refsection>
1516 </refentry>