rollback %ls to %s, it was better before
[scilab.git] / scilab / modules / randlib / sci_gateway / cpp / sci_grand.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - 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 "randlib_gw.hxx"
14 #include "function.hxx"
15 #include "double.hxx"
16 #include "string.hxx"
17 #include "configvariable.hxx"
18
19 extern "C"
20 {
21 #include "MALLOC.h"
22 #include "localization.h"
23 #include "Scierror.h"
24 #include "sciprint.h"
25 #include "grand.h"
26 #include "clcg4.h"
27 #include "others_generators.h"
28 }
29 /*--------------------------------------------------------------------------*/
30
31 types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     enum {MT, KISS, CLCG4, CLCG2, URAND, FSULTRA};
34
35     //  names at the scilab level
36     wchar_t* names_gen[6] = {L"mt", L"kiss", L"clcg4", L"clcg2", L"urand", L"fsultra"};
37
38     types::String* pStrMethod   = NULL;
39     types::String* pStrGenOrPhr = NULL;
40
41     std::vector<types::Double*> vectpDblInput;
42
43     int iStrPos     = 0;
44     int iPos        = 0;
45     int meth        = -1;
46     int iRows       = -1;
47     int iCols       = -1;
48     int iNumIter    = -1;
49
50     int current_base_gen = ConfigVariable::getCurrentBaseGen();
51
52 // *** check the maximal number of input args. ***
53     if(in.size() > 6)
54     {
55         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "grand", 1, 6);
56         return types::Function::Error;
57     }
58
59 // *** check number of output args. ***
60     if(_iRetCount > 1)
61     {
62         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "grand", 1);
63         return types::Function::Error;
64     }
65
66 // *** find the mothod string. ***
67     for(int i = 0; i < in.size(); i++)
68     {
69         if(in[i]->isString())
70         {
71             pStrMethod = in[i]->getAs<types::String>();
72             iStrPos = i;
73             break;
74         }
75     }
76
77     if(pStrMethod == NULL)
78     {
79         Scierror(78, _("%s: Wrong number of output argument(s): At least %d string expected.\n"), "grand", 1);
80         return types::Function::Error;
81     }
82
83     wchar_t* wcsMeth = pStrMethod->get(0);
84     int iNumInputArg = 5;
85     if(wcscmp(wcsMeth, L"bet") == 0) // beta
86     {
87         meth = 1;
88     }
89     else if(wcscmp(wcsMeth,L"bin") == 0) // binomial
90     {
91         meth = 2;
92     }
93     else if(wcscmp(wcsMeth, L"nbn") == 0) // negative binomial
94     {
95         meth = 3;
96     }
97     else if(wcscmp(wcsMeth, L"chi") == 0) // chisquare
98     {
99         iNumInputArg = 4;
100         meth = 4;
101     }
102     else if(wcscmp(wcsMeth, L"nch") == 0) // non central chisquare
103     {
104         meth = 5;
105     }
106     else if(wcscmp(wcsMeth, L"exp") == 0) // exponential
107     {
108         iNumInputArg = 4;
109         meth = 6;
110     }
111     else if(wcscmp(wcsMeth, L"f") == 0) // F variance ratio
112     {
113         meth = 7;
114     }
115     else if(wcscmp(wcsMeth, L"nf") == 0) // non central F variance ratio
116     {
117         iNumInputArg = 6;
118         meth = 8;
119     }
120     else if(wcscmp(wcsMeth, L"gam") == 0) // gamma
121     {
122         meth = 9;
123     }
124     else if(wcscmp(wcsMeth, L"nor") == 0) // Gauss Laplace (normal)
125     {
126         meth = 10;
127     }
128     else if(wcscmp(wcsMeth, L"mn") == 0) // multivariate gaussian (multivariate normal)
129     {
130         iNumInputArg = 4;
131         meth = 11;
132     }
133     else if(wcscmp(wcsMeth, L"geom") == 0) // geometric
134     {
135         iNumInputArg = 4;
136         meth = 12;
137     }
138     else if(wcscmp(wcsMeth, L"markov") == 0) // markov
139     {
140         iNumInputArg = 4;
141         meth = 13;
142     }
143     else if(wcscmp(wcsMeth, L"mul") == 0) // multinomial
144     {
145         iNumInputArg = 4;
146         meth = 14;
147     }
148     else if(wcscmp(wcsMeth, L"poi") == 0) // Poisson
149     {
150         iNumInputArg = 4;
151         meth = 15;
152     }
153     else if(wcscmp(wcsMeth, L"prm") == 0) // random permutations
154     {
155         iNumInputArg = 3;
156         meth = 16;
157     }
158     else if(wcscmp(wcsMeth, L"def") == 0) // uniform (def)
159     {
160         iNumInputArg = 3;
161         meth = 17;
162     }
163     else if(wcscmp(wcsMeth, L"unf") == 0) // uniform (unf)
164     {
165         meth = 18;
166     }
167     else if(wcscmp(wcsMeth, L"uin") == 0) // uniform (uin)
168     {
169         meth = 19;
170     }
171     else if(wcscmp(wcsMeth, L"lgi") == 0) // uniform (lgi)
172     {
173         iNumInputArg = 3;
174         meth = 20;
175     }
176     else if(wcscmp(wcsMeth, L"getgen") == 0) // getgen
177     {
178         iNumInputArg = 1;
179         meth = 21;
180     }
181     else if(wcscmp(wcsMeth, L"setgen") == 0) // setgen
182     {
183         iNumInputArg = 2;
184         meth = 22;
185     }
186     else if(wcscmp(wcsMeth, L"getsd") == 0) // getsd
187     {
188         iNumInputArg = 1;
189         meth = 23;
190     }
191     else if(wcscmp(wcsMeth, L"setsd") == 0) // setsd
192     {
193         switch(current_base_gen)
194         {
195             case MT:        iNumInputArg = 2; break;
196             case KISS:      iNumInputArg = 5; break;
197             case CLCG2:     iNumInputArg = 3; break;
198             case CLCG4:     iNumInputArg = 5; break;
199             case URAND:     iNumInputArg = 2; break;
200             case FSULTRA:
201             {
202                 if(in.size() != 2 && in.size() != 3)
203                 {
204                     char* pstMeth = wide_string_to_UTF8(wcsMeth);
205                     Scierror(77, _("%s: Wrong number of input argument(s) for method %s: %d or %d expected.\n"), "grand", pstMeth, 2, 3);
206                     FREE(pstMeth);
207                     return types::Function::Error;
208                 }
209
210                 iNumInputArg = in.size();
211                 break;
212             }
213         }
214         meth = 24;
215     }
216     else if(wcscmp(wcsMeth, L"phr2sd") == 0) // phr2sd
217     {
218         iNumInputArg = 2;
219         meth = 25;
220     }
221     else if(wcscmp(wcsMeth, L"setcgn") == 0) // setcgn
222     {
223         iNumInputArg = 2;
224         meth = 26;
225     }
226     else if(wcscmp(wcsMeth, L"getcgn") == 0) // getcgn
227     {
228         iNumInputArg = 1;
229         meth = 27;
230     }
231     else if(wcscmp(wcsMeth, L"initgn") == 0) // initgn
232     {
233         iNumInputArg = 2;
234         meth = 28;
235     }
236     else if(wcscmp(wcsMeth, L"setall") == 0) // setall
237     {
238         iNumInputArg = 5;
239         meth = 29;
240     }
241     else if(wcscmp(wcsMeth, L"advnst") == 0) // advnst
242     {
243         iNumInputArg = 2;
244         meth = 30;
245     }
246     else
247     {
248         char* pstMeth = wide_string_to_UTF8(wcsMeth);
249         Scierror(78, _("%s: Wrong method for input argument #%d: %s unknown.\n"), "grand", iStrPos+1, pstMeth);
250         FREE(pstMeth);
251         return types::Function::Error;
252     }
253
254 // *** get arguments before methode string. ***
255     if(meth == 11 || meth == 13 || meth == 14 || meth == 16) // grand(n, "...", ...
256     {
257         if(iStrPos != 1)
258         {
259             Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d.\n"), "grand", iStrPos+1, 2);
260             return types::Function::Error;
261         }
262
263         if(in[iPos]->isDouble() == false)
264         {
265             Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "grand", iPos+1);
266             return types::Function::Error;
267         }
268
269         types::Double* pDblTemp = in[iPos]->getAs<types::Double>();
270
271         if(pDblTemp->isScalar() == false)
272         {
273             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos+1);
274             return types::Function::Error;
275         }
276
277         iNumIter = pDblTemp->get(0);
278         iPos++;
279     }
280     else if(meth < 21) // grand(m, n, "...", ... || grand(matrix, "...", ...
281     {
282         if(iStrPos != 1 && iStrPos != 2)
283         {
284             Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d or %d.\n"), "grand", iStrPos+1, 2, 3);
285             return types::Function::Error;
286         }
287
288         std::vector<types::Double*> vectpDblTemp;
289         for(int i = 0; i < iStrPos; i++)
290         {
291             if(in[iPos]->isDouble() == false)
292             {
293                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos+1);
294                 return types::Function::Error;
295             }
296
297             vectpDblTemp.push_back(in[iPos]->getAs<types::Double>());
298
299             if(iStrPos == 3 && vectpDblTemp[i]->isScalar() == false)
300             {
301                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos+1);
302                 return types::Function::Error;
303             }
304             iPos++;
305         }
306
307         if(iStrPos == 2)
308         {
309             iRows = static_cast<int>(vectpDblTemp[0]->get(0));
310             iCols = static_cast<int>(vectpDblTemp[1]->get(0));;
311         }
312         else // iStrPos == 2
313         {
314             iRows = vectpDblTemp[0]->getRows();
315             iCols = vectpDblTemp[0]->getCols();
316             iNumInputArg--;
317         }
318     }
319
320     if(in[iPos]->isString() == false)
321     {
322         Scierror(999, _("%s: Wrong type for input argument #%d : A string expected.\n"), "grand", iPos+1);
323         return types::Function::Error;
324     }
325     iPos++; // method string has been already got.
326
327 // *** check the number of input args according the methode. ***
328     if(in.size() != iNumInputArg)
329     {
330         char* pstMeth = wide_string_to_UTF8(wcsMeth);
331         Scierror(77, _("%s: Wrong number of input argument(s) for method %s: %d expected.\n"), "grand", pstMeth, iNumInputArg);
332         FREE(pstMeth);
333         return types::Function::Error;
334     }
335
336 // *** get arguments after methode string. ***
337     if(meth == 22 || meth == 25) // setgen || phr2sd
338     {
339         if(in[iPos]->isString() == false)
340         {
341             Scierror(999, _("%s: Wrong type for input argument #%d : A string expected.\n"), "grand", iPos+1);
342             return types::Function::Error;
343         }
344
345         pStrGenOrPhr = in[iPos]->getAs<types::String>();
346
347         if(pStrGenOrPhr->isScalar() == false)
348         {
349             Scierror(999, _("%s: Wrong size for input argument #%d : Only one string expected.\n"), "grand", iPos+1);
350             return types::Function::Error;
351         }
352     }
353     else
354     {
355         for(int i = iPos; i < in.size(); i++)
356         {
357             if(in[i]->isDouble() == false)
358             {
359                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos+1);
360                 return types::Function::Error;
361             }
362
363             vectpDblInput.push_back(in[i]->getAs<types::Double>());
364         }
365     }
366
367 // *** perform operation in according method string and return result. ***
368
369     switch(meth)
370     {
371         case 1: // beta
372         {
373             types::Double* pDblOut = new types::Double(iRows, iCols);
374             double minlog   = 1.e-37;
375
376             for(int i = 0; i < 2; i++)
377             {
378                 if(vectpDblInput[i]->isScalar() == false)
379                 {
380                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
381                     return types::Function::Error;
382                 }
383
384                 if(vectpDblInput[i]->get(0) < minlog)
385                 {
386                     Scierror(999, _("%s: Wrong value for input argument #%d : At least %lf expected.\n"), "grand", iPos+1, minlog);
387                     return types::Function::Error;
388                 }
389             }
390
391             for(int i = 0; i < pDblOut->getSize(); i++)
392             {
393                 pDblOut->set(i, C2F(genbet)(vectpDblInput[0]->get(),vectpDblInput[1]->get()));
394             }
395
396             out.push_back(pDblOut);
397             break;
398         }
399         case 2: // binomial
400         case 3: // negative binomial
401         {
402             types::Double* pDblOut = new types::Double(iRows, iCols);
403
404             for(int i = 0; i < 2; i++)
405             {
406                 if(vectpDblInput[i]->isScalar() == false)
407                 {
408                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
409                     return types::Function::Error;
410                 }
411             }
412
413             if(vectpDblInput[0]->get(0) < 0.0) // N
414             {
415                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive integer expected.\n"), "grand", 4);
416                 return types::Function::Error;
417             }
418
419             if(vectpDblInput[1]->get(0) < 0.0 || vectpDblInput[1]->get(0) > 1.0) // p
420             {
421                 Scierror(999, _("%s: Wrong value for input argument #%d : A value  expected.\n"), "grand", 5);
422                 return types::Function::Error;
423             }
424
425             int N = static_cast<int>(vectpDblInput[0]->get(0));
426             if(meth == 2)
427             {
428                 for(int i = 0; i < pDblOut->getSize(); i++)
429                 {
430                     pDblOut->set(i, static_cast<double>(C2F(ignbin)(&N, vectpDblInput[1]->get())));
431                 }
432             }
433             else // meth == 3
434             {
435                 for(int i = 0; i < pDblOut->getSize(); i++)
436                 {
437                     pDblOut->set(i, static_cast<double>(C2F(ignnbn)(&N, vectpDblInput[1]->get())));
438                 }
439             }
440
441             out.push_back(pDblOut);
442             break;
443         }
444         case 4: // chisquare
445         {
446             types::Double* pDblOut = new types::Double(iRows, iCols);
447
448             if(vectpDblInput[0]->isScalar() == false)
449             {
450                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"),"grand", 4);
451                 return types::Function::Error;
452             }
453
454             if(vectpDblInput[0]->get(0) <= 0.0) // Df
455             {
456                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", 4);
457                 return types::Function::Error;
458             }
459
460             for(int i = 0; i < pDblOut->getSize(); i++)
461             {
462                 pDblOut->set(i, C2F(genchi)(vectpDblInput[0]->get()));
463             }
464
465             out.push_back(pDblOut);
466             break;
467         }
468         case 5: // non central chisquare
469         {
470             types::Double* pDblOut = new types::Double(iRows, iCols);
471
472             for(int i = 0; i < 2; i++)
473             {
474                 if(vectpDblInput[i]->isScalar() == false)
475                 {
476                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
477                     return types::Function::Error;
478                 }
479             }
480
481             if(vectpDblInput[0]->get(0) < 1.0) // Df
482             {
483                 Scierror(999, _("%s: Wrong value for input argument #%d : value >1.0 expected.\n"), "grand", 4);
484                 return types::Function::Error;
485             }
486
487             if(vectpDblInput[1]->get(0) < 0.0) // Xnon
488             {
489                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
490                 return types::Function::Error;
491             }
492
493             for(int i = 0; i < pDblOut->getSize(); i++)
494             {
495                 pDblOut->set(i, C2F(gennch)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
496             }
497
498             out.push_back(pDblOut);
499             break;
500         }
501         case 6: // exponential
502         {
503             types::Double* pDblOut = new types::Double(iRows, iCols);
504
505             if(vectpDblInput[0]->isScalar() == false)
506             {
507                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
508                 return types::Function::Error;
509             }
510
511             if(vectpDblInput[0]->get(0) < 0.0) // Av
512             {
513                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
514                 return types::Function::Error;
515             }
516
517             for(int i = 0; i < pDblOut->getSize(); i++)
518             {
519                 pDblOut->set(i, C2F(genexp)(vectpDblInput[0]->get()));
520             }
521
522             out.push_back(pDblOut);
523             break;
524         }
525         case 7: // F variance ratio
526         {
527             types::Double* pDblOut = new types::Double(iRows, iCols);
528
529             for(int i = 0; i < 2; i++)
530             {
531                 if(vectpDblInput[i]->isScalar() == false)
532                 {
533                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
534                     return types::Function::Error;
535                 }
536
537                 if(vectpDblInput[i]->get(0) <= 0.0) // Dfn Dfd
538                 {
539                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", i+4);
540                     return types::Function::Error;
541                 }
542             }
543
544             for(int i = 0; i < pDblOut->getSize(); i++)
545             {
546                 pDblOut->set(i, C2F(genf)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
547             }
548
549             out.push_back(pDblOut);
550             break;
551         }
552         case 8: // non F variance ratio
553         {
554             types::Double* pDblOut = new types::Double(iRows, iCols);
555
556             for(int i = 0; i < 3; i++)
557             {
558                 if(vectpDblInput[i]->isScalar() == false)
559                 {
560                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
561                     return types::Function::Error;
562                 }
563             }
564
565             if(vectpDblInput[0]->get(0) < 1.0) // Dfn
566             {
567                 Scierror(999, _("%s: Wrong value for input argument #%d : value > 1.0 expected.\n"), "grand", 4);
568                 return types::Function::Error;
569             }
570
571             if(vectpDblInput[1]->get(0) <= 0.0) // Dfd
572             {
573                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", 5);
574                 return types::Function::Error;
575             }
576
577             if(vectpDblInput[2]->get(0) < 0.0) // Xnon
578             {
579                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 6);
580                 return types::Function::Error;
581             }
582
583             for(int i = 0; i < pDblOut->getSize(); i++)
584             {
585                 pDblOut->set(i, C2F(gennf)(vectpDblInput[0]->get(), vectpDblInput[1]->get(), vectpDblInput[2]->get()));
586             }
587
588             out.push_back(pDblOut);
589             break;
590         }
591         case 9: // gamma
592         {
593             types::Double* pDblOut = new types::Double(iRows, iCols);
594
595             for(int i = 0; i < 2; i++)
596             {
597                 if(vectpDblInput[i]->isScalar() == false)
598                 {
599                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
600                     return types::Function::Error;
601                 }
602
603                 if(vectpDblInput[i]->get(0) <= 0.0)
604                 {
605                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", i+4);
606                     return types::Function::Error;
607                 }
608             }
609
610             for(int i = 0; i < pDblOut->getSize(); i++)
611             {
612                 //** WARNING ** : order is changed in parameters for
613                 //compatibility between Rand(...'gam',..) and cdfgam
614
615                 pDblOut->set(i, C2F(gengam)(vectpDblInput[1]->get(), vectpDblInput[0]->get()));
616             }
617
618             out.push_back(pDblOut);
619             break;
620         }
621         case 10: // Gauss Laplace (normal)
622         {
623             types::Double* pDblOut = new types::Double(iRows, iCols);
624
625             for(int i = 0; i < 2; i++)
626             {
627                 if(vectpDblInput[i]->isScalar() == false)
628                 {
629                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
630                     return types::Function::Error;
631                 }
632             }
633
634             if(vectpDblInput[1]->get(0) < 0.0)
635             {
636                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
637                 return types::Function::Error;
638             }
639
640             for(int i = 0; i < pDblOut->getSize(); i++)
641             {
642                 pDblOut->set(i, C2F(gennor)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
643             }
644
645             out.push_back(pDblOut);
646             break;
647         }
648         case 11: // multivariate gaussian (multivariate normal)
649         {
650             if(vectpDblInput[0]->getCols() != 1 || vectpDblInput[0]->getSize() == 0)
651             {
652                 Scierror(999, _("%s: Wrong size for input argument #%d : A matrix of size m x 1 expected.(m > 0)\n"), "grand", 4);
653                 return types::Function::Error;
654             }
655
656             if(vectpDblInput[0]->getRows() != vectpDblInput[1]->getRows())
657             {
658                 Scierror(999, _("%s: Wrong size for input argument #%d and #%d: Mean and Cov have incompatible dimensions.\n"), "grand", 4, 5);
659                 return types::Function::Error;
660             }
661
662             if(vectpDblInput[1]->getRows() != vectpDblInput[1]->getCols())
663             {
664                 Scierror(999, _("%s: Wrong size for input argument #%d : A square symmetric positive definite matrix expected.\n"), "grand", 5);
665                 return types::Function::Error;
666             }
667
668             int ierr = 0;
669             int size = vectpDblInput[0]->getRows();
670             int mp   = (size * (size + 3)) / 2 + 1;
671
672             types::Double* pDblOut = new types::Double(size, iNumIter);
673             double* work  = (double*)malloc(size * sizeof(double));
674             double* param = (double*)malloc(mp   * sizeof(double));
675
676             types::Double* pDblMean = vectpDblInput[0]->clone()->getAs<types::Double>();
677             types::Double* pDblCov  = vectpDblInput[1]->clone()->getAs<types::Double>();
678
679             C2F(setgmn)(pDblMean->get(), pDblCov->get(), &size, &size, param, &ierr);
680
681             delete pDblMean;
682             delete pDblCov;
683
684             if(ierr == 1)
685             {
686                 free(work);
687                 free(param);
688                 Scierror(999, _("%s: setgmn return with state %d.\n"), "grand", ierr);
689                 return types::Function::Error;
690             }
691
692             for(int i = 0; i < iNumIter; i++)
693             {
694                 C2F(genmn)(param, pDblOut->get() + (size * i), work);
695             }
696
697             free(work);
698             free(param);
699
700             out.push_back(pDblOut);
701             break;
702         }
703         case 12: // geometric
704         {
705             types::Double* pDblOut = new types::Double(iRows, iCols);
706             double pmin = 1.3e-307;
707
708             if(vectpDblInput[0]->isScalar() == false)
709             {
710                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
711                 return types::Function::Error;
712             }
713
714             if(vectpDblInput[0]->get(0) < pmin || vectpDblInput[0]->get(0) > 1.0)
715             {
716                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %lf and %d.\n"), "grand", 4, pmin, 1);
717                 return types::Function::Error;
718             }
719
720             for(int i = 0; i < pDblOut->getSize(); i++)
721             {
722                 pDblOut->set(i, igngeom(vectpDblInput[0]->get(0)));
723             }
724
725             out.push_back(pDblOut);
726             break;
727         }
728         case 13: // markov
729         {
730             if( vectpDblInput[0]->getRows() != vectpDblInput[0]->getCols() &&
731                 vectpDblInput[0]->getRows() != 1)
732             {
733                 Scierror(999, _("%s: Wrong size for input argument #%d : A square matrix or a row vector expected.\n"), "grand", 4);
734                 return types::Function::Error;
735             }
736
737             if(vectpDblInput[1]->getSize() == 0)
738             {
739                 Scierror(999, _("%s: Wrong size for input argument #%d: No empty matrix expected.\n"), "grand", 5);
740                 return types::Function::Error;
741             }
742
743             int sizeOfX0 = vectpDblInput[1]->getSize();
744
745             for( int i = 0; i < sizeOfX0; i++)
746             {
747                 if(vectpDblInput[1]->get(i) < 1 || (vectpDblInput[1]->get(i) - 1) >= vectpDblInput[0]->getCols())
748                 {
749                     Scierror(999, _("%s: X0(%d) must be in the range [1,%d[.\n"), "grand", i+1, vectpDblInput[0]->getCols() + 1);
750                     return types::Function::Error;
751                 }
752             }
753
754             types::Double* pDblOut = new types::Double(sizeOfX0, iNumIter);
755             int size = vectpDblInput[0]->getSize() + vectpDblInput[0]->getCols();
756             double* work = (double*)malloc(size * sizeof(double));
757
758             for(int i = 0; i < vectpDblInput[0]->getRows(); i++)
759             {
760                 double ptot = 0.0;
761                 for(int j = 0; j < vectpDblInput[0]->getCols(); j++)
762                 {
763                     int position = vectpDblInput[0]->getRows() * j + i;
764
765                     if(vectpDblInput[0]->get(position) < 0 || vectpDblInput[0]->get(position) > 1)
766                     {
767                         Scierror(999, _("%s: Wrong value for input argument #%d: P(%d,%d) must be in the range [0 1].\n"), "grand", i+1, j+1);
768                         return types::Function::Error;
769                     }
770
771                     ptot += vectpDblInput[0]->get(position);
772                 }
773
774                 if(fabs(ptot - 1.0) > 1e-8)
775                 {
776                     Scierror(999, _("%s: Sum of P(%d,1:%d)=%lf ~= 1.\n"), "grand", i+1, vectpDblInput[0]->getCols(), ptot);
777                     return types::Function::Error;
778                 }
779             }
780             // ** Computing the cumulative sum of the P matrix **
781
782             for(int i = 0; i < vectpDblInput[0]->getRows(); i++)
783             {
784                 double cumsum = 0.0;
785                 work[i] = 0.0;
786                 for(int j = 1; j < vectpDblInput[0]->getCols() + 1; j++)
787                 {
788                     cumsum += vectpDblInput[0]->get(vectpDblInput[0]->getRows() * (j - 1) + i);
789                     work[vectpDblInput[0]->getRows() * j + i] = cumsum;
790                 }
791             }
792
793             for(int i = 0; i < vectpDblInput[1]->getSize(); i++)
794             {
795                 int iCur = static_cast<int>(vectpDblInput[1]->get(i)) - 1;
796                 for(int j = 0; j < iNumIter; j++)
797                 {
798                     int niv = 0;
799                     double rr = C2F(ranf)();
800
801                     if(vectpDblInput[0]->getRows() == 1)
802                     {
803                         iCur = 0;
804                     }
805
806                     while(  rr >= work[iCur + vectpDblInput[0]->getRows() * niv] &&
807                             niv < (vectpDblInput[0]->getRows() + 1))
808                     {
809                         niv++;
810                     }
811
812                     // ** projection to avoid boundaries **
813                     niv = Max(Min(niv, vectpDblInput[0]->getCols()), 1);
814                     pDblOut->set(vectpDblInput[1]->getSize() * j + i, static_cast<double>(niv));
815                     iCur = niv - 1;
816                 }
817             }
818
819             out.push_back(pDblOut);
820             break;
821         }
822         case 14: // multinomial
823         {
824             if(vectpDblInput[0]->isScalar() == false)
825             {
826                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "grand", 4);
827                 return types::Function::Error;
828             }
829
830             if(vectpDblInput[0]->get(0) < 0)
831             {
832                 Scierror(999, _("%s: Wrong value for input argument #%d: A positive scalar expected.\n"), "grand", 4);
833                 return types::Function::Error;
834             }
835
836             if(vectpDblInput[1]->getCols() != 1 || vectpDblInput[1]->getRows() <= 0)
837             {
838                 Scierror(999, _("%s: Wrong size for input argument #%d: A colomn vector expected.\n"), "grand", 5);
839                 return types::Function::Error;
840             }
841
842             double ptot = 0.0;
843             int ncat = vectpDblInput[1]->getRows() + 1;
844             types::Double* pDblOut = new types::Double(ncat, iNumIter);
845
846             for(int i = 0; i < vectpDblInput[1]->getCols(); i++)
847             {
848                 if(vectpDblInput[1]->get(i) < 0.0 || vectpDblInput[1]->get(i) > 1.0)
849                 {
850                     Scierror(999, _("%s: Wrong value for input argument #%d: P(%d) must be in the range [0 1].\n"), "grand", 5, i+1);
851                     return types::Function::Error;
852                 }
853                 ptot += vectpDblInput[1]->get(i);
854             }
855
856             if(ptot > 1.0)
857             {
858                 Scierror(999, _("%s: Wrong value for input argument #%d: Sum of P(i) > 1.\n"), "grand", 5);
859                 return types::Function::Error;
860             }
861
862             int* piP    = (int*)malloc(vectpDblInput[0]->getSize() * sizeof(int));
863             int* piOut  = (int*)malloc(pDblOut->getSize() * sizeof(int));
864
865             for(int i = 0; i < vectpDblInput[0]->getSize(); i++)
866             {
867                 piP[i] = static_cast<int>(vectpDblInput[0]->get(i));
868             }
869
870             for(int i = 0; i < iNumIter; i++)
871             {
872                 C2F(genmul)(piP, vectpDblInput[1]->get(), &ncat,  + &piOut[ncat * i]);
873             }
874
875             for(int i = 0; i < pDblOut->getSize(); i++)
876             {
877                 pDblOut->set(i, static_cast<double>(piOut[i]));
878             }
879
880             free(piP);
881             free(piOut);
882
883             out.push_back(pDblOut);
884             break;
885         }
886         case 15: // Poisson
887         {
888             types::Double* pDblOut = new types::Double(iRows, iCols);
889
890             if(vectpDblInput[0]->isScalar() == false)
891             {
892                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
893                 return types::Function::Error;
894             }
895
896             if(vectpDblInput[0]->get(0) < 0.0)
897             {
898                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
899                 return types::Function::Error;
900             }
901
902             for(int i = 0; i < pDblOut->getSize(); i++)
903             {
904                 pDblOut->set(i, static_cast<double>(C2F(ignpoi)(vectpDblInput[0]->get())));
905             }
906
907             out.push_back(pDblOut);
908             break;
909         }
910         case 16: // random permutations
911         {
912             if(vectpDblInput[0]->getCols() != 1)
913             {
914                 Scierror(999, _("%s: Wrong type for input argument #%d : A colomn vector expected.\n"), "grand", 4);
915                 return types::Function::Error;
916             }
917
918             int iVectRows = vectpDblInput[0]->getRows();
919             types::Double* pDblOut = new types::Double(iVectRows, iNumIter);
920
921             for(int i = 0; i < iNumIter; i++)
922             {
923                 for(int j = 0; j < iVectRows; j++)
924                 {
925                     pDblOut->set(iVectRows * i + j, vectpDblInput[0]->get(j));
926                 }
927                 C2F(genprm)(pDblOut->get() + (iVectRows * i), &iVectRows);
928             }
929
930             out.push_back(pDblOut);
931             break;
932         }
933         case 17: // uniform (def)
934         {
935             types::Double* pDblOut = new types::Double(iRows, iCols);
936
937             for(int i = 0; i < pDblOut->getSize(); i++)
938             {
939                 pDblOut->set(i, C2F(ranf)());
940             }
941
942             out.push_back(pDblOut);
943             break;
944         }
945         case 18: // uniform (unf)
946         {
947             types::Double* pDblOut = new types::Double(iRows, iCols);
948
949             for(int i = 0; i < 2; i++)
950             {
951                 if(vectpDblInput[i]->isScalar() == false)
952                 {
953                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
954                     return types::Function::Error;
955                 }
956             }
957
958             double low  = vectpDblInput[0]->get(0);
959             double high = vectpDblInput[1]->get(0);
960
961             if(low > high)
962             {
963                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
964                 return types::Function::Error;
965             }
966
967             for(int i = 0; i < pDblOut->getSize(); i++)
968             {
969                 pDblOut->set(i, low + (high - low) * C2F(ranf)());
970             }
971
972             out.push_back(pDblOut);
973             break;
974         }
975         case 19: // uniform (uin)
976         {
977             types::Double* pDblOut = new types::Double(iRows, iCols);
978
979             for(int i = 0; i < 2; i++)
980             {
981                 if(vectpDblInput[i]->isScalar() == false)
982                 {
983                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+4);
984                     return types::Function::Error;
985                 }
986             }
987
988             int low  = static_cast<int>(vectpDblInput[0]->get(0));
989             int high = static_cast<int>(vectpDblInput[1]->get(0));
990
991             if(low > high)
992             {
993                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
994                 return types::Function::Error;
995             }
996
997             if( low  != vectpDblInput[0]->get(0) ||
998                 high != vectpDblInput[1]->get(0) ||
999                 (high - low + 1) > 2147483561)
1000             {
1001                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low and High must be integers and (high - low + 1) <=  2147483561.\n"), "grand", 4, 5);
1002                 return types::Function::Error;
1003             }
1004
1005             for(int i = 0; i < pDblOut->getSize(); i++)
1006             {
1007                 pDblOut->set(i, C2F(ignuin)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
1008             }
1009
1010             out.push_back(pDblOut);
1011             break;
1012         }
1013         case 20: // uniform (lgi)
1014         {
1015             types::Double* pDblOut = new types::Double(iRows, iCols);
1016
1017             for(int i = 0; i < pDblOut->getSize(); i++)
1018             {
1019                 pDblOut->set(i, static_cast<double>(ignlgi()));
1020             }
1021
1022             out.push_back(pDblOut);
1023             break;
1024         }
1025         case 21: // getgen
1026         {
1027             types::String* pStrOut = new types::String(names_gen[current_base_gen]);
1028             out.push_back(pStrOut);
1029             break;
1030         }
1031         case 22: // setgen
1032         {
1033             wchar_t* wcsGen = pStrGenOrPhr->get(0);
1034
1035             if(wcscmp(wcsGen, L"mt") == 0)
1036             {
1037                 ConfigVariable::setCurrentBaseGen(MT);
1038             }
1039             else if(wcscmp(wcsGen, L"kiss") == 0)
1040             {
1041                 ConfigVariable::setCurrentBaseGen(KISS);
1042             }
1043             else if(wcscmp(wcsGen, L"clcg4") == 0)
1044             {
1045                 ConfigVariable::setCurrentBaseGen(CLCG4);
1046             }
1047             else if(wcscmp(wcsGen, L"clcg2") == 0)
1048             {
1049                 ConfigVariable::setCurrentBaseGen(CLCG2);
1050             }
1051             else if(wcscmp(wcsGen, L"urand") == 0)
1052             {
1053                 ConfigVariable::setCurrentBaseGen(URAND);
1054             }
1055             else if(wcscmp(wcsGen, L"fsultra") == 0)
1056             {
1057                 ConfigVariable::setCurrentBaseGen(FSULTRA);
1058             }
1059             else
1060             {
1061                 Scierror(999, _("%s: Wrong value for input argument #%d: '%s', '%s', '%s', '%s', '%s' or '%s' expected.\n"), "grand", 2, "mt", "kiss", "clcg4", "clcg2", "urand", "fsultra");
1062                 return types::Function::Error;
1063             }
1064
1065             int current_gen = ConfigVariable::getCurrentBaseGen();
1066             types::String* pStrOut = new types::String(names_gen[current_gen]);
1067             out.push_back(pStrOut);
1068             break;
1069         }
1070         case 23: // getsd
1071         {
1072             types::Double* pDblOut = NULL;
1073
1074             switch(current_base_gen)
1075             {
1076                 case MT:
1077                 {
1078                     pDblOut = new types::Double(625, 1);
1079                     get_state_mt(pDblOut->get());
1080                     break;
1081                 }
1082                 case KISS:
1083                 {
1084                     pDblOut = new types::Double(4, 1);
1085                     get_state_kiss(pDblOut->get());
1086                     break;
1087                 }
1088                 case CLCG4:
1089                 {
1090                     pDblOut = new types::Double(4, 1);
1091                     int current_clcg4 = ConfigVariable::getCurrentClcg4();
1092                     get_state_clcg4(current_clcg4, pDblOut->get());
1093                     break;
1094                 }
1095                 case CLCG2:
1096                 {
1097                     pDblOut = new types::Double(2, 1);
1098                     get_state_clcg2(pDblOut->get());
1099                     break;
1100                 }
1101                 case URAND:
1102                 {
1103                     pDblOut = new types::Double(1, 1);
1104                     get_state_urand(pDblOut->get());
1105                     break;
1106                 }
1107                 case FSULTRA:
1108                 {
1109                     pDblOut = new types::Double(40, 1);
1110                     get_state_fsultra(pDblOut->get());
1111                     break;
1112                 }
1113             }
1114
1115             out.push_back(pDblOut);
1116             break;
1117         }
1118         case 24: // setsd
1119         {
1120             int ierr = 0;
1121             switch(current_base_gen)
1122             {
1123                 case MT:
1124                 {
1125                     if(vectpDblInput[0]->isScalar())
1126                     {
1127                         ierr = set_state_mt_simple(vectpDblInput[0]->get(0));
1128                     }
1129                     else if(vectpDblInput[0]->getSize() != 625)
1130                     {
1131                         ierr = set_state_mt(vectpDblInput[0]->get());
1132                     }
1133                     else
1134                     {
1135                         Scierror(999, _("%s: Wrong size for input argument #%d : A scalar or a vector of size %d expected.\n"), "grand", 4, 625);
1136                         return types::Function::Error;
1137                     }
1138
1139                     break;
1140                 }
1141                 case KISS:
1142                 case CLCG4:
1143                 {
1144                     for(int i = 0; i < 4; i++)
1145                     {
1146                         if(vectpDblInput[i]->isScalar() == false)
1147                         {
1148                             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+2);
1149                             return types::Function::Error;
1150                         }
1151                     }
1152
1153                     if(current_base_gen == KISS)
1154                     {
1155                         ierr = set_state_kiss(  vectpDblInput[0]->get(0),
1156                                                 vectpDblInput[1]->get(0),
1157                                                 vectpDblInput[2]->get(0),
1158                                                 vectpDblInput[3]->get(0));
1159                     }
1160                     else // CLCG4
1161                     {
1162                         ierr = set_seed_clcg4(  ConfigVariable::getCurrentClcg4(),
1163                                                 vectpDblInput[0]->get(0),
1164                                                 vectpDblInput[1]->get(0),
1165                                                 vectpDblInput[2]->get(0),
1166                                                 vectpDblInput[3]->get(0));
1167                     }
1168
1169                     break;
1170                 }
1171                 case CLCG2:
1172                 {
1173                     for(int i = 0; i < 2; i++)
1174                     {
1175                         if(vectpDblInput[i]->isScalar() == false)
1176                         {
1177                             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+2);
1178                             return types::Function::Error;
1179                         }
1180                     }
1181
1182                     ierr = set_state_clcg2(vectpDblInput[0]->get(0), vectpDblInput[1]->get(0));
1183                     break;
1184                 }
1185                 case URAND:
1186                 {
1187                     if(vectpDblInput[0]->isScalar() == false)
1188                     {
1189                         Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1190                         return types::Function::Error;
1191                     }
1192
1193                     ierr = set_state_urand(vectpDblInput[0]->get(0));
1194                     break;
1195                 }
1196                 case FSULTRA:
1197                 {
1198                     if(in.size() == 2)
1199                     {
1200                         if(vectpDblInput[0]->getRows() != 40 || vectpDblInput[0]->getCols() != 1)
1201                         {
1202                             Scierror(999, _("%s: Wrong size for input argument #%d : A vector of size %d x %d expected.\n"), "grand", 2, 40, 1);
1203                             return types::Function::Error;
1204                         }
1205
1206                         ierr = set_state_fsultra(vectpDblInput[0]->get());
1207                     }
1208                     else // in.size() == 3
1209                     {
1210                         for(int i = 0; i < 2; i++)
1211                         {
1212                             if(vectpDblInput[i]->isScalar() == false)
1213                             {
1214                                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+2);
1215                                 return types::Function::Error;
1216                             }
1217                         }
1218
1219                         ierr = set_state_fsultra_simple(vectpDblInput[0]->get(0), vectpDblInput[1]->get(0));
1220                     }
1221
1222                     break;
1223                 }
1224             }
1225
1226             if(ierr == 0)
1227             {
1228                 Scierror(999, _("%s: Wrong value for the last %d input argument(s).\n"), "grand", in.size() - 1);
1229                 return types::Function::Error;
1230             }
1231
1232             break;
1233         }
1234         case 25: // phr2sd
1235         {
1236             if(pStrGenOrPhr->isScalar() == false)
1237             {
1238                 Scierror(999, _("%s: Wrong type for input argument #%d : One string expected.\n"), "grand", 2);
1239                 return types::Function::Error;
1240             }
1241
1242             types::Double* pDblOut = new types::Double(1, 2);
1243             int size = wcslen(pStrGenOrPhr->get(0));
1244             int piOut[2];
1245             char* strPhr = wide_string_to_UTF8(pStrGenOrPhr->get(0));
1246
1247             C2F(phrtsd)(strPhr, &size, piOut, piOut + 1, size);
1248
1249             pDblOut->set(0, static_cast<double>(piOut[0]));
1250             pDblOut->set(1, static_cast<double>(piOut[1]));
1251
1252             FREE(strPhr);
1253             out.push_back(pDblOut);
1254             break;
1255         }
1256         case 26: // setcgn
1257         {
1258             if(current_base_gen != CLCG4)
1259             {
1260                 sciprint(_("The %s option affects only the %s generator\n"), "setcgn", "clcg4");
1261             }
1262
1263             if(vectpDblInput[0]->isScalar() == false)
1264             {
1265                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1266                 return types::Function::Error;
1267             }
1268
1269             if(vectpDblInput[0]->get(0) < 0 || vectpDblInput[0]->get(0) > Maxgen)
1270             {
1271                Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %d and %d.\n"), "grand", 0, Maxgen);
1272                 return types::Function::Error;
1273             }
1274
1275             ConfigVariable::setCurrentClcg4(static_cast<int>(vectpDblInput[0]->get(0)));
1276             double dOut = static_cast<double>(ConfigVariable::getCurrentClcg4());
1277             out.push_back(new types::Double(dOut));
1278             break;
1279         }
1280         case 27: // getcgn
1281         {
1282             double dOut = static_cast<double>(ConfigVariable::getCurrentClcg4());
1283             out.push_back(new types::Double(dOut));
1284             break;
1285         }
1286         case 28: // initgn
1287         {
1288             SeedType where;
1289             if(current_base_gen != CLCG4)
1290             {
1291                 sciprint(_("The %s option affects only the %s generator\n"), "initgn", "clcg4");
1292             }
1293
1294             if(vectpDblInput[0]->isScalar() == false)
1295             {
1296                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1297                 return types::Function::Error;
1298             }
1299
1300             if( vectpDblInput[0]->get(0) != 0 &&
1301                 vectpDblInput[0]->get(0) != -1 &&
1302                 vectpDblInput[0]->get(0) != 1)
1303             {
1304                Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %d, %d or %d.\n"), "grand", 2, -1, 0, 1);
1305                 return types::Function::Error;
1306             }
1307
1308             where = (SeedType)(int)(vectpDblInput[0]->get(0) + 1);
1309             init_generator_clcg4(ConfigVariable::getCurrentClcg4(), where);
1310             out.push_back(vectpDblInput[0]);
1311             break;
1312         }
1313         case 29: // setall
1314         {
1315             if(current_base_gen != CLCG4)
1316             {
1317                 sciprint(_("The %s option affects only the %s generator\n"), "setall", "clcg4");
1318             }
1319
1320             for(int i = 0; i < 4; i++)
1321             {
1322                 if(vectpDblInput[i]->isScalar() == false)
1323                 {
1324                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i+2);
1325                     return types::Function::Error;
1326                 }
1327             }
1328
1329             int ierr = set_initial_seed_clcg4(  vectpDblInput[0]->get(0),
1330                                                 vectpDblInput[1]->get(0),
1331                                                 vectpDblInput[2]->get(0),
1332                                                 vectpDblInput[3]->get(0));
1333             if(ierr == 0)
1334             {
1335                 Scierror(999, _("%s: Wrong value for the last %d input argument(s).\n"), "grand", 4);
1336                 return types::Function::Error;
1337             }
1338
1339             out.push_back(pStrMethod);
1340             break;
1341         }
1342         case 30: // advnst
1343         {
1344             if(current_base_gen != CLCG4)
1345             {
1346                 sciprint(_("The %s option affects only the %s generator\n"), "advnst", "clcg4");
1347             }
1348
1349             if(vectpDblInput[0]->isScalar() == false)
1350             {
1351                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1352                 return types::Function::Error;
1353             }
1354
1355             int k = static_cast<int>(vectpDblInput[0]->get(0));
1356
1357             if(k < 1)
1358             {
1359                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be > %d.\n"), "grand", 2, 0);
1360                 return types::Function::Error;
1361             }
1362
1363             advance_state_clcg4(ConfigVariable::getCurrentClcg4(), k);
1364             out.push_back(new types::Double(static_cast<double>(k)));
1365             break;
1366         }
1367     }
1368
1369     return types::Function::OK;
1370 }
1371 /*--------------------------------------------------------------------------*/
1372