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