[optimization] windows link fixed && valgrind error fixed
[scilab.git] / scilab / modules / optimization / sci_gateway / cpp / sci_optim.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - Cedric DELAMARRE
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16 #include "optimization_gw.hxx"
17 #include "function.hxx"
18 #include "double.hxx"
19 #include "string.hxx"
20 #include "polynom.hxx"
21 #include "list.hxx"
22 #include "optimizationfunctions.hxx"
23 #include "numericconstants.hxx"
24 #include <limits>
25
26 extern "C"
27 {
28 #include "localization.h"
29 #include "Scierror.h"
30 #include "sciprint.h"
31 #include "common_structure_optimization.h"
32 #include "sci_malloc.h"
33 #include "elem_common.h"
34 #include "scioptimfunctions.h"
35 #include "checkoptimerror.h"
36 }
37 /*--------------------------------------------------------------------------*/
38 types::Function::ReturnValue sci_optim(types::typed_list &in, types::optional_list &opt, int _iRetCount, types::typed_list &out)
39 {
40     types::Function::ReturnValue ret = types::Function::Error;
41     OptimizationFunctions* opFunctionsManager = NULL;
42
43     types::Double* pDblX0   = NULL;
44     types::Double* pDblBinf = NULL;
45     types::Double* pDblBsub = NULL;
46     types::Double* pDblTi   = NULL;
47     types::Double* pDblTd   = NULL;
48     types::Double* pDblNap  = NULL;
49     types::Double* pDblIter = NULL;
50     types::Double* pDblEpsg = NULL;
51     types::Double* pDblEpsf = NULL;
52     types::Double* pDblEpsx = NULL;
53     types::Double* pDblWork = NULL;
54     //    types::Polynom* pPolyX0     = NULL;
55     //    types::Polynom* pPolyBinf   = NULL;
56     //    types::Polynom* pPolyBsub   = NULL;
57
58     int* piIzs          = NULL;
59     int* piWork         = NULL;
60     float* pfRzs        = NULL;
61     double* pdblDzs     = NULL;
62     double* pdblWork    = NULL;
63     double* pdblWork2   = NULL;
64     double* pdblX0      = NULL;
65     double* pdblEpsx    = NULL;
66     double* pdblG       = NULL;
67     double* pdblVar     = NULL;
68     double* pdblBsub    = NULL;
69     double* pdblBinf    = NULL;
70
71     bool bMem       = false;
72     int iEpsx       = 1;
73     int iDzs        = 1;
74     int iPos        = 0;
75     int iContr      = 1;
76     int iSizeX0     = 0;
77     int iSizeBinf   = 0;
78     int iSizeBsub   = 0;
79     int iSim        = 0; // 1 : c function || 2 : macro
80     int iAlgo       = 1; // 1 : qn || 2 : gc || 10 : nd
81     int iMem        = 0;
82     int iGetArgs    = 0;
83     int iIndSim     = 0;
84     int iIndOpt     = 0;
85     int iSaveI      = 0;
86     int iSaveD      = 0;
87     int iArret      = 0;
88     int iMode       = 1;
89     int iWorkSize   = 0;
90     int iWorkSizeI  = 0;
91     int iNitv       = 0;
92     int io          = 0; // not used in scilab 6 and more
93     int iPrint        = 0;
94     int iZero       = 0;
95     int iOne        = 1;
96     double df0      = 1;
97     double dF       = 0;
98
99     // stop arguments "ar"
100     int iNap        = 100;
101     int iItMax      = 100;
102     double dEpsg    = NumericConstants::eps_machine; // p : eps*base
103     double dTol     = dEpsg;
104     double dEpsf    = 0;
105
106     try
107     {
108         if (in.size() < 2 || in.size() > 18)
109         {
110             Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "optim", 2, 18);
111             throw ast::ScilabException();
112         }
113
114         if (_iRetCount > 7)
115         {
116             Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "optim", 1, 7);
117             throw ast::ScilabException();
118         }
119
120         /*** get inputs arguments ***/
121         // get optionals
122         for (const auto & o : opt)
123         {
124             // "iprint"
125             if (o.first == L"iprint")
126             {
127                 if (o.second->isDouble() == false)
128                 {
129                     Scierror(999, _("%s: Wrong type for input argument #%s: A scalar expected.\n"), "optim", "iprint");
130                     throw ast::ScilabException();
131                 }
132
133                 types::Double* pDblPrint = o.second->getAs<types::Double>();
134
135                 if (pDblPrint->isScalar() == false)
136                 {
137                     Scierror(999, _("%s: Wrong type for input argument #%s: A scalar expected.\n"), "optim", "iprint");
138                     throw ast::ScilabException();
139                 }
140
141                 iPrint = (int)pDblPrint->get(0);
142             }
143             // "nap"
144             else if (o.first == L"nap")
145             {
146                 if (o.second->isDouble() == false)
147                 {
148                     Scierror(999, _("%s: Wrong type for input argument #%s: A real scalar expected.\n"), "optim", "nap");
149                     throw ast::ScilabException();
150                 }
151
152                 pDblNap = o.second->getAs<types::Double>();
153                 if (pDblNap->isScalar() == false || pDblNap->isComplex())
154                 {
155                     Scierror(999, _("%s: Wrong size for input argument #%s: A real scalar expected.\n"), "optim", "nap");
156                     throw ast::ScilabException();
157                 }
158
159                 iNap = (int)pDblNap->get(0);
160             }
161             // "iter"
162             else if (o.first == L"iter")
163             {
164                 if (o.second->isDouble() == false)
165                 {
166                     Scierror(999, _("%s: Wrong type for input argument #%s: A real scalar expected.\n"), "optim", "iter");
167                     throw ast::ScilabException();
168                 }
169
170                 pDblIter = o.second->getAs<types::Double>();
171                 if (pDblIter->isScalar() == false || pDblIter->isComplex())
172                 {
173                     Scierror(999, _("%s: Wrong size for input argument #%s: A real scalar expected.\n"), "optim", "iter");
174                     throw ast::ScilabException();
175                 }
176
177                 iItMax = (int)pDblIter->get(0);
178             }
179             // "epsg"
180             else if (o.first == L"epsg")
181             {
182                 if (o.second->isDouble() == false)
183                 {
184                     Scierror(999, _("%s: Wrong type for input argument #%s: A real scalar expected.\n"), "optim", "epsg");
185                     throw ast::ScilabException();
186                 }
187
188                 pDblEpsg = o.second->getAs<types::Double>();
189                 if (pDblEpsg->isScalar() == false || pDblEpsg->isComplex())
190                 {
191                     Scierror(999, _("%s: Wrong size for input argument #%s: A real scalar expected.\n"), "optim", "epsg");
192                     throw ast::ScilabException();
193                 }
194
195                 dEpsg = pDblEpsg->get(0);
196             }
197             // "epsf"
198             else if (o.first == L"epsf")
199             {
200                 if (o.second->isDouble() == false)
201                 {
202                     Scierror(999, _("%s: Wrong type for input argument #%s: A real scalar expected.\n"), "optim", "epsf");
203                     throw ast::ScilabException();
204                 }
205
206                 pDblEpsf = o.second->getAs<types::Double>();
207                 if (pDblEpsf->isScalar() == false || pDblEpsf->isComplex())
208                 {
209                     Scierror(999, _("%s: Wrong size for input argument #%s: A real scalar expected.\n"), "optim", "epsf");
210                     throw ast::ScilabException();
211                 }
212
213                 dEpsf = pDblEpsf->get(0);
214             }
215             // "epsx"
216             else if (o.first == L"epsx")
217             {
218                 if (o.second->isDouble() == false)
219                 {
220                     Scierror(999, _("%s: Wrong type for input argument #%s: A real scalar expected.\n"), "optim", "epsx");
221                     throw ast::ScilabException();
222                 }
223
224                 pDblEpsx = o.second->getAs<types::Double>();
225                 iEpsx = 0;
226                 pdblEpsx = pDblEpsx->get();
227             }
228         }
229         // get costf
230         opFunctionsManager = new OptimizationFunctions(L"optim");
231         Optimization::addOptimizationFunctions(opFunctionsManager);
232
233
234         if (in[iPos]->isCallable())
235         {
236             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
237             opFunctionsManager->setOptimCostfFunction(pCall);
238             iSim = 2;
239         }
240         else if (in[iPos]->isString())
241         {
242             types::String* pStr = in[iPos]->getAs<types::String>();
243             char* pst = wide_string_to_UTF8(pStr->get(0));
244             bool bOK = opFunctionsManager->setOptimCostfFunction(pStr);
245             iSim = 1;
246
247             if (bOK == false)
248             {
249                 Scierror(50, _("%s: Subroutine not found: %s\n"), "optim", pst);
250                 FREE(pst);
251                 throw ast::ScilabException();
252             }
253
254             memcpy(C2F(optim).nomsub, pst, std::max(size_t(6), strlen(pst)));
255             FREE(pst);
256         }
257         else if (in[iPos]->isList())
258         {
259             types::List* pList = in[iPos]->getAs<types::List>();
260             if (pList->getSize() == 0)
261             {
262                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "optim", iPos + 1, "(string empty)");
263                 throw ast::ScilabException();
264             }
265
266             if (pList->get(0)->isString())
267             {
268                 types::String* pStr = pList->get(0)->getAs<types::String>();
269                 char* pst = wide_string_to_UTF8(pStr->get(0));
270                 bool bOK = opFunctionsManager->setOptimCostfFunction(pStr);
271                 iSim = 1;
272
273                 if (bOK == false)
274                 {
275                     Scierror(50, _("%s: Subroutine not found: %s\n"), "optim", pst);
276                     FREE(pst);
277                     throw ast::ScilabException();
278                 }
279
280                 memcpy(C2F(optim).nomsub, pst, std::max(size_t(6), strlen(pst)));
281                 FREE(pst);
282             }
283             else if (pList->get(0)->isCallable())
284             {
285                 types::Callable* pCall = pList->get(0)->getAs<types::Callable>();
286                 opFunctionsManager->setOptimCostfFunction(pCall);
287                 iSim = 2;
288                 for (int iter = 1; iter < pList->getSize(); iter++)
289                 {
290                     opFunctionsManager->setCostfArgs(pList->get(iter)->getAs<types::InternalType>());
291                 }
292             }
293             else
294             {
295                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "optim", iPos + 1);
296                 throw ast::ScilabException();
297             }
298         }
299         else
300         {
301             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "optim", iPos + 1);
302             throw ast::ScilabException();
303         }
304
305         iPos++;
306
307         // if contr, get binf and bsup
308         if (in[iPos]->isString())
309         {
310             types::String* pStrContr = in[iPos]->getAs<types::String>();
311             if (pStrContr->isScalar() == false || wcscmp(pStrContr->get(0), L"b"))
312             {
313                 Scierror(999, _("%s: Wrong type for input argument #%d: String \"b\" expected.\n"), "optim", iPos + 1);
314                 throw ast::ScilabException();
315             }
316
317             if (in.size() < 5)
318             {
319                 Scierror(77, _("%s: Wrong number of input argument(s): %d or more expected.\n"), "optim", 5);
320                 throw ast::ScilabException();
321             }
322
323             iPos++;
324
325             if (in[iPos]->isDouble())
326             {
327                 pDblBinf = in[iPos]->getAs<types::Double>();
328                 iSizeBinf = pDblBinf->getSize();
329                 pdblBinf = pDblBinf->get();
330             }
331             //        else if(in[iPos]->isPoly())
332             //        {
333             //            pPolyBinf = in[iPos]->getAs<types::Polynom>();
334             //            iSizeBinf = pPolyBinf->getSize();
335             //        }
336             else
337             {
338                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or polynom expected.\n"), "optim", iPos + 1);
339                 throw ast::ScilabException();
340             }
341
342             iPos++;
343
344             if (in[iPos]->isDouble())
345             {
346                 pDblBsub = in[iPos]->getAs<types::Double>();
347                 iSizeBsub = pDblBsub->getSize();
348                 pdblBsub = pDblBsub->get();
349             }
350             //        else if(in[iPos]->isPoly())
351             //        {
352             //            pPolyBsub = in[iPos]->getAs<types::Polynom>();
353             //            iSizeBsub = pPolyBsub->getSize();
354             //        }
355             else
356             {
357                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or polynom expected.\n"), "optim", iPos + 1);
358                 throw ast::ScilabException();
359             }
360
361             iContr = 2;
362             iPos++;
363         }
364
365         // get x0
366
367         if (in[iPos]->isDouble())
368         {
369             pDblX0 = in[iPos]->getAs<types::Double>();
370             iSizeX0 = pDblX0->getSize();
371             if (pDblX0->isComplex())
372             {
373                 iSizeX0 *= 2;
374                 pdblX0 = new double[iSizeX0];
375                 memcpy(pdblX0, pDblX0->get(), pDblX0->getSize() * sizeof(double));
376                 memcpy(pdblX0 + pDblX0->getSize(), pDblX0->getImg(), pDblX0->getSize() * sizeof(double));
377             }
378             else
379             {
380                 pdblX0 = new double[iSizeX0];
381                 memcpy(pdblX0, pDblX0->get(), pDblX0->getSize() * sizeof(double));
382             }
383
384             opFunctionsManager->setXRows(pDblX0->getRows());
385             opFunctionsManager->setXCols(pDblX0->getCols());
386         }
387         //    else if(in[iPos]->isPoly())
388         //    {
389         //
390         //    }
391         else
392         {
393             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or polynom expected.\n"), "optim", iPos + 1);
394             throw ast::ScilabException();
395         }
396
397         if (iContr == 2 && (iSizeX0 != iSizeBinf || iSizeX0 != iSizeBsub))
398         {
399             Scierror(999, _("%s: Bounds and initial guess are incompatible.\n"), "optim");
400             throw ast::ScilabException();
401         }
402
403         if (pDblEpsx != NULL && (pDblEpsx->getSize() != iSizeX0))
404         {
405             Scierror(999, _("%s: Wrong value for input argument #%s: Incorrect stopping parameters.\n"), "optim", "epsx");
406             throw ast::ScilabException();
407         }
408
409         // alloc g output data
410         pdblG = new double[iSizeX0];
411
412         iPos++;
413
414         // get algo
415         if (iPos < in.size() && in[iPos]->isString())
416         {
417             types::String* pStr = in[iPos]->getAs<types::String>();
418             if (pStr->isScalar() == false)
419             {
420                 Scierror(999, _("%s: Wrong type for input argument #%d: Scalar string expected.\n"), "optim", iPos + 1);
421                 throw ast::ScilabException();
422             }
423
424             if (wcscmp(pStr->get(0), L"qn") == 0) // default case
425             {
426                 iAlgo = 1;
427                 iPos++;
428             }
429             else if (wcscmp(pStr->get(0), L"gc") == 0)
430             {
431                 iAlgo = 2;
432                 iPos++;
433             }
434             else if (wcscmp(pStr->get(0), L"nd") == 0)
435             {
436                 iAlgo = 10;
437                 iPos++;
438             }
439             //        else
440             //        {
441             //            Scierror(999, _("%s: Wrong value for input argument #%d: \"qn\", \"gc\", \"nd\" expected.\n"), "optim", iPos + 1);
442             //            throw ast::ScilabException();
443             //        }
444         }
445
446         // get df0 and mem
447         if (iPos < in.size() && in[iPos]->isDouble() && in[iPos]->getAs<types::Double>()->isScalar())
448         {
449             df0 = in[iPos]->getAs<types::Double>()->get(0);
450             iPos++;
451
452             // get mem
453             if (iPos < in.size() && iContr == 1 && (iAlgo == 2 || iAlgo == 10) && in[iPos]->isDouble())
454             {
455                 types::Double* pDbl = in[iPos]->getAs<types::Double>();
456                 if (in[iPos]->isDouble() && pDbl->isScalar() == false)
457                 {
458                     Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "optim", iPos + 1);
459                     throw ast::ScilabException();
460                 }
461
462                 iMem = (int)pDbl->get(0);
463                 bMem = true;
464                 iPos++;
465             }
466         }
467
468         // management of work table
469         if (iAlgo == 1)
470         {
471             // compute size
472
473             if (iContr == 1)
474             {
475                 iWorkSize = (int)(iSizeX0 * ((iSizeX0 + 13) / 2.0));
476                 iWorkSizeI = 0;
477             }
478             else // iContr == 2
479             {
480
481                 iWorkSize = (int)(iSizeX0 * (iSizeX0 + 1) / 2.0 + 4 * iSizeX0 + 1);
482                 iWorkSizeI = 2 * iSizeX0;
483             }
484             /* See bug #9701 for this hard-coded value */
485             /* Fortran underlying algorithm does not support values higher than 46333 */
486             if (iSizeX0 > 46333)
487             {
488                 Scierror(999, _("Can not allocate %.2f MB memory.\n"), (double)(iWorkSize * sizeof(double)) / 1.e6);
489                 throw ast::ScilabException();
490             }
491             try
492             {
493                 // alloc data
494                 pdblWork = new double[iWorkSize];
495
496                 if (iContr == 2)
497                 {
498                     piWork = new int[iWorkSizeI];
499                 }
500             }
501             catch (std::bad_alloc& /*ba*/)
502             {
503                 Scierror(999, _("Can not allocate %.2f MB memory.\n"), (double)(iWorkSize * sizeof(double)) / 1.e6);
504                 throw ast::ScilabException();
505             }
506         }
507
508         // get work
509         if (iPos < in.size() && in[iPos]->isDouble())
510         {
511             if (iAlgo != 1)
512             {
513                 Scierror(137, _("%s: NO hot restart available in this method.\n"), "optim");
514                 throw ast::ScilabException();
515             }
516
517             pDblWork = in[iPos]->getAs<types::Double>();
518
519             if (pDblWork->getSize() != iWorkSize + iWorkSizeI)
520             {
521                 Scierror(137, _("Hot restart: dimension of working table (argument n:%d).\n"), "optim", iPos + 1);
522                 throw ast::ScilabException();
523             }
524
525             double* pdbl = pDblWork->get();
526             if (iContr == 1)
527             {
528                 C2F(dcopy)(&iWorkSize, pdbl, &iOne, pdblWork, &iOne);
529             }
530             else
531             {
532                 C2F(dcopy)(&iWorkSize, pdbl, &iOne, pdblWork, &iOne);
533
534                 for (int i = 0; i < iWorkSizeI; i++)
535                 {
536                     piWork[i] = (int)pdbl[i];
537                 }
538             }
539
540             iMode = 3;
541             iPos++;
542         }
543
544         // get stop
545         while (iPos < in.size() && in[iPos]->isString())
546         {
547             types::String* pStr = in[iPos]->getAs<types::String>();
548             if (wcscmp(pStr->get(0), L"ar") == 0)
549             {
550                 iPos++;
551                 for (int i = iPos; i < in.size(); i++)
552                 {
553                     // get nap, iter, epsg, apsf, epsx
554                     if (in[i]->isDouble())
555                     {
556                         if (pDblNap == NULL)
557                         {
558                             pDblNap = in[i]->getAs<types::Double>();
559                             if (pDblNap->isScalar() == false || pDblNap->isComplex())
560                             {
561                                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "optim", i + 1);
562                                 throw ast::ScilabException();
563                             }
564
565                             iNap = (int)pDblNap->get(0);
566                         }
567                         else if (pDblIter == NULL)
568                         {
569                             pDblIter = in[i]->getAs<types::Double>();
570                             if (pDblIter->isScalar() == false || pDblIter->isComplex())
571                             {
572                                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "optim", i + 1);
573                                 throw ast::ScilabException();
574                             }
575
576                             iItMax = (int)pDblIter->get(0);
577                         }
578                         else if (pDblEpsg == NULL)
579                         {
580                             pDblEpsg = in[i]->getAs<types::Double>();
581                             if (pDblEpsg->isScalar() == false || pDblEpsg->isComplex())
582                             {
583                                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "optim", i + 1);
584                                 throw ast::ScilabException();
585                             }
586
587                             dEpsg = pDblEpsg->get(0);
588                         }
589                         else if (pDblEpsf == NULL)
590                         {
591                             pDblEpsf = in[i]->getAs<types::Double>();
592                             if (pDblEpsf->isScalar() == false || pDblEpsf->isComplex())
593                             {
594                                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "optim", i + 1);
595                                 throw ast::ScilabException();
596                             }
597
598                             dEpsf = pDblEpsf->get(0);
599                         }
600                         else if (pDblEpsx == NULL)
601                         {
602                             pDblEpsx = in[i]->getAs<types::Double>();
603                             if (pDblEpsx->getSize() != iSizeX0)
604                             {
605                                 Scierror(999, _("%s: Wrong value for input argument #%d: Incorrect stopping parameters.\n"), "optim", i + 1);
606                                 throw ast::ScilabException();
607                             }
608
609                             iEpsx = 0;
610                             pdblEpsx = pDblEpsx->get();
611                         }
612                         else
613                         {
614                             Scierror(999, _("%s: Wrong type for input argument #%d: A String expected.\n"), "optim", i + 1);
615                             throw ast::ScilabException();
616                         }
617                     }
618                     else if (in[i]->isString())
619                     {
620                         iPos = i - 1;
621                         break;
622                     }
623                     else
624                     {
625                         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar of a string expected.\n"), "optim", i + 1);
626                         throw ast::ScilabException();
627                     }
628                 }
629             }
630             else if (wcscmp(pStr->get(0), L"in") == 0)
631             {
632                 if (iSim != 1)
633                 {
634                     Scierror(999, _("%s: \"in\" not allowed with simulator defined by a function.\n"), "optim");
635                     throw ast::ScilabException();
636                 }
637
638                 iIndSim = 10;
639                 costf(&iIndSim, &iSizeX0, pDblX0->get(), &dF, pdblG, NULL, NULL, NULL);
640
641                 if (iIndSim == 0)
642                 {
643                     Scierror(131, _("%s: Stop requested by simulator (ind=0).\n"), "optim");
644                     throw ast::ScilabException();
645                 }
646                 else if (iIndSim < 0)
647                 {
648                     Scierror(134, _("%s: Problem with initial constants in simul.\n"), "optim");
649                     throw ast::ScilabException();
650                 }
651
652                 if (piIzs)
653                 {
654                     delete[] piIzs;
655                 }
656                 if (pfRzs)
657                 {
658                     delete[] pfRzs;
659                 }
660                 if (pdblDzs && iDzs)
661                 {
662                     delete[] pdblDzs;
663                 }
664
665                 piIzs   = new int[C2F(nird).nizs];
666                 pfRzs   = new float[C2F(nird).nrzs];
667                 pdblDzs = new double[C2F(nird).ndzs];
668
669                 iIndSim = 11;
670                 costf(&iIndSim, &iSizeX0, pDblX0->get(), &dF, pdblG, piIzs, pfRzs, pdblDzs);
671
672                 if (iIndSim == 0)
673                 {
674                     Scierror(131, _("%s: Stop requested by simulator (ind=0).\n"), "optim");
675                     throw ast::ScilabException();
676                 }
677                 else if (iIndSim < 0)
678                 {
679                     Scierror(134, _("%s: Problem with initial constants in simul.\n"), "optim");
680                     throw ast::ScilabException();
681                 }
682             }
683             else if (wcscmp(pStr->get(0), L"ti") == 0)
684             {
685                 iPos++;
686                 if (in[iPos]->isDouble() == false)
687                 {
688                     Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "optim", iPos + 1);
689                     throw ast::ScilabException();
690                 }
691
692                 pDblTi = in[iPos]->getAs<types::Double>();
693                 C2F(nird).nizs = pDblTi->getSize();
694                 if (piIzs)
695                 {
696                     delete[] piIzs;
697                 }
698                 piIzs = new int[pDblTi->getSize()];
699
700                 for (int i = 0; i < pDblTi->getSize(); i++)
701                 {
702                     piIzs[i] = (int)pDblTi->get(i);
703                 }
704             }
705             else if (wcscmp(pStr->get(0), L"td") == 0)
706             {
707                 iPos++;
708                 if (in[iPos]->isDouble() == false)
709                 {
710                     Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "optim", iPos + 1);
711                     throw ast::ScilabException();
712                 }
713
714                 if (pdblDzs && iDzs)
715                 {
716                     delete[] pdblDzs;
717                 }
718
719                 pDblTd = in[iPos]->getAs<types::Double>();
720                 C2F(nird).ndzs = pDblTd->getSize();
721                 pdblDzs = pDblTd->get();
722                 iDzs = 0;
723             }
724             else if (wcscmp(pStr->get(0), L"si") == 0)
725             {
726                 iSaveI = 1;
727             }
728             else if (wcscmp(pStr->get(0), L"sd") == 0)
729             {
730                 iSaveD = 1;
731             }
732             else
733             {
734                 Scierror(999, _("%s: Wrong value for input argument #%d: \"ar\", \"in\", \"ti\" or \"td\" not allowed.\n"), "optim", iPos + 1);
735                 throw ast::ScilabException();
736             }
737
738             iPos++;
739         }
740
741         // initialisation eventuelle de f et g
742         if (iNap < 2 || iItMax < 1)
743         {
744             iArret = 1;
745         }
746
747         if (iContr == 1 && (iAlgo == 2 || iAlgo == 10) ||
748                 (iContr == 2 && iAlgo == 1 && pdblWork) ||
749                 (iArret == 1))
750         {
751             iIndSim = 4;
752             costf(&iIndSim, &iSizeX0, pDblX0->get(), &dF, pdblG, piIzs, pfRzs, pdblDzs);
753
754             if (iIndSim == 0)
755             {
756                 Scierror(131, _("%s: Stop requested by simulator (ind=0).\n"), "optim");
757                 throw ast::ScilabException();
758             }
759             else if (iIndSim < 0)
760             {
761                 Scierror(134, _("%s: Problem with initial constants in simul.\n"), "optim");
762                 throw ast::ScilabException();
763             }
764
765             if (iNap < 2 || iItMax < 1)
766             {
767                 // skip perform operation part
768                 iContr = 3;
769             }
770         }
771
772         /*** perform operations ***/
773         // n1qn1 : Quasi-Newton without constraints
774         if (iContr == 1 && iAlgo == 1) // bounds not setted && algo qn
775         {
776             pdblVar = new double[iSizeX0];
777             for (int i = 0; i < iSizeX0; i++)
778             {
779                 pdblVar[i] = 0.10;
780             }
781
782             int iItmax1   = iItMax;
783             int iNap1     = iNap;
784             double dEpsg1 = dEpsg;
785
786             C2F(n1qn1)(costf, &iSizeX0, pdblX0, &dF, pdblG,
787                        pdblVar, &dEpsg, &iMode, &iItMax, &iNap, &iPrint, &io, pdblWork,
788                        piIzs, pfRzs, pdblDzs);
789
790             dEpsg = sqrt(dEpsg);
791
792             iIndOpt = 9;
793             if (dEpsg1 >= dEpsg)
794             {
795                 iIndOpt = 1;
796             }
797             else if (iNap >= iNap1)
798             {
799                 iIndOpt = 4;
800             }
801             else if (iItMax >= iItmax1)
802             {
803                 iIndOpt = 5;
804             }
805
806             if (checkOptimError(iArret, iIndOpt, iPrint, dEpsg))
807             {
808                 throw ast::ScilabException();
809             }
810         }
811         // algorithme n1qn3 : Gradient Conjugate without constraints
812         else if (iContr == 1 && iAlgo == 2) // bounds not setted && algo gc
813         {
814             double dxmin = dEpsg;
815             double dZng = 0;
816
817             if (bMem == false)
818             {
819                 iMem = 10;
820             }
821
822             // compute epsrel
823             for (int i = 0; i < iSizeX0; i++)
824             {
825                 dZng += (pdblG[i] * pdblG[i]);
826             }
827
828             dZng = sqrt(dZng);
829
830             if (dZng > 0)
831             {
832                 dEpsg /= dZng;
833             }
834
835             // compute dxmin
836             if (iEpsx == 0)
837             {
838                 dxmin = pdblEpsx[0];
839                 if (iSizeX0 > 1)
840                 {
841                     for (int i = 0; i < iSizeX0; i++)
842                     {
843                         dxmin = std::min(dxmin, pdblEpsx[i]);
844                     }
845                 }
846             }
847
848             // work table
849             iWorkSize = 4 * iSizeX0 + iMem * (2 * iSizeX0 + 1);
850             pdblWork = new double[iWorkSize];
851
852             iIndSim = 4;
853             costf(&iIndSim, &iSizeX0, pDblX0->get(), &dF, pdblG, piIzs, pfRzs, pdblDzs);
854             C2F(n1qn3)( costf, C2F(fuclid), C2F(ctonb), C2F(ctcab), &iSizeX0, pdblX0,
855                         &dF, pdblG, &dxmin, &df0, &dEpsg, &iPrint, &io, &iMode, &iItMax,
856                         &iNap, pdblWork, &iWorkSize, piIzs, pfRzs, pdblDzs);
857
858             switch (iMode)
859             {
860                 case 0 :
861                     iIndOpt = 0;
862                     break;
863                 case 1 :
864                     iIndOpt = 1;
865                     break;
866                 case 2 :
867                     iIndOpt = -10;
868                     break;
869                 case 7 :
870                 case 3 :
871                     iIndOpt = 7;
872                     break;
873                 case 4 :
874                     iIndOpt = 5;
875                     break;
876                 case 5 :
877                     iIndOpt = 4;
878                     break;
879                 case 6 :
880                     iIndOpt = 3;
881                     break;
882                 default :
883                     iIndOpt = 9;
884             }
885
886             if (checkOptimError(iArret, iIndOpt, iPrint, dEpsg))
887             {
888                 throw ast::ScilabException();
889             }
890
891         }
892         // optimiseur n1fc1 : non smooth without constraints
893         else if (iContr == 1 && iAlgo == 10) // bounds not setted && algo nd
894         {
895             if (bMem == false)
896             {
897                 iMem = 10;
898             }
899
900             int iNitv   = 2 * iMem + 2;
901             int iNtv1   = 5 * iSizeX0 + (iSizeX0 + 4) * iMem;
902             int iNtv2   = (iMem + 9) * iMem + 8;
903
904             piWork      = new int[iNitv];
905             pdblWork    = new double[iNtv1];
906             pdblWork2   = new double[iNtv2];
907
908             if (iEpsx == 1)
909             {
910                 pdblEpsx = new double[iSizeX0];
911                 C2F(dcopy)(&iSizeX0, &dTol, &iZero, pdblEpsx, &iOne);
912             }
913
914             C2F(n1fc1)(costf, C2F(fuclid), &iSizeX0, pdblX0, &dF, pdblG,
915                        pdblEpsx, &df0, &dEpsf, &dTol, &iPrint, &io, &iMode,
916                        &iItMax, &iNap, &iMem, piWork, pdblWork, pdblWork2,
917                        piIzs, pfRzs, pdblDzs);
918
919             switch (iMode)
920             {
921                 case 0 :
922                     iIndOpt = 0;
923                     break;
924                 case 1 :
925                     iIndOpt = 2;
926                     break;
927                 case 2 :
928                     iIndOpt = -10;
929                     break;
930                 case 4 :
931                     iIndOpt = 5;
932                     break;
933                 case 5 :
934                     iIndOpt = 4;
935                     break;
936                 case 6 :
937                     iIndOpt = 3;
938                     break;
939                 default :
940                     iIndOpt = 9;
941             }
942
943             if (checkOptimError(iArret, iIndOpt, iPrint, dEpsg))
944             {
945                 throw ast::ScilabException();
946             }
947         }
948         // optimiseur qnbd : quasi-newton with bound constraints
949         else if (iContr == 2 && iAlgo == 1) // bounds setted && algo qn
950         {
951             int iNfac = 0;
952             if (iEpsx == 1)
953             {
954                 pdblEpsx = new double[iSizeX0];
955                 C2F(dcopy)(&iSizeX0, &dTol, &iZero, pdblEpsx, &iOne);
956             }
957
958             iIndOpt = 1 + pDblWork ? 1 : 0;
959             C2F(qnbd)(&iIndOpt, costf, &iSizeX0, pdblX0, &dF, pdblG, &iPrint, &io, &dTol, &iNap,
960                       &iItMax, &dEpsf, &dEpsg, pdblEpsx, &df0, pdblBinf, pdblBsub,
961                       &iNfac, pdblWork, &iWorkSize, piWork, &iWorkSizeI, piIzs, pfRzs, pdblDzs);
962
963             if (checkOptimError(iArret, iIndOpt, iPrint, dEpsg))
964             {
965                 throw ast::ScilabException();
966             }
967         }
968         // optimiseur gcbd : Gradient Conjugate with bound constraints
969         else if (iContr == 2 && iAlgo == 2) // bounds setted && algo gc
970         {
971             int iNfac = 0;
972             int nt = 2;
973             if (bMem)
974             {
975                 nt = std::max(1, iMem / 3);
976             }
977
978             iWorkSize   = iSizeX0 * (5 + 3 * nt) + 2 * nt + 1;
979             iWorkSizeI  = iSizeX0 + nt + 1;
980             pdblWork    = new double[iWorkSize];
981             piWork      = new int[iWorkSizeI];
982
983             if (iEpsx == 1)
984             {
985                 pdblEpsx = new double[iSizeX0];
986                 C2F(dcopy)(&iSizeX0, &dTol, &iZero, pdblEpsx, &iOne);
987             }
988
989             iIndOpt = 1;
990             C2F(gcbd)(&iIndOpt, costf, C2F(optim).nomsub, &iSizeX0, pdblX0, &dF, pdblG,
991                       &iPrint, &io, &dTol, &iNap, &iItMax, &dEpsf, &dEpsg, pdblEpsx, &df0,
992                       pdblBinf, pdblBsub, &iNfac, pdblWork, &iWorkSize, piWork, &iWorkSizeI,
993                       piIzs, pfRzs, pdblDzs);
994
995             if (checkOptimError(iArret, iIndOpt, iPrint, dEpsg))
996             {
997                 throw ast::ScilabException();
998             }
999         }
1000         else if (iContr != 3) // bad algo
1001         {
1002             Scierror(136, _("%s: This method is NOT implemented.\n"), "optim");
1003             throw ast::ScilabException();
1004         }
1005
1006         /*** return output arguments ***/
1007         int iRetCount1 = _iRetCount - iSaveI - iSaveD;
1008         if (iRetCount1 == 0)
1009         {
1010             Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "optim", iSaveI + iSaveD, iSaveI + iSaveD + 1);
1011             throw ast::ScilabException();
1012         }
1013
1014         // return f
1015         out.push_back(new types::Double(dF));
1016
1017         // return x
1018         if (iRetCount1 >= 2)
1019         {
1020             if (pDblX0)
1021             {
1022                 types::Double* pDbl = new types::Double(pDblX0->getDims(), pDblX0->getDimsArray(), pDblX0->isComplex());
1023                 double* pdblReal = pDbl->get();
1024                 memcpy(pdblReal, pdblX0, pDbl->getSize() * sizeof(double));
1025
1026                 if (pDbl->isComplex())
1027                 {
1028                     double* pdblImg = pDbl->getImg();
1029                     memcpy(pdblImg, pdblX0 + pDbl->getSize(), pDbl->getSize() * sizeof(double));
1030                 }
1031
1032                 out.push_back(pDbl);
1033             }
1034             //        else // if (pPolyX0)
1035             //        {
1036             //
1037             //        }
1038         }
1039
1040         // return g
1041         if (iRetCount1 >= 3)
1042         {
1043             if (pdblG)
1044             {
1045                 types::Double* pDbl = new types::Double(pDblX0->getDims(), pDblX0->getDimsArray(), pDblX0->isComplex());
1046                 double* pdblReal = pDbl->get();
1047                 memcpy(pdblReal, pdblG, pDbl->getSize() * sizeof(double));
1048
1049                 if (pDbl->isComplex())
1050                 {
1051                     double* pdblImg = pDbl->getImg();
1052                     memcpy(pdblImg, pdblG + pDbl->getSize(), pDbl->getSize() * sizeof(double));
1053                 }
1054
1055                 out.push_back(pDbl);
1056             }
1057             //        else // if (pPolyX0)
1058             //        {
1059             //
1060             //        }
1061         }
1062
1063         // return work
1064         if (iRetCount1 >= 4)
1065         {
1066             if (iAlgo != 1)
1067             {
1068                 Scierror(137, _("%s: NO hot restart available in this method.\n"), "optim");
1069                 throw ast::ScilabException();
1070             }
1071
1072             if (iContr == 1)
1073             {
1074                 types::Double* pDbl = new types::Double(1, iWorkSize);
1075                 double* pdbl = pDbl->get();
1076                 C2F(dcopy)(&iWorkSize, pdblWork, &iOne, pdbl, &iOne);
1077                 out.push_back(pDbl);
1078             }
1079             else if (iContr == 2)
1080             {
1081                 types::Double* pDbl = new types::Double(1, iWorkSize + iWorkSizeI);
1082                 double* pdbl = pDbl->get();
1083                 C2F(dcopy)(&iWorkSize, pdblWork, &iOne, pdbl, &iOne);
1084                 for (int i = 0; i < iWorkSizeI; i++)
1085                 {
1086                     pdbl[iWorkSize + i] = (double)(piWork[i]);
1087                 }
1088
1089                 out.push_back(pDbl);
1090             }
1091             else // iContr == 3
1092             {
1093                 out.push_back(pDblWork->clone());
1094             }
1095         }
1096
1097         if (iRetCount1 >= 5)
1098         {
1099             out.push_back(new types::Double((double)iItMax));
1100         }
1101
1102         if (iRetCount1 >= 6)
1103         {
1104             out.push_back(new types::Double((double)iNap));
1105         }
1106
1107         if (iRetCount1 >= 7)
1108         {
1109             out.push_back(new types::Double((double)iIndOpt));
1110         }
1111
1112         if (iSaveI)
1113         {
1114             if (C2F(nird).nizs == 0)
1115             {
1116                 out.push_back(types::Double::Empty());
1117             }
1118             else
1119             {
1120                 types::Double* pDbl = new types::Double(1, C2F(nird).nizs);
1121                 double* pdbl = pDbl->get();
1122                 for (int i = 0; i < C2F(nird).nizs; i++)
1123                 {
1124                     pdbl[i] = (double)piIzs[i];
1125                 }
1126
1127                 out.push_back(pDbl);
1128             }
1129         }
1130
1131         if (iSaveD)
1132         {
1133             if (pdblDzs == NULL || C2F(nird).ndzs == 0)
1134             {
1135                 out.push_back(types::Double::Empty());
1136             }
1137             else
1138             {
1139                 types::Double* pDbl = new types::Double(1, C2F(nird).ndzs);
1140                 memcpy(pDbl->get(), pdblDzs, C2F(nird).ndzs * sizeof(double));
1141                 out.push_back(pDbl);
1142             }
1143         }
1144
1145         ret = types::Function::OK;
1146     }
1147     catch (const ast::InternalError& e)
1148     {
1149         char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1150         Scierror(999, pstrMsg);
1151         FREE(pstrMsg);
1152     }
1153     catch (const ast::ScilabException& /* e */)
1154     {
1155         // free memory, then return error
1156     }
1157
1158     /*** free memory ***/
1159     if (opFunctionsManager)
1160     {
1161         Optimization::removeOptimizationFunctions();
1162         delete opFunctionsManager;
1163     }
1164
1165     if (piIzs)
1166     {
1167         delete[] piIzs;
1168     }
1169
1170     if (pfRzs)
1171     {
1172         delete[] pfRzs;
1173     }
1174
1175     if (pdblG)
1176     {
1177         delete[] pdblG;
1178     }
1179
1180     if (pdblDzs && iDzs)
1181     {
1182         delete[] pdblDzs;
1183     }
1184
1185     if (pdblWork)
1186     {
1187         delete[] pdblWork;
1188     }
1189
1190     if (pdblWork2)
1191     {
1192         delete[] pdblWork2;
1193     }
1194
1195     if (piWork)
1196     {
1197         delete[] piWork;
1198     }
1199
1200     if (pdblX0)
1201     {
1202         delete[] pdblX0;
1203     }
1204
1205     if (pdblVar)
1206     {
1207         delete[] pdblVar;
1208     }
1209
1210     if (pdblEpsx && iEpsx)
1211     {
1212         delete[] pdblEpsx;
1213     }
1214
1215     return ret;
1216 }
1217