3cbf0b6425a85890dd76100d55c692b72b0d0251
[scilab.git] / scilab / modules / optimization / sci_gateway / cpp / sci_lsqrsolve.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
17 #include "optimization_gw.hxx"
18 #include "function.hxx"
19 #include "double.hxx"
20 #include "string.hxx"
21 #include "list.hxx"
22 #include "optimizationfunctions.hxx"
23
24 extern "C"
25 {
26 #include "localization.h"
27 #include "Scierror.h"
28 #include "sciprint.h"
29 #include "scioptimfunctions.h"
30 #include "sci_malloc.h"
31 }
32 /*--------------------------------------------------------------------------*/
33
34 types::Function::ReturnValue sci_lsqrsolve(types::typed_list &in, int _iRetCount, types::typed_list &out)
35 {
36     OptimizationFunctions* opFunctionsManager = NULL;
37
38     types::Double* pDblX    = NULL;
39     types::Double* pDblM    = NULL;
40     types::Double* pDblV    = NULL;
41
42     double* pdblDiag = NULL;
43
44     int iSizeX      = 0;
45     int iWorkSize   = 0;
46     int iInfo       = 0;
47     int iM          = 0;
48     int iPos        = 0;
49     int iNprint     = 0;
50     int iMode       = 1;
51     int iNfev       = 0;
52     int iNjev       = 0;
53
54     bool bJac       = false;
55
56     // stop
57     double dFtol    = 1.0e-8;
58     double dXtol    = 1.0e-8;
59     double dGtol    = 1.0e-5;
60     int iMaxfev     = 1000;
61     double dEpsfcn  = 0;
62     double dFactor  = 100;
63
64     if (in.size() < 3 || in.size() > 6)
65     {
66         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "lsqrsolve", 3, 6);
67         return types::Function::Error;
68     }
69
70     if (_iRetCount > 3)
71     {
72         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lsqrsolve", 1, 3);
73         return types::Function::Error;
74     }
75
76     /*** get inputs arguments ***/
77     // get X
78     if (in[0]->isDouble() == false)
79     {
80         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "lsqrsolve", 1);
81         return types::Function::Error;
82     }
83
84     pDblX = in[0]->clone()->getAs<types::Double>();
85
86     if (pDblX->isComplex())
87     {
88         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "lsqrsolve", 1);
89         return types::Function::Error;
90     }
91
92     iSizeX = pDblX->getSize();
93
94     // get function
95     opFunctionsManager = new OptimizationFunctions(L"lsqrsolve");
96     Optimization::addOptimizationFunctions(opFunctionsManager);
97     opFunctionsManager->setXRows(pDblX->getRows());
98     opFunctionsManager->setXCols(pDblX->getCols());
99
100     if (in[1]->isCallable())
101     {
102         types::Callable* pCall = in[1]->getAs<types::Callable>();
103         opFunctionsManager->setFsolveFctFunction(pCall);
104     }
105     else if (in[1]->isString())
106     {
107         types::String* pStr = in[1]->getAs<types::String>();
108         char* pst = wide_string_to_UTF8(pStr->get(0));
109         bool bOK = opFunctionsManager->setFsolveFctFunction(pStr);
110
111         if (bOK == false)
112         {
113             Scierror(50, _("%s: Subroutine not found: %s\n"), "lsqrsolve", pst);
114             FREE(pst);
115             return types::Function::Error;
116         }
117
118         FREE(pst);
119     }
120     else if (in[1]->isList())
121     {
122         types::List* pList = in[1]->getAs<types::List>();
123         if (pList->getSize() == 0)
124         {
125             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "lsqrsolve", 2, "(string empty)");
126             return types::Function::Error;
127         }
128
129         if (pList->get(0)->isString())
130         {
131             types::String* pStr = pList->get(0)->getAs<types::String>();
132             char* pst = wide_string_to_UTF8(pStr->get(0));
133             bool bOK = opFunctionsManager->setFsolveFctFunction(pStr);
134
135             if (bOK == false)
136             {
137                 Scierror(50, _("%s: Subroutine not found: %s\n"), "lsqrsolve", pst);
138                 FREE(pst);
139                 return types::Function::Error;
140             }
141
142             FREE(pst);
143         }
144         else if (pList->get(0)->isCallable())
145         {
146             types::Callable* pCall = pList->get(0)->getAs<types::Callable>();
147             opFunctionsManager->setFsolveFctFunction(pCall);
148             for (int iter = 1; iter < pList->getSize(); iter++)
149             {
150                 opFunctionsManager->setFsolveFctArgs(pList->get(iter)->getAs<types::InternalType>());
151             }
152         }
153         else
154         {
155             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "lsqrsolve", 2);
156             return types::Function::Error;
157         }
158     }
159     else
160     {
161         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "lsqrsolve", 2);
162         return types::Function::Error;
163     }
164
165     // get M
166     if (in[2]->isDouble() == false)
167     {
168         Scierror(999, _("%s: Wrong type for input argument #%d: A real Scalar expected.\n"), "lsqrsolve", 3);
169         return types::Function::Error;
170     }
171
172     pDblM = in[2]->getAs<types::Double>();
173
174     if (pDblM->isComplex() || pDblM->isScalar() == false)
175     {
176         Scierror(999, _("%s: Wrong type for input argument #%d: A real Scalar expected.\n"), "lsqrsolve", 3);
177         return types::Function::Error;
178     }
179
180     iM = (int)pDblM->get(0);
181
182     if (iM < iSizeX)
183     {
184         Scierror(999, _("%s: Wrong value for input argument #%d: At most %d expected.\n"), "lsqrsolve", 3, iM);
185         return types::Function::Error;
186     }
187
188     iPos = 3;
189
190     // get jacobian
191     if (in.size() >= 4)
192     {
193         if (in[iPos]->isCallable())
194         {
195             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
196             opFunctionsManager->setFsolveJacFunction(pCall);
197
198             bJac = true;
199             iPos++;
200         }
201         else if (in[iPos]->isString())
202         {
203             types::String* pStr = in[iPos]->getAs<types::String>();
204             char* pst = wide_string_to_UTF8(pStr->get(0));
205             bool bOK = opFunctionsManager->setFsolveJacFunction(pStr);
206
207             if (bOK == false)
208             {
209                 Scierror(50, _("%s: Subroutine not found: %s\n"), "lsqrsolve", pst);
210                 FREE(pst);
211                 return types::Function::Error;
212             }
213
214             bJac = true;
215             iPos++;
216             FREE(pst);
217         }
218         else if (in[iPos]->isList())
219         {
220             types::List* pList = in[iPos]->getAs<types::List>();
221             if (pList->getSize() == 0)
222             {
223                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "lsqrsolve", 4, "(string empty)");
224                 return types::Function::Error;
225             }
226
227             if (pList->get(0)->isString())
228             {
229                 types::String* pStr = pList->get(0)->getAs<types::String>();
230                 char* pst = wide_string_to_UTF8(pStr->get(0));
231                 bool bOK = opFunctionsManager->setFsolveJacFunction(pStr);
232
233                 if (bOK == false)
234                 {
235                     Scierror(50, _("%s: Subroutine not found: %s\n"), "lsqrsolve", pst);
236                     return types::Function::Error;
237                 }
238
239                 bJac = true;
240                 iPos++;
241                 FREE(pst);
242             }
243             else if (pList->get(0)->isCallable())
244             {
245                 types::Callable* pCall = pList->get(0)->getAs<types::Callable>();
246                 opFunctionsManager->setFsolveJacFunction(pCall);
247                 for (int iter = 1; iter < pList->getSize(); iter++)
248                 {
249                     opFunctionsManager->setFsolveJacArgs(pList->get(iter)->getAs<types::InternalType>());
250                 }
251
252                 bJac = true;
253                 iPos++;
254             }
255             else
256             {
257                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "lsqrsolve", 4);
258                 return types::Function::Error;
259             }
260         }
261     }
262
263     if (bJac == false && in.size() > 5)
264     {
265         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "lsqrsolve", 3, 5);
266         return types::Function::Error;
267     }
268
269     // get stop
270     if (iPos < in.size())
271     {
272         if (in[iPos]->isDouble())
273         {
274             types::Double* pDblStop = in[iPos]->getAs<types::Double>();
275             if (pDblStop->getRows() != 1 || pDblStop->getCols() != 6)
276             {
277                 Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d x %d expected.\n"), "lsqrsolve", iPos + 1, 1, 6);
278                 return types::Function::Error;
279             }
280
281             dFtol   = pDblStop->get(0);
282             dXtol   = pDblStop->get(1);
283             dGtol   = pDblStop->get(2);
284             iMaxfev = static_cast<int>(pDblStop->get(3));
285             dEpsfcn = pDblStop->get(4);
286             dFactor = pDblStop->get(5);
287
288             if (dFtol < 0 || dXtol < 0 || dGtol < 0 || iMaxfev <= 0 || dFactor <= 0)
289             {
290                 Scierror(999, _("%s: Wrong value for input argument #%d.\n"), "lsqrsolve", iPos + 1, 1, 6);
291                 return types::Function::Error;
292             }
293         }
294         else
295         {
296             Scierror(999, _("%s: Wrong type for input argument #%d: A vector expected.\n"), "lsqrsolve", iPos + 1);
297             return types::Function::Error;
298         }
299
300         iPos++;
301     }
302
303     if (iPos < in.size())
304     {
305         if (in[iPos]->isDouble())
306         {
307             types::Double* pDblDiag = in[iPos]->getAs<types::Double>();
308             if (pDblDiag->getSize() != iSizeX)
309             {
310                 Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d expected.\n"), "lsqrsolve", iPos + 1, iSizeX);
311                 return types::Function::Error;
312             }
313
314             pdblDiag = pDblDiag->get();
315
316             for (int i = 0; i < iSizeX; i++)
317             {
318                 if (pdblDiag[i] <= 0)
319                 {
320                     Scierror(999, _("%s: Wrong value for input argument #%d: A strictly positive vector expected.\n"), "lsqrsolve", iPos + 1);
321                     return types::Function::Error;
322                 }
323             }
324
325             iMode = 2;
326         }
327         else
328         {
329             Scierror(999, _("%s: Wrong type for input argument #%d: A vector expected.\n"), "lsqrsolve", iPos + 1);
330             return types::Function::Error;
331         }
332     }
333
334     /*** perform operations ***/
335     // alloc working table
336     double* pdblJac     = new double[iSizeX * iM];
337     int* piPvt          = new int[iSizeX];
338     double* pDblQtf     = new double[iSizeX];
339     double* pdblWork1   = new double[iSizeX];
340     double* pdblWork2   = new double[iSizeX];
341     double* pdblWork3   = new double[iSizeX];
342     double* pdblWork4   = new double[iM];
343
344     if (pdblDiag == NULL) // iMode != 2
345     {
346         pdblDiag = new double[iSizeX];
347     }
348
349     // alloc output data
350     pDblV = new types::Double(iM, 1);
351
352     char const* pstrFunc = "fct";
353     try
354     {
355         if (bJac)
356         {
357             pstrFunc = "jac";
358             C2F(lmder)( lsqrjac, &iM, &iSizeX, pDblX->get(), pDblV->get(), pdblJac, &iM, &dFtol,
359                         &dXtol, &dGtol, &iMaxfev, pdblDiag, &iMode, &dFactor, &iNprint, &iInfo,
360                         &iNfev, &iNjev, piPvt, pDblQtf, pdblWork1, pdblWork2, pdblWork3, pdblWork4);
361         }
362         else
363         {
364             C2F(lmdif)( lsqrfct, &iM, &iSizeX, pDblX->get(), pDblV->get(), &dFtol, &dXtol, &dGtol,
365                         &iMaxfev, &dEpsfcn, pdblDiag, &iMode, &dFactor, &iNprint, &iInfo, &iNfev,
366                         pdblJac, &iM, piPvt, pDblQtf, pdblWork1, pdblWork2, pdblWork3, pdblWork4);
367         }
368     }
369     catch (const ast::InternalError &e)
370     {
371         char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
372         sciprint(_("%s: exception caught in '%s' subroutine.\n"), "lsqrsolve", pstrFunc);
373         Scierror(999, pstrMsg);
374         FREE(pstrMsg);
375         delete pDblX;
376         delete piPvt;
377         delete pDblQtf;
378         delete pdblWork1;
379         delete pdblWork2;
380         delete pdblWork3;
381         delete pdblWork4;
382         delete pdblJac;
383         if (iMode != 2)
384         {
385             delete pdblDiag;
386         }
387
388         return types::Function::Error;
389     }
390
391     delete piPvt;
392     delete pDblQtf;
393     delete pdblWork1;
394     delete pdblWork2;
395     delete pdblWork3;
396     delete pdblWork4;
397     delete pdblJac;
398     if (iMode != 2)
399     {
400         delete pdblDiag;
401     }
402
403     /*** return output arguments ***/
404     out.push_back(pDblX);
405
406     if (_iRetCount >= 2)
407     {
408         out.push_back(pDblV);
409     }
410     else
411     {
412         delete pDblV;
413     }
414
415     // info = 0   improper input parameters.
416     // info = 1   algorithm estimates that the relative error
417     // between x and the solution is at most tol.
418     // info = 2   number of calls to fcn with iflag = 1 has
419     // reached 100*(n+1).
420     // info = 3   tol is too small. no further improvement in
421     // the approximate solution x is possible.
422     // info = 4   iteration is not making good progress.
423     if (_iRetCount == 3) // return
424     {
425         out.push_back(new types::Double((double)iInfo));
426     }
427     else // printf warning
428     {
429         switch (iInfo)
430         {
431             case 5:
432                 sciprint(_("%s: %s: Number of calls to %s has reached or exceeded %s.\n"), _("Warning"), "lsqrsolve", "fct", "maxfev");
433                 break;
434             case 6:
435                 sciprint(_("%s: %s: %s is too small. No further reduction in the criterion is possible.\n"), _("Warning"), "lsqrsolve", "ftol");
436                 break;
437             case 7:
438                 sciprint(_("%s: %s: %s is too small. No further reduction in the criterion is possible.\n"), _("Warning"), "lsqrsolve", "xtol");
439                 break;
440             case 8:
441                 sciprint(_("%s: %s: %s is too small. %s is orthogonal to the columns of the Jacobian to machine precision.\n"), _("Warning"), "lsqrsolve", "gtol", "fvec");
442                 break;
443         }
444     }
445
446     return types::Function::OK;
447 }
448