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