1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "monte_carlo.h"
22 #include "gretl_func.h"
23 #include "objstack.h"
24 #include "cmd_private.h"
25 #include "libset.h"
26 #include "gretl_mt.h"
27 #include "uservar.h"
28 #include "gretl_panel.h"
29 #include "gretl_string_table.h"
30 #include "gretl_xml.h"
31 #include "forecast.h"
32 #include "gretl_typemap.h"
33
34 #ifdef USE_CURL
35 # include "gretl_www.h"
36 #endif
37
38 #ifdef WIN32
39 # include "gretl_win32.h"
40 #endif
41
42 #ifdef _OPENMP
43 # include <omp.h>
44 #endif
45
46 #ifdef HAVE_MPI
47 # include "gretl_mpi.h"
48 #endif
49
50 #if defined(HAVE_MPI) || defined(USE_RLIB)
51 # include "gretl_foreign.h"
52 #endif
53
54 #if HAVE_GMP
55 # include <gmp.h>
56 #endif
57
58 #include <errno.h>
59
60 #if defined(_OPENMP) || defined(HAVE_MPI) || defined(WIN32)
61 # define WANT_XTIMER
62 #endif
63 #if !defined(_OPENMP) && !defined(WIN32)
64 # define NEED_ITIMER
65 #endif
66
67 static void gretl_tests_cleanup (void);
68
69 /**
70 * date_as_double:
71 * @t: observation number (zero-based).
72 * @pd: data periodicity or frequency.
73 * @sd0: floating point representation of starting date.
74 *
75 * Returns: the date corresponding to @t, as a double-precision number.
76 */
77
date_as_double(int t,int pd,double sd0)78 double date_as_double (int t, int pd, double sd0)
79 {
80 int ysd, yy, pp, yp;
81 int p10 = 10;
82
83 ysd = (int) sd0;
84
85 if (pd == 1) {
86 return (double) (ysd + t);
87 }
88
89 pp = pd;
90 while ((pp = pp / 10)) {
91 p10 *= 10;
92 }
93
94 pp = t % pd + p10 * (sd0 - ysd) + .5;
95 if (pp != pd) {
96 yy = ysd + t/pd + pp/pd + .5;
97 yp = pp % pd;
98 } else {
99 yy = ysd + t/pd + .5;
100 yp = pp;
101 }
102
103 return yy + (double) yp / p10;
104 }
105
106 /**
107 * ijton:
108 * @i: row number (0-based)
109 * @j: column number (0-based)
110 * @nrows: number of rows (and columns) in symmetric matrix.
111 *
112 * Given a (row, column) reference into a symmetric 2-dimensional
113 * matrix A, finds the index into a 1-dimensional array x
114 * composed of the non-redundant (lower) elements of A.
115 *
116 * E.g. for the 3 x 3 case with 6 non-redundant elements, 0 to 5,
117 *
118 * A(0,0) = x[0] A(0,1) = x[1] A(0,2) = x[2]
119 * A(1,0) = x[1] A(1,1) = x[3] A(1,2) = x[4]
120 * A(2,0) = x[2] A(2,1) = x[4] A(2,1) = x[5]
121 *
122 * Returns: 0-based index into flat array.
123 */
124
ijton(int i,int j,int nrows)125 int ijton (int i, int j, int nrows)
126 {
127 if (i > j) {
128 int tmp = i;
129
130 i = j;
131 j = tmp;
132 }
133
134 return nrows * i + j - i - ((i - 1) * i / 2);
135 }
136
137 /**
138 * transcribe_array:
139 * @targ: arrat to which to write.
140 * @src: array from which to read.
141 * @dset: data information struct.
142 *
143 * Copy from @src to @targ, skipping any missing values,
144 * over the sample range defined in @dset.
145 *
146 * Returns: the number of valid observations put into @targ.
147 */
148
transcribe_array(double * targ,const double * src,const DATASET * dset)149 int transcribe_array (double *targ, const double *src,
150 const DATASET *dset)
151 {
152 int t, n = 0;
153
154 for (t=dset->t1; t<=dset->t2; t++) {
155 if (!na(src[t])) {
156 targ[n++] = src[t];
157 }
158 }
159
160 return n;
161 }
162
163 /**
164 * gretl_isdummy:
165 * @t1: starting observation.
166 * @t2: ending observation.
167 * @x: data series to examine.
168 *
169 * Check whether series @x has only 0 or 1 values over the
170 * given sample range (or possibly missing values).
171 *
172 * Returns: 0 if the variable is not a 0/1 dummy, otherwise the
173 * number of 1s in the series.
174 */
175
gretl_isdummy(int t1,int t2,const double * x)176 int gretl_isdummy (int t1, int t2, const double *x)
177 {
178 int t, m = 0, goodobs = 0;
179
180 for (t=t1; t<=t2; t++) {
181 if (na(x[t])) {
182 continue;
183 }
184 if (x[t] != 0.0 && x[t] != 1.0) {
185 return 0;
186 }
187 if (x[t] == 1.0) {
188 m++;
189 }
190 goodobs++;
191 }
192
193 if (m < goodobs) {
194 return m;
195 }
196
197 return 0;
198 }
199
200 /**
201 * gretl_iszero:
202 * @x: data series to examine.
203 * @t1: starting observation.
204 * @t2: ending observation.
205 *
206 * Check whether series @x has only zero values over the
207 * given sample range (or possibly missing values).
208 *
209 * Returns: 1 if the variable is all zeros, otherwise 0.
210 */
211
gretl_iszero(int t1,int t2,const double * x)212 int gretl_iszero (int t1, int t2, const double *x)
213 {
214 double sum = 0.0;
215 int t;
216
217 for (t=t1; t<=t2; t++) {
218 if (!na(x[t])) {
219 sum += x[t] * x[t];
220 }
221 }
222
223 return floateq(sum, 0.0);
224 }
225
226 /**
227 * gretl_isconst:
228 * @x: data series to examine.
229 * @t1: starting observation.
230 * @t2: ending observation.
231 *
232 * Check whether series @x is constant over the
233 * given sample range (aside from any missing values).
234 *
235 * Returns: 1 if the variable is constant, otherwise 0.
236 */
237
gretl_isconst(int t1,int t2,const double * x)238 int gretl_isconst (int t1, int t2, const double *x)
239 {
240 int t, ret = 1;
241
242 while (na(x[t1]) && t1 <= t2) {
243 t1++;
244 }
245
246 if (t1 >= t2) {
247 return 0;
248 }
249
250 for (t=t1+1; t<=t2; t++) {
251 if (na(x[t])) {
252 continue;
253 }
254 if (floatneq(x[t], x[t1])) {
255 ret = 0;
256 break;
257 }
258 }
259
260 return ret;
261 }
262
263 /**
264 * gretl_isstoch:
265 * @x: data series to examine.
266 * @t1: starting observation.
267 * @t2: ending observation.
268 *
269 * Check whether series @x is stochastic, and contains a
270 * contiguous set of valid values within the given sample span.
271 * The simple and fallible heuristic for a stochastic series is
272 * that the differences between successive values are all the
273 * same.
274 *
275 * Returns: 1 if the variable appears to be stochastic on
276 * the specified criterion, otherwise 0.
277 */
278
gretl_isstoch(int t1,int t2,const double * x)279 int gretl_isstoch (int t1, int t2, const double *x)
280 {
281 double dx0;
282 int na_count = 0;
283 int multidiff = 0;
284 int s1 = -1, s2 = -1;
285 int t;
286
287 if (t1 >= t2) {
288 return 0;
289 }
290
291 for (t=t1; t<=t2; t++) {
292 if (!na(x[t])) {
293 s1 = t;
294 break;
295 }
296 }
297 for (t=t2; t>=s1; t--) {
298 if (!na(x[t])) {
299 s2 = t;
300 break;
301 }
302 }
303 if (s1 < 0 || s2 < 0 || s1 >= s2) {
304 return 0;
305 }
306
307 dx0 = x[s1+1] - x[s1];
308
309 for (t=s1; t<=s2; t++) {
310 if (na(x[t])) {
311 na_count++;
312 break;
313 }
314 if (t > s1 && !multidiff) {
315 if (x[t] - x[t-1] != dx0) {
316 multidiff = 1;
317 }
318 }
319 }
320
321 return (na_count == 0 && multidiff);
322 }
323
324 /**
325 * gretl_isunits:
326 * @x: data series to examine.
327 * @t1: starting observation.
328 * @t2: ending observation.
329 *
330 * Check whether series @x equals 1 over the
331 * given sample range (aside from any missing values).
332 *
333 * Returns: 1 if so, otherwise 0.
334 */
335
gretl_isunits(int t1,int t2,const double * x)336 int gretl_isunits (int t1, int t2, const double *x)
337 {
338 int t, ret = 1;
339
340 for (t=t1; t<=t2; t++) {
341 if (!na(x[t]) && x[t] != 1.0) {
342 ret = 0;
343 break;
344 }
345 }
346
347 return ret;
348 }
349
350 /**
351 * gretl_iscount:
352 * @x: data series to examine.
353 * @t1: starting observation.
354 * @t2: ending observation.
355 *
356 * Check whether series @x contains nothing but non-negative
357 * integer values (some of which are > 1) over the
358 * given sample range.
359 *
360 * Returns: 1 if so, otherwise 0.
361 */
362
gretl_iscount(int t1,int t2,const double * x)363 int gretl_iscount (int t1, int t2, const double *x)
364 {
365 int t, xi;
366 int g1 = 0;
367 int ret = 1;
368
369 for (t=t1; t<=t2; t++) {
370 if (na(x[t])) {
371 continue;
372 }
373 if (x[t] < 0.0) {
374 ret = 0;
375 break;
376 }
377 xi = x[t];
378 if (x[t] != (double) xi) {
379 ret = 0;
380 break;
381 }
382 if (x[t] > 1.0) {
383 g1 = 1;
384 }
385 }
386
387 if (g1 == 0) {
388 ret = 0;
389 }
390
391 return ret;
392 }
393
394 #define FEW_VALS 32
395 #define FEWER_VALS 8
396
few_vals(int t1,int t2,const double * x,double * ratio)397 static int few_vals (int t1, int t2, const double *x, double *ratio)
398 {
399 double test[FEW_VALS];
400 int match;
401 int i, t, n = 0, nv = 0;
402
403 for (t=t1; t<=t2; t++) {
404 if (!na(x[t])) {
405 match = 0;
406 for (i=0; i<nv; i++) {
407 if (x[t] == test[i]) {
408 /* we've already seen this value */
409 match = 1;
410 break;
411 }
412 }
413 if (!match) {
414 /* x[t] is a "new" value */
415 if (nv == FEW_VALS) {
416 /* hit the limit */
417 nv++;
418 break;
419 }
420 test[nv++] = x[t];
421 }
422 n++;
423 }
424 }
425
426 /* ratio of distinct values to observations */
427 *ratio = nv / (double) n;
428
429 return nv;
430 }
431
432 /**
433 * gretl_isdiscrete:
434 * @x: data series to examine.
435 * @t1: starting observation.
436 * @t2: ending observation.
437 *
438 * Checks the variable @x over the range @t1 to @t2 for discreteness.
439 * This is a heuristic whose components are (a) whether the values
440 * are "fairly round" (multiples of 0.25) or not, and, if test (a) is
441 * passed, (b) whether the variable takes on only "few" distinct
442 * values.
443 *
444 * Returns: 0 if test (a) is not passed or the number of distinct values
445 * is > 32; else 1 if the number of distinct values is <= 32; else 2 if
446 * the number of distinct values is <= 8. A return of 1 is supposed
447 * to indicate that it's "reasonable" to treat @x as discrete, while
448 * a return of 2 indicates that it's probably unreasonable _not_ to
449 * treat @x as discrete for the purpose of drawing up a frequency
450 * distribution.
451 */
452
gretl_isdiscrete(int t1,int t2,const double * x)453 int gretl_isdiscrete (int t1, int t2, const double *x)
454 {
455 int t, n = 0, disc = 1;
456 int allints = 1;
457 double r = 0;
458
459 for (t=t1; t<=t2; t++) {
460 if (na(x[t])) {
461 continue;
462 }
463 n++;
464 if (!ok_int(x[t])) {
465 allints = disc = 0;
466 break;
467 }
468 r = x[t] - floor(x[t]);
469 if (allints && r != 0) {
470 allints = 0;
471 }
472 if (r != 0.0 && r != 0.25 && r != 0.5 && r != 0.75) {
473 disc = 0;
474 break;
475 }
476 }
477
478 if (n == 0) {
479 disc = 0;
480 }
481
482 if (disc) {
483 n = few_vals(t1, t2, x, &r);
484 if (allints) {
485 if (n <= FEW_VALS && r < 0.2) {
486 disc = 2;
487 } else {
488 disc = (n <= FEWER_VALS)? 2 : 1;
489 }
490 } else if (n > FEW_VALS) {
491 disc = 0;
492 } else if (r > 0.9 && n > 30) {
493 /* somewhat arbitrary: but if r (= ratio of distinct
494 values to number of cases) is "too high", and the
495 number of observations is not tiny, perhaps we should
496 not automatically take the var as discrete
497 */
498 disc = 0;
499 } else if (n <= FEWER_VALS) {
500 disc = 2;
501 }
502 }
503
504 return disc;
505 }
506
accept_as_discrete(const DATASET * dset,int v,int strict)507 int accept_as_discrete (const DATASET *dset, int v, int strict)
508 {
509 if (series_is_discrete(dset, v)) {
510 /* the series has been explicitly marked as discrete */
511 return 1;
512 } else {
513 /* check for plausibility of discreteness */
514 int d = gretl_isdiscrete(dset->t1, dset->t2, dset->Z[v]);
515
516 return strict ? d > 1 : d > 0;
517 }
518 }
519
520 /**
521 * gretl_ispositive:
522 * @x: data series to examine.
523 * @t1: starting observation.
524 * @t2: ending observation.
525 * @strict: boolean, strict inequality.
526 *
527 * Returns: 1 if all non-missing values in @x, over the range @t1 to
528 * @t2, are positive (if strict==1) or non-negative (if strict==0),
529 * otherwise 0.
530 */
531
gretl_ispositive(int t1,int t2,const double * x,int strict)532 int gretl_ispositive (int t1, int t2, const double *x, int strict)
533 {
534 int t;
535
536 if (strict) {
537 for (t=t1; t<=t2; t++) {
538 if (x[t] <= 0) {
539 return 0;
540 }
541 }
542 } else {
543 for (t=t1; t<=t2; t++) {
544 if (x[t] < 0) {
545 return 0;
546 }
547 }
548 }
549
550 return 1;
551 }
552
553 /**
554 * gretl_is_oprobit_ok:
555 * @x: data series to examine.
556 * @t1: starting observation.
557 * @t2: ending observation.
558 *
559 * Checks the variable @x over the range @t1 to @t2 for its
560 * suitability as the dependent variable in an ordered probit
561 * analysis. The criterion used is that the variable has
562 * only non-negative integer values.
563 *
564 * Returns: 1 if the test succeeds, otherwise 0.
565 */
566
gretl_is_oprobit_ok(int t1,int t2,const double * x)567 int gretl_is_oprobit_ok (int t1, int t2, const double *x)
568 {
569 int t, n = 0, d = 1;
570
571 for (t=t1; t<=t2; t++) {
572 if (na(x[t])) {
573 continue;
574 }
575 n++;
576 if (x[t] != floor(x[t]) || x[t] < 0) {
577 d = 0;
578 break;
579 }
580 }
581
582 return (d > 0 && n > 0);
583 }
584
585 /**
586 * true_const:
587 * @v: index number of variable to test.
588 * @dset: dataset struct.
589 *
590 * Check whether variable Z[v] equals 1 over the sample
591 * range given in @dset (aside from any missing values).
592 *
593 * Returns: 1 if so, otherwise 0.
594 */
595
true_const(int v,const DATASET * dset)596 int true_const (int v, const DATASET *dset)
597 {
598 if (v < 0 || v >= dset->v) {
599 return 0;
600 }
601
602 return gretl_isunits(dset->t1, dset->t2, dset->Z[v]);
603 }
604
605 /**
606 * gretl_compare_doubles:
607 * @a: pointer to first element to compare.
608 * @b: pointer to second element to compare.
609 *
610 * Comparison function for use with qsort. Sorts doubles in
611 * ascending order.
612 *
613 * Returns: appropriate value for qsort.
614 */
615
gretl_compare_doubles(const void * a,const void * b)616 int gretl_compare_doubles (const void *a, const void *b)
617 {
618 const double *da = (const double *) a;
619 const double *db = (const double *) b;
620
621 if (isnan(*da) || isnan(*db)) {
622 if (!isnan(*da)) {
623 return -1;
624 } else if (!isnan(*db)) {
625 return 1;
626 } else {
627 return 0;
628 }
629 } else {
630 return (*da > *db) - (*da < *db);
631 }
632 }
633
634 /**
635 * gretl_inverse_compare_doubles:
636 * @a: pointer to first element to compare.
637 * @b: pointer to second element to compare.
638 *
639 * Comparison function for use with qsort. Sorts doubles in
640 * descending order.
641 *
642 * Returns: appropriate value for qsort.
643 */
644
gretl_inverse_compare_doubles(const void * a,const void * b)645 int gretl_inverse_compare_doubles (const void *a, const void *b)
646 {
647 const double *da = (const double *) a;
648 const double *db = (const double *) b;
649
650 if (isnan(*da) || isnan(*db)) {
651 if (!isnan(*da)) {
652 return 1;
653 } else if (!isnan(*db)) {
654 return -1;
655 } else {
656 return 0;
657 }
658 } else {
659 return (*da < *db) - (*da > *db);
660 }
661 }
662
663 /**
664 * gretl_compare_ints:
665 * @a: pointer to first element to compare.
666 * @b: pointer to second element to compare.
667 *
668 * Comparison function for use with qsort. Sorts integers in
669 * ascending order.
670 *
671 * Returns: appropriate value for qsort.
672 */
673
gretl_compare_ints(const void * a,const void * b)674 int gretl_compare_ints (const void *a, const void *b)
675 {
676 const int *ia = (const int *) a;
677 const int *ib = (const int *) b;
678
679 return *ia - *ib;
680 }
681
682 /**
683 * count_distinct_values:
684 * @x: sorted array of doubles.
685 * @n: number of elements in array.
686 *
687 * Returns: the number of distinct values in array @x,
688 * provided that @x is already sorted.
689 */
690
count_distinct_values(const double * x,int n)691 int count_distinct_values (const double *x, int n)
692 {
693 int i, c = 1;
694
695 for (i=1; i<n; i++) {
696 if (x[i] != x[i-1]) {
697 c++;
698 }
699 }
700
701 return c;
702 }
703
704 /**
705 * count_distinct_int_values:
706 * @x: sorted array of ints.
707 * @n: number of elements in array.
708 *
709 * Returns: the number of distinct values in array @x,
710 * provided that @x is already sorted.
711 */
712
count_distinct_int_values(const int * x,int n)713 int count_distinct_int_values (const int *x, int n)
714 {
715 int i, c = 1;
716
717 for (i=1; i<n; i++) {
718 if (x[i] != x[i-1]) c++;
719 }
720
721 return c;
722 }
723
724 /**
725 * rearrange_id_array:
726 * @x: sorted array of doubles.
727 * @m: number of distinct values in array.
728 * @n: number of elements in array.
729 *
730 * Rearranges the sorted array @x such that the first @m
731 * elements contain the @m distinct values in sorted order.
732 *
733 * Returns: 0 on success, 1 on error (in case @m is greater
734 * than @n).
735 */
736
rearrange_id_array(double * x,int m,int n)737 int rearrange_id_array (double *x, int m, int n)
738 {
739 int i, k = 1;
740
741 if (m >= n || m == 1) {
742 return 1;
743 }
744
745 for (i=1; i<n && k<m; i++) {
746 if (x[i] != x[i-1]) {
747 x[k++] = x[i];
748 }
749 }
750
751 return 0;
752 }
753
754 /**
755 * printlist:
756 * @list: array of integers.
757 * @msg: message to print along with @list (or NULL).
758 *
759 * Prints to stderr the given @list of integers along with a message.
760 */
761
printlist(const int * list,const char * msg)762 void printlist (const int *list, const char *msg)
763 {
764 int i;
765
766 if (msg) {
767 fprintf(stderr, "%s:\n", msg);
768 } else {
769 fprintf(stderr, "list: ");
770 }
771
772 if (list == NULL) {
773 fputs("list is NULL", stderr);
774 } else {
775 fprintf(stderr, "%d : ", list[0]);
776 for (i=1; i<=list[0]; i++) {
777 if (list[i] == LISTSEP) {
778 fputs("; ", stderr);
779 } else {
780 fprintf(stderr, "%d ", list[i]);
781 }
782 }
783 }
784
785 fputc('\n', stderr);
786 }
787
788 /**
789 * gretl_calculate_criteria:
790 * @ess: error sum of squares.
791 * @n: number of observations.
792 * @k: number of parameters estimated.
793 * @ll: pointer to recieve loglikelihood.
794 * @aic: pointer to recieve Akaike criterion.
795 * @bic: pointer to recieve Schwarz Bayesian criterion.
796 * @hqc: pointer to recieve Hannan-Quinn criterion.
797 *
798 * Calculates model selection criteria based on @ess, @nobs and
799 * @k, for a model estimated via least squares.
800 *
801 * Returns: 0 on success, non-zero on error.
802 */
803
gretl_calculate_criteria(double ess,int n,int k,double * ll,double * aic,double * bic,double * hqc)804 int gretl_calculate_criteria (double ess, int n, int k,
805 double *ll, double *aic, double *bic,
806 double *hqc)
807 {
808 double lnl, c[3];
809 int err = 0;
810
811 if (na(ess) || ess <= 0.0 || n <= k) {
812 err = 1;
813 } else {
814 errno = 0;
815 lnl = -.5 * n * log(ess);
816 if (errno == EDOM || errno == ERANGE) {
817 err = 1;
818 } else {
819 lnl += -.5 * n * (1 + LN_2_PI - log(n));
820 c[0] = -2.0 * lnl + 2 * k;
821 c[1] = -2.0 * lnl + k * log(n);
822 c[2] = -2.0 * lnl + 2 * k * log(log(n));
823 }
824 }
825
826 if (err) {
827 *ll = *aic = *bic = *hqc = NADBL;
828 } else {
829 *ll = lnl;
830 *aic = c[0];
831 *bic = c[1];
832 *hqc = c[2];
833 }
834
835 return err;
836 }
837
838 /**
839 * ls_criteria:
840 * @pmod: pointer to gretl model structure.
841 *
842 * Fills out the model selection criteria members of @pmod, using
843 * gretl_calculate_criteria().
844 *
845 * Returns: 0 on success, non-zero on error.
846 */
847
ls_criteria(MODEL * pmod)848 int ls_criteria (MODEL *pmod)
849 {
850 double ll, aic, bic, hqc;
851 int err;
852
853 err = gretl_calculate_criteria(pmod->ess, pmod->nobs, pmod->ncoeff,
854 &ll, &aic, &bic, &hqc);
855
856 pmod->lnL = ll;
857 pmod->criterion[C_AIC] = aic;
858 pmod->criterion[C_BIC] = bic;
859 pmod->criterion[C_HQC] = hqc;
860
861 return err;
862 }
863
864 static char *
real_format_obs(char * obs,int maj,int min,int pd,char sep)865 real_format_obs (char *obs, int maj, int min, int pd, char sep)
866 {
867 if (pd >= 10) {
868 int pdp = pd / 10, minlen = 2;
869 char fmt[18];
870
871 while ((pdp = pdp / 10)) minlen++;
872 sprintf(fmt, "%%d%c%%0%dd", sep, minlen);
873 sprintf(obs, fmt, maj, min);
874 } else {
875 sprintf(obs, "%d%c%d", maj, sep, min);
876 }
877
878 return obs;
879 }
880
881 /**
882 * format_obs:
883 * @obs: target string (should be of length #OBSLEN).
884 * @maj: major period (e.g. year).
885 * @min: minor period (e.g. quarter, month).
886 * @pd: data frequency.
887 *
888 * Prints to @obs the gretl-type date string representing
889 * the observation given by @maj, @min and @pd.
890 *
891 * Returns: @obs.
892 */
893
format_obs(char * obs,int maj,int min,int pd)894 char *format_obs (char *obs, int maj, int min, int pd)
895 {
896 return real_format_obs(obs, maj, min, pd, ':');
897 }
898
get_stobs_maj_min(char * stobs,int structure,int * maj,int * min)899 static int get_stobs_maj_min (char *stobs, int structure,
900 int *maj, int *min)
901 {
902 int dotc = 0;
903 char *p = stobs;
904 int err = 0;
905
906 while (*p) {
907 if (*p == ':') {
908 *p = '.';
909 dotc++;
910 } else if (*p == '.') {
911 dotc++;
912 } else if (!isdigit((unsigned char) *p)) {
913 err = 1;
914 break;
915 }
916 p++;
917 }
918
919 if (!err) {
920 if (dotc > 1 || *stobs == '.' ||
921 stobs[strlen(stobs) - 1] == '.') {
922 err = 1;
923 }
924 }
925
926 if (!err) {
927 if (dotc > 0) {
928 sscanf(stobs, "%d.%d", maj, min);
929 if (*maj <= 0 || *min <= 0) {
930 err = 1;
931 }
932 } else {
933 sscanf(stobs, "%d", maj);
934 if (*maj <= 0 && structure != SPECIAL_TIME_SERIES) {
935 err = 1;
936 }
937 }
938 }
939
940 return err;
941 }
942
943 static int
catch_setobs_errors(const char * stobs,int pd,int min,int structure)944 catch_setobs_errors (const char *stobs, int pd, int min, int structure)
945 {
946 int panel = structure == STACKED_TIME_SERIES ||
947 structure == STACKED_CROSS_SECTION;
948 int err = 0;
949
950 if (pd == 1) {
951 if (min > 0) {
952 gretl_errmsg_set(_("no ':' allowed in starting obs with "
953 "frequency 1"));
954 err = 1;
955 } else if (panel) {
956 gretl_errmsg_set(_("panel data must have frequency > 1"));
957 err = 1;
958 }
959 } else {
960 if (min == 0) {
961 gretl_errmsg_set(_("starting obs must contain a ':' with "
962 "frequency > 1"));
963 err = 1;
964 } else if (min > pd) {
965 gretl_errmsg_sprintf(_("starting obs '%s' is incompatible with frequency"),
966 stobs);
967 err = 1;
968 } else if (structure == CROSS_SECTION) {
969 gretl_errmsg_set(_("cross-sectional data: frequency must be 1"));
970 err = 1;
971 }
972 }
973
974 return err;
975 }
976
invalid_stobs(const char * s)977 static int invalid_stobs (const char *s)
978 {
979 gretl_errmsg_sprintf(_("starting obs '%s' is invalid"), s);
980 return E_DATA;
981 }
982
maybe_fix_daily_start(guint32 * ed,int pd)983 static void maybe_fix_daily_start (guint32 *ed, int pd)
984 {
985 int wday = weekday_from_epoch_day(*ed);
986 int fix = 0;
987
988 if (wday == 0) {
989 /* 5- or 6-day data: sunday not valid */
990 fix = 1;
991 } else if (wday == 6 && pd == 5) {
992 /* 5-day data: saturday not valid */
993 fix = 2;
994 }
995
996 if (fix) {
997 char *fixed, *msg;
998
999 *ed += fix;
1000 fixed = ymd_extended_from_epoch_day(*ed, 0, NULL);
1001 msg = gretl_strdup_printf("the starting date was corrected to Monday %s",
1002 fixed);
1003 gretl_warnmsg_set(msg);
1004 free(msg);
1005 free(fixed);
1006 }
1007 }
1008
1009 #define likely_calendar_obs_string(s) (strchr(s, '-') || strchr(s, '/'))
1010
1011 #define recognized_ts_frequency(f) (f == 4 || f == 12 || f == 24)
1012
process_starting_obs(const char * stobs_in,int pd,int * pstructure,double * psd0,guint32 * ped0,gretlopt opt)1013 static int process_starting_obs (const char *stobs_in, int pd,
1014 int *pstructure, double *psd0,
1015 guint32 *ped0, gretlopt opt)
1016 {
1017 char stobs[OBSLEN];
1018 int structure = *pstructure;
1019 double sd0 = 0.0;
1020 int maybe_tseries = 1;
1021 int dated = 0;
1022 int err = 0;
1023
1024 if (structure == CROSS_SECTION ||
1025 structure == STACKED_TIME_SERIES ||
1026 structure == STACKED_CROSS_SECTION) {
1027 maybe_tseries = 0;
1028 }
1029
1030 *stobs = '\0';
1031 strncat(stobs, stobs_in, OBSLEN - 1);
1032
1033 /* truncate stobs if not a calendar date */
1034
1035 if (likely_calendar_obs_string(stobs)) {
1036 if (maybe_tseries) {
1037 dated = 1;
1038 } else {
1039 return invalid_stobs(stobs);
1040 }
1041 } else {
1042 stobs[8] = '\0';
1043 }
1044
1045 if (dated) {
1046 if (pd == 5 || pd == 6 || pd == 7 || pd == 52) {
1047 /* calendar-dated data, daily or weekly */
1048 guint32 ed0 = get_epoch_day(stobs);
1049
1050 if (ed0 <= 0) {
1051 return invalid_stobs(stobs);
1052 } else {
1053 if (pd < 7) {
1054 maybe_fix_daily_start(&ed0, pd);
1055 }
1056 sd0 = ed0;
1057 *ped0 = ed0;
1058 structure = TIME_SERIES;
1059 }
1060 } else {
1061 return invalid_stobs(stobs);
1062 }
1063 } else if (structure == TIME_SERIES && pd == 10) {
1064 /* decennial data */
1065 sd0 = (double) atoi(stobs);
1066 } else {
1067 int maj = 0, min = 0;
1068
1069 if (get_stobs_maj_min(stobs, structure, &maj, &min)) {
1070 return invalid_stobs(stobs);
1071 }
1072
1073 if ((pd == 5 || pd == 6 || pd == 7 || pd == 52) &&
1074 min == 0 && maybe_tseries) {
1075 /* catch undated daily or weekly data */
1076 structure = TIME_SERIES;
1077 } else {
1078 if (catch_setobs_errors(stobs, pd, min, structure)) {
1079 return E_DATA;
1080 } else if (pd == 1) {
1081 sprintf(stobs, "%d", maj);
1082 if (structure == STRUCTURE_UNKNOWN) {
1083 if (maj > 1) {
1084 structure = TIME_SERIES; /* annual? */
1085 } else {
1086 structure = CROSS_SECTION;
1087 }
1088 }
1089 } else {
1090 if (structure == TIME_SERIES && min > 0 &&
1091 !recognized_ts_frequency(pd)) {
1092 if (opt & OPT_I) {
1093 return invalid_stobs(stobs);
1094 } else {
1095 structure = SPECIAL_TIME_SERIES;
1096 }
1097 }
1098 real_format_obs(stobs, maj, min, pd, '.');
1099 if (structure == STRUCTURE_UNKNOWN &&
1100 recognized_ts_frequency(pd)) {
1101 structure = TIME_SERIES;
1102 }
1103 }
1104 }
1105
1106 /* for non-calendar data */
1107 sd0 = dot_atof(stobs);
1108 }
1109
1110 if (!err) {
1111 *pstructure = structure;
1112 *psd0 = sd0;
1113 }
1114
1115 return err;
1116 }
1117
1118 /**
1119 * set_obs:
1120 * @parm1: first parameter.
1121 * @parm2: second parameter.
1122 * @dset: dataset struct.
1123 * @opt: %OPT_S for stacked time-series, %OPT_C for stacked cross-section,
1124 * %OPT_T for time series, %OPT_X for cross section, %OPT_P to set
1125 * panel structure via two variables representing unit and period
1126 * respectively. For data already set as panel, %OPT_G to set panel
1127 * group names or %OPT_I to set panel time-dimension information.
1128 *
1129 * Set the frequency and initial observation for a dataset, or
1130 * in the case of a panel dataset, extra group or time information.
1131 *
1132 * Returns: 0 on successful completion, non-zero on error.
1133 */
1134
set_obs(const char * parm1,const char * parm2,DATASET * dset,gretlopt opt)1135 int set_obs (const char *parm1, const char *parm2,
1136 DATASET *dset, gretlopt opt)
1137 {
1138 const char *stobs = NULL;
1139 int structure = STRUCTURE_UNKNOWN;
1140 double sd0 = dset->sd0;
1141 guint32 ed0 = 0;
1142 int pd, panel = 0;
1143 int err = 0;
1144
1145 if (dset == NULL) {
1146 return E_NODATA;
1147 }
1148
1149 if ((opt & (OPT_R | OPT_P)) && dset->Z == NULL) {
1150 return E_NODATA;
1151 }
1152
1153 if ((opt & (OPT_G | OPT_I)) && !dataset_is_panel(dset)) {
1154 return E_DATA;
1155 }
1156
1157 /* prevent substantive reorganization of the dataset
1158 within a function */
1159 if ((opt & (OPT_P | OPT_C)) && gretl_function_depth() > 0) {
1160 gretl_errmsg_set("You cannot do this within a function");
1161 return E_DATA;
1162 }
1163
1164 gretl_error_clear();
1165
1166 if (opt & OPT_R) {
1167 /* restructure panel: "hidden" option */
1168 return switch_panel_orientation(dset);
1169 }
1170
1171 if (opt == OPT_NONE && parm1 == NULL && parm2 == NULL) {
1172 /* we'll just print current obs info */
1173 return 0;
1174 }
1175
1176 if (parm1 == NULL || (!(opt & OPT_G) && parm2 == NULL)) {
1177 /* all cases need at least one param, most need two */
1178 return E_ARGS;
1179 }
1180
1181 if (opt & OPT_P) {
1182 return set_panel_structure_from_varnames(parm1, parm2, dset);
1183 } else if (opt & OPT_G) {
1184 /* --panel-groups */
1185 return set_panel_group_strings(parm1, parm2, dset);
1186 }
1187
1188 /* now we get down to business */
1189
1190 pd = gretl_int_from_string(parm1, &err);
1191 if (!err && pd < 1) {
1192 gretl_errmsg_sprintf(_("frequency (%d) does not seem to make sense"), pd);
1193 return E_DATA;
1194 }
1195
1196 stobs = parm2;
1197
1198 /* if an explicit structure option was passed in, respect it */
1199 if (opt == OPT_X) {
1200 structure = CROSS_SECTION;
1201 } else if (opt == OPT_T) {
1202 structure = TIME_SERIES;
1203 } else if (opt == OPT_S) {
1204 structure = STACKED_TIME_SERIES;
1205 panel = 1;
1206 } else if (opt == OPT_C) {
1207 structure = STACKED_CROSS_SECTION;
1208 panel = 1;
1209 } else if (opt == OPT_N) {
1210 structure = SPECIAL_TIME_SERIES;
1211 } else if (opt == OPT_I) {
1212 /* --panel-time */
1213 structure = TIME_SERIES;
1214 }
1215
1216 if (panel && dset->n > 0 && pd > dset->n) {
1217 gretl_errmsg_sprintf(_("frequency (%d) does not seem to make sense"), pd);
1218 return 1;
1219 }
1220
1221 err = process_starting_obs(stobs, pd, &structure, &sd0, &ed0, opt);
1222
1223 if (err) {
1224 return err;
1225 }
1226
1227 if (opt == OPT_I) {
1228 dset->panel_pd = pd;
1229 dset->panel_sd0 = sd0;
1230 return 0;
1231 }
1232
1233 if (panel && dset->n % pd != 0) {
1234 int sts = structure == STACKED_TIME_SERIES;
1235
1236 gretl_errmsg_sprintf(_("Panel datasets must be balanced.\n"
1237 "The number of observations (%d) is not a multiple\n"
1238 "of the number of %s (%d)."),
1239 dset->n, sts ? _("periods") : _("units"), pd);
1240 return E_DATA;
1241 }
1242
1243 if (ed0 > 0) {
1244 /* replace any existing markers with date strings */
1245 dataset_destroy_obs_markers(dset);
1246 } else if (structure == TIME_SERIES && (pd == 1 || pd == 4 || pd == 12)) {
1247 /* force use of regular time-series obs labels */
1248 dataset_destroy_obs_markers(dset);
1249 }
1250
1251 dset->pd = pd;
1252 dset->structure = structure;
1253 dset->sd0 = sd0;
1254
1255 if (ed0 > 0) {
1256 calendar_date_string(dset->stobs, 0, dset);
1257 calendar_date_string(dset->endobs, dset->n - 1, dset);
1258 } else {
1259 ntolabel(dset->stobs, 0, dset);
1260 ntolabel(dset->endobs, dset->n - 1, dset);
1261 }
1262
1263 /* pre-process stacked cross-sectional panels: put into canonical
1264 stacked time series form
1265 */
1266 if (dset->structure == STACKED_CROSS_SECTION) {
1267 err = switch_panel_orientation(dset);
1268 }
1269
1270 #if 0
1271 fprintf(stderr, "setobs: pd=%d, stobs=%s, sd0=%g, markers=%d, S=%p\n",
1272 dset->pd, dset->stobs, dset->sd0, dset->markers, (void *) dset->S);
1273 #endif
1274
1275 return err;
1276 }
1277
1278 /**
1279 * simple_set_obs:
1280 * @dset: pointer to data information struct.
1281 * @pd: data frequency.
1282 * @stobs: string representation of starting observation.
1283 * @opt: options flags.
1284 *
1285 * See the "setobs" command.
1286 *
1287 * Returns: 0 on success, non-zero code on error.
1288 */
1289
simple_set_obs(DATASET * dset,int pd,const char * stobs,gretlopt opt)1290 int simple_set_obs (DATASET *dset, int pd, const char *stobs,
1291 gretlopt opt)
1292 {
1293 int structure = STRUCTURE_UNKNOWN;
1294 double sd0 = dset->sd0;
1295 guint32 ed0 = 0;
1296 int panel = 0;
1297 int err = 0;
1298
1299 if (opt == OPT_X) {
1300 structure = CROSS_SECTION;
1301 } else if (opt == OPT_T) {
1302 structure = TIME_SERIES;
1303 } else if (opt == OPT_S) {
1304 structure = STACKED_TIME_SERIES;
1305 panel = 1;
1306 } else if (opt == OPT_C) {
1307 structure = STACKED_CROSS_SECTION;
1308 panel = 1;
1309 } else if (opt == OPT_N) {
1310 structure = SPECIAL_TIME_SERIES;
1311 } else if (opt == OPT_I) {
1312 /* --panel-time */
1313 structure = TIME_SERIES;
1314 }
1315
1316 err = process_starting_obs(stobs, pd, &structure, &sd0, &ed0, OPT_NONE);
1317
1318 if (err) {
1319 return err;
1320 }
1321
1322 if (panel && dset->n % pd != 0) {
1323 int sts = structure == STACKED_TIME_SERIES;
1324
1325 gretl_errmsg_sprintf(_("Panel datasets must be balanced.\n"
1326 "The number of observations (%d) is not a multiple\n"
1327 "of the number of %s (%d)."),
1328 dset->n, sts ? _("periods") : _("units"), pd);
1329 return E_DATA;
1330 }
1331
1332 if (ed0 > 0) {
1333 /* replace any existing markers with date strings */
1334 dataset_destroy_obs_markers(dset);
1335 } else if (structure == TIME_SERIES && (pd == 1 || pd == 4 || pd == 12)) {
1336 /* force use of regular time-series obs labels */
1337 dataset_destroy_obs_markers(dset);
1338 }
1339
1340 dset->pd = pd;
1341 dset->structure = structure;
1342 dset->sd0 = sd0;
1343
1344 if (ed0 > 0) {
1345 calendar_date_string(dset->stobs, 0, dset);
1346 calendar_date_string(dset->endobs, dset->n - 1, dset);
1347 } else {
1348 ntolabel(dset->stobs, 0, dset);
1349 ntolabel(dset->endobs, dset->n - 1, dset);
1350 }
1351
1352 if (dset->structure == STACKED_CROSS_SECTION) {
1353 err = switch_panel_orientation(dset);
1354 }
1355
1356 return err;
1357 }
1358
1359 /**
1360 * gretl_double_from_string:
1361 * @s: string to examine.
1362 * @err: location to receive error code.
1363 *
1364 * If @s is a valid string representation of a double,
1365 * return its value, otherwise if @s is the name of a
1366 * scalar variable, return the value of that variable,
1367 * otherwise set the content of @err to a non-zero value.
1368 *
1369 * Returns: double value.
1370 */
1371
gretl_double_from_string(const char * s,int * err)1372 double gretl_double_from_string (const char *s, int *err)
1373 {
1374 char *test;
1375 double x;
1376
1377 if (s == NULL || *s == '\0') {
1378 *err = E_DATA;
1379 return NADBL;
1380 }
1381
1382 if (isalpha(*s)) {
1383 return get_scalar_value_by_name(s, err);
1384 }
1385
1386 gretl_push_c_numeric_locale();
1387 errno = 0;
1388 x = strtod(s, &test);
1389 gretl_pop_c_numeric_locale();
1390
1391 if (*test != '\0' || errno == ERANGE) {
1392 *err = E_DATA;
1393 errno = 0;
1394 return NADBL;
1395 }
1396
1397 return x;
1398 }
1399
1400 /**
1401 * gretl_int_from_string:
1402 * @s: string to examine.
1403 * @err: location to receive error code.
1404 *
1405 * If @s is a valid string representation of an integer,
1406 * return that integer, otherwise if @s is the name of a
1407 * scalar variable, return the value of that variable,
1408 * provided it can be converted to an integer, otherwise
1409 * set the content of @err to a non-zero value.
1410 *
1411 * Returns: integer value.
1412 */
1413
gretl_int_from_string(const char * s,int * err)1414 int gretl_int_from_string (const char *s, int *err)
1415 {
1416 char *test;
1417 int n = 0;
1418
1419 if (s == NULL || *s == '\0') {
1420 *err = E_DATA;
1421 return 0;
1422 }
1423
1424 if (isalpha(*s)) {
1425 double x = get_scalar_value_by_name(s, err);
1426
1427 if (!*err) {
1428 n = gretl_int_from_double(x, err);
1429 }
1430 return n;
1431 }
1432
1433 errno = 0;
1434 n = strtol(s, &test, 10);
1435
1436 if (*test != '\0' || errno == ERANGE) {
1437 *err = E_DATA;
1438 errno = 0;
1439 return 0;
1440 }
1441
1442 return n;
1443 }
1444
1445 /**
1446 * positive_int_from_string:
1447 * @s: string to examine.
1448 *
1449 * If @s is a valid string representation of a positive integer,
1450 * return that integer, otherwise return -1.
1451 *
1452 * Returns: integer value.
1453 */
1454
positive_int_from_string(const char * s)1455 int positive_int_from_string (const char *s)
1456 {
1457 int ret = -1;
1458
1459 if (s != NULL && *s != '\0') {
1460 char *test;
1461
1462 errno = 0;
1463
1464 ret = strtol(s, &test, 10);
1465 if (*test != '\0' || !strcmp(s, test) || errno == ERANGE) {
1466 ret = -1;
1467 }
1468 }
1469
1470 return ret;
1471 }
1472
letter_to_int(char c)1473 static int letter_to_int (char c)
1474 {
1475 const char *s = "abcdefghij";
1476 int i = 0;
1477
1478 while (*s) {
1479 if (c == *s) {
1480 return i;
1481 }
1482 s++;
1483 i++;
1484 }
1485
1486 return 0;
1487 }
1488
1489 /**
1490 * gretl_version_number:
1491 * @version: gretl program version in string form.
1492 *
1493 * Returns: the integer gretl version number.
1494 */
1495
gretl_version_number(const char * version)1496 int gretl_version_number (const char *version)
1497 {
1498 int vnum = 0;
1499
1500 if (atoi(version) >= 2015) {
1501 /* as in "2015d" and subsequent releases */
1502 int Y;
1503 char c;
1504
1505 sscanf(version, "%d%c", &Y, &c);
1506 vnum = 10 * Y + letter_to_int(c);
1507 } else {
1508 /* old style: "1.9.14" or whatever */
1509 int x, y, z;
1510
1511 sscanf(version, "%d.%d.%d", &x, &y, &z);
1512 vnum = 10000 * x + 100 * y + z;
1513 }
1514
1515 return vnum;
1516 }
1517
1518 /* Mapping from new-style version numbers, based on
1519 year, to old style based on maj.min.pl. This
1520 covers gretl 1.9.4 to 1.10.2.
1521 */
1522
old_style_gretl_version_number(int v)1523 static int old_style_gretl_version_number (int v)
1524 {
1525 const int vtrans[18][2] = {
1526 {10904, 20110},
1527 {10905, 20111},
1528 {10906, 20112},
1529 {10907, 20113},
1530 {10908, 20120},
1531 {10909, 20121},
1532 {10910, 20122},
1533 {10911, 20123},
1534 {10912, 20130},
1535 {10913, 20131},
1536 {10914, 20132},
1537 {10990, 20140},
1538 {10991, 20141},
1539 {10992, 20142},
1540 {11000, 20150},
1541 {11001, 20151},
1542 {11002, 20152},
1543 {11003, 20153}
1544 };
1545 int i;
1546
1547 for (i=0; i<17; i++) {
1548 if (v == vtrans[i][1] || v < vtrans[i+1][1]) {
1549 return vtrans[i][0];
1550 }
1551 }
1552
1553 /* default to 1.9.4 */
1554
1555 return vtrans[0][0];
1556 }
1557
1558 /**
1559 * gretl_version_string:
1560 * @targ: string into which to write (9 bytes minimum).
1561 * @vnum: integer program version number.
1562 *
1563 * Returns: the string representation of @vnum, in @targ.
1564 */
1565
gretl_version_string(char * targ,int vnum)1566 char *gretl_version_string (char *targ, int vnum)
1567 {
1568 if (vnum >= 20153) {
1569 /* "2105d" or higher */
1570 const char *s = "abcdefghij";
1571 int y, idx;
1572 char c;
1573
1574 y = vnum / 10;
1575 idx = vnum - 10 * y;
1576
1577 if (idx >= 0 && idx < 10) {
1578 c = s[idx];
1579 } else {
1580 c = 'a';
1581 }
1582
1583 sprintf(targ, "%d%c", y, c);
1584 } else {
1585 int x, y, z;
1586
1587 if (vnum > 20000) {
1588 /* translate back to old-style */
1589 vnum = old_style_gretl_version_number(vnum);
1590 }
1591
1592 x = vnum / 10000;
1593 y = (vnum - x * 10000) / 100;
1594 z = vnum % 100;
1595
1596 sprintf(targ, "%d.%d.%d", x, y, z);
1597 }
1598
1599 return targ;
1600 }
1601
1602 /**
1603 * varnum_from_string:
1604 * @str: string representation of an integer ID number.
1605 * @dset: dataset information.
1606 *
1607 * Returns: integer ID number, or -1 on failure.
1608 */
1609
varnum_from_string(const char * str,DATASET * dset)1610 int varnum_from_string (const char *str, DATASET *dset)
1611 {
1612 int v = positive_int_from_string(str);
1613
1614 if (v <= 0 || v >= dset->v) {
1615 v = -1;
1616 }
1617
1618 return v;
1619 }
1620
1621 /**
1622 * gretl_int_from_double:
1623 * @x: double-precision floating point value
1624 * @err: location to receive error code.
1625 *
1626 * Returns: the value of @x converted to an integer, if
1627 * possible. Otherwise returns -1 with @err set to a
1628 * non-zero value. Note that it is considered an
1629 * error if @x is "too far" from the nearest integer;
1630 * it must be "almost integral", with tolerance 0.001.
1631 */
1632
gretl_int_from_double(double x,int * err)1633 int gretl_int_from_double (double x, int *err)
1634 {
1635 int k = -1;
1636
1637 if (na(x) || fabs(x) > INT_MAX || fabs(x - nearbyint(x)) > 0.001) {
1638 *err = E_INVARG;
1639 } else {
1640 k = (int) lrint(x);
1641 }
1642
1643 return k;
1644 }
1645
1646 /**
1647 * gretl_unsigned_from_double:
1648 * @x: double-precision floating point value
1649 * @err: location to receive error code.
1650 *
1651 * Returns: the value of @x converted to an integer, if
1652 * possible. Otherwise returns -1 with @err set to a
1653 * non-zero value. Note that it is considered an
1654 * error if @x is "too far" from the nearest integer;
1655 * it must be "almost integral", with tolerance 0.001.
1656 */
1657
gretl_unsigned_from_double(double x,int * err)1658 guint32 gretl_unsigned_from_double (double x, int *err)
1659 {
1660 guint32 u = 0;
1661
1662 if (na(x) || x < 0 || fabs(x) > UINT_MAX) {
1663 *err = E_INVARG;
1664 } else {
1665 double f = floor(x);
1666 double c = ceil(x);
1667
1668 if (x - f < 1e-6) {
1669 u = f;
1670 } else if (c - x < 1e-6) {
1671 u = c;
1672 } else {
1673 *err = E_INVARG;
1674 }
1675 }
1676
1677 return u;
1678 }
1679
1680 /**
1681 * gretl_type_from_name:
1682 * @s: the name to check.
1683 * @dset: dataset information.
1684 *
1685 * Returns: non-zero if @s is the name of an existing
1686 * series, matrix, scalar, list, string or bundle,
1687 * otherwise 0.
1688 */
1689
gretl_type_from_name(const char * s,const DATASET * dset)1690 GretlType gretl_type_from_name (const char *s, const DATASET *dset)
1691 {
1692 if (dset != NULL && gretl_is_series(s, dset)) {
1693 return GRETL_TYPE_SERIES;
1694 } else {
1695 return user_var_get_type_by_name(s);
1696 }
1697 }
1698
1699 /**
1700 * copyvec:
1701 * @src: array of doubles.
1702 * @n: number of elements to copy.
1703 *
1704 * Returns: an allocated copy of the first @n elements of
1705 * array @src, or %NULL on failure.
1706 */
1707
copyvec(const double * src,int n)1708 double *copyvec (const double *src, int n)
1709 {
1710 double *targ = NULL;
1711 size_t sz = n * sizeof *targ;
1712
1713 if (n > 0 && src != NULL) {
1714 targ = malloc(sz);
1715 }
1716
1717 if (targ != NULL) {
1718 memcpy(targ, src, sz);
1719 }
1720
1721 return targ;
1722 }
1723
1724 /**
1725 * doubles_array_free:
1726 * @X: 2-dimensional array of doubles.
1727 * @m: number of sub-arrays.
1728 *
1729 * Frees a 2-dimensional array of doubles, first freeing
1730 * each sub-array.
1731 */
1732
doubles_array_free(double ** X,int m)1733 void doubles_array_free (double **X, int m)
1734 {
1735 if (X != NULL) {
1736 int i;
1737
1738 for (i=0; i<m; i++) {
1739 free(X[i]);
1740 }
1741 free(X);
1742 }
1743 }
1744
1745 /**
1746 * doubles_array_new:
1747 * @m: number of sub-arrays.
1748 * @n: length of each sub-array.
1749 *
1750 * Allocates a 2-dimensional array of doubles, that is,
1751 * @m arrays each containing @n elements. If @n is
1752 * zero the sub-arrays are just set to %NULL.
1753 *
1754 * Returns: the allocated array, or %NULL on failure.
1755 */
1756
doubles_array_new(int m,int n)1757 double **doubles_array_new (int m, int n)
1758 {
1759 double **X;
1760 int i;
1761
1762 if (m == 0) {
1763 return NULL;
1764 }
1765
1766 X = malloc(m * sizeof *X);
1767 if (X == NULL) {
1768 return X;
1769 }
1770
1771 for (i=0; i<m; i++) {
1772 X[i] = NULL;
1773 }
1774
1775 if (n > 0) {
1776 for (i=0; i<m; i++) {
1777 X[i] = malloc(n * sizeof **X);
1778 if (X[i] == NULL) {
1779 doubles_array_free(X, m);
1780 X = NULL;
1781 break;
1782 }
1783 }
1784 }
1785
1786 return X;
1787 }
1788
1789 /**
1790 * doubles_array_new0:
1791 * @m: number of sub-arrays.
1792 * @n: length of each sub-array.
1793 *
1794 * Works just as doubles_array_new(), except that on
1795 * successful allocation all values in the arrays are
1796 * set to zero.
1797 *
1798 * Returns: the allocated array, or %NULL on failure.
1799 */
1800
doubles_array_new0(int m,int n)1801 double **doubles_array_new0 (int m, int n)
1802 {
1803 double **X = doubles_array_new(m, n);
1804
1805 if (X != NULL && n > 0) {
1806 int i, j;
1807
1808 for (i=0; i<m; i++) {
1809 for (j=0; j<n; j++) {
1810 X[i][j] = 0.0;
1811 }
1812 }
1813 }
1814
1815 return X;
1816 }
1817
1818 /**
1819 * doubles_array_adjust_length:
1820 * @X: original two-dimensional array.
1821 * @m: number of sub-arrays in @X.
1822 * @new_n: new length of each sub-array.
1823 *
1824 * For each of the @m sub-arrays in @X, reallocate to
1825 * a length of @new_n.
1826 *
1827 * Returns: 0 on success, non-zero on error.
1828 */
1829
doubles_array_adjust_length(double ** X,int m,int new_n)1830 int doubles_array_adjust_length (double **X, int m, int new_n)
1831 {
1832 int err = 0;
1833
1834 if (X == NULL || m == 0) {
1835 ; /* no-op */
1836 } else {
1837 double *xi;
1838 int i;
1839
1840 for (i=0; i<m && !err; i++) {
1841 if (new_n == 0) {
1842 free(X[i]);
1843 X[i] = NULL;
1844 } else {
1845 xi = realloc(X[i], new_n * sizeof *xi);
1846 if (xi == NULL) {
1847 err = E_ALLOC;
1848 } else {
1849 X[i] = xi;
1850 }
1851 }
1852 }
1853 }
1854
1855 return err;
1856 }
1857
1858 /**
1859 * data_array_from_model:
1860 * @pmod: reference model.
1861 * @Z: main data array.
1862 * @missv: should equal 1 if there are missing values to be
1863 * skipped, else 0.
1864 *
1865 * Constructs a dataset containing all the variables referenced in
1866 * @pmod. The arrays start at the correct sample offset for @pmod,
1867 * and are contiguous. If @missvals equals 0, this is done by creating
1868 * a set of pointers into the main dataset, but if there are missing
1869 * values to be handled, the sub-arrays are newly allocated and purged
1870 * of NAs.
1871 *
1872 * Returns: two-dimensional array, or %NULL on failure.
1873 */
1874
data_array_from_model(const MODEL * pmod,double ** Z,int missv)1875 double **data_array_from_model (const MODEL *pmod, double **Z, int missv)
1876 {
1877 double **X;
1878 int nv = pmod->list[0];
1879 int offset = pmod->t1;
1880 int v, i;
1881
1882 if (missv) {
1883 X = doubles_array_new(nv, pmod->nobs);
1884 } else {
1885 X = malloc(nv * sizeof *X);
1886 }
1887
1888 if (X == NULL) return NULL;
1889
1890 if (missv) {
1891 int t, s;
1892
1893 for (t=0; t<pmod->nobs; t++) {
1894 X[0][t] = 1.0;
1895 }
1896
1897 for (i=1; i<nv; i++) {
1898 v = (i == 1)? pmod->list[1] : pmod->list[i + 1];
1899 s = 0;
1900 for (t=pmod->t1; t<=pmod->t2; t++) {
1901 if (!na(pmod->uhat[t])) {
1902 X[i][s++] = Z[v][t];
1903 }
1904 }
1905 }
1906 } else {
1907 /* constant in slot 0 */
1908 X[0] = Z[0] + offset;
1909
1910 /* dependent var in slot 1 */
1911 X[1] = Z[pmod->list[1]] + offset;
1912
1913 /* independent vars in slots 2, 3, ... */
1914 for (i=2; i<nv; i++) {
1915 v = pmod->list[i + 1];
1916 X[i] = Z[v] + offset;
1917 }
1918 }
1919
1920 return X;
1921 }
1922
1923 #ifndef WIN32
1924
gretl_spawn(char * cmdline)1925 int gretl_spawn (char *cmdline)
1926 {
1927 GError *error = NULL;
1928 gchar *errout = NULL;
1929 gchar *sout = NULL;
1930 int ok, status;
1931 int ret = 0;
1932
1933 gretl_error_clear();
1934
1935 ok = g_spawn_command_line_sync(cmdline,
1936 &sout, /* standard output */
1937 &errout, /* standard error */
1938 &status, /* exit status */
1939 &error);
1940
1941 if (!ok) {
1942 gretl_errmsg_set(error->message);
1943 fprintf(stderr, "gretl_spawn: '%s'\n", error->message);
1944 g_error_free(error);
1945 if (errout != NULL) {
1946 fprintf(stderr, " stderr = '%s'\n", errout);
1947 }
1948 ret = 1;
1949 } else if (status != 0) {
1950 if (errout != NULL && *errout) {
1951 gretl_errmsg_set(errout);
1952 fprintf(stderr, "gretl_spawn: status = %d: '%s'\n", status, errout);
1953 } else if (sout != NULL && *sout) {
1954 gretl_errmsg_set(sout);
1955 fprintf(stderr, "gretl_spawn: status = %d: '%s'\n", status, sout);
1956 } else {
1957 gretl_errmsg_set(_("Command failed"));
1958 fprintf(stderr, "gretl_spawn: status = %d\n", status);
1959 }
1960 ret = 1;
1961 }
1962
1963 if (errout != NULL) g_free(errout);
1964 if (sout != NULL) g_free(sout);
1965
1966 if (ret) {
1967 fprintf(stderr, "Failed command: '%s'\n", cmdline);
1968 }
1969
1970 return ret;
1971 }
1972
1973 #endif /* !WIN32 */
1974
1975 /* file copying */
1976
gretl_copy_file(const char * src,const char * dest)1977 int gretl_copy_file (const char *src, const char *dest)
1978 {
1979 FILE *srcfd, *destfd;
1980 char buf[8192];
1981 size_t n;
1982
1983 if (!strcmp(src, dest)) {
1984 /* @src and @dest are the same: no-op */
1985 return 0;
1986 }
1987
1988 if ((srcfd = gretl_fopen(src, "rb")) == NULL) {
1989 gretl_errmsg_sprintf(_("Couldn't open %s"), src);
1990 return E_FOPEN;
1991 }
1992
1993 if ((destfd = gretl_fopen(dest, "wb")) == NULL) {
1994 gretl_errmsg_sprintf(_("Couldn't write to %s"), dest);
1995 fclose(srcfd);
1996 return E_FOPEN;
1997 }
1998
1999 while ((n = fread(buf, 1, sizeof buf, srcfd)) > 0) {
2000 fwrite(buf, 1, n, destfd);
2001 }
2002
2003 fclose(srcfd);
2004 fclose(destfd);
2005
2006 return 0;
2007 }
2008
maybe_unload_function_package(const char * s,PRN * prn)2009 static int maybe_unload_function_package (const char *s,
2010 PRN *prn)
2011 {
2012 char *p, pkgname[FN_NAMELEN+4];
2013 const char *path;
2014 fnpkg *pkg;
2015 int done = 0;
2016
2017 *pkgname = '\0';
2018 strncat(pkgname, s, FN_NAMELEN+3);
2019 p = strrchr(pkgname, '.');
2020 *p = '\0';
2021
2022 pkg = get_function_package_by_name(pkgname);
2023
2024 if (pkg != NULL) {
2025 path = get_function_package_path_by_name(pkgname);
2026 if (path != NULL) {
2027 function_package_unload_full_by_filename(path);
2028 done = 1;
2029 }
2030 }
2031
2032 if (done) {
2033 pprintf(prn, "Unloaded package %s\n", pkgname);
2034 } else {
2035 pprintf(prn, "Package %s was not loaded\n", pkgname);
2036 }
2037
2038 return 0;
2039 }
2040
2041 /* note: this is not used for series */
2042
gretl_delete_var_by_name(const char * s,PRN * prn)2043 int gretl_delete_var_by_name (const char *s, PRN *prn)
2044 {
2045 int err = 0;
2046
2047 if (s == NULL || *s == '\0') {
2048 return E_PARSE;
2049 }
2050
2051 if (object_is_function_arg(s)) {
2052 gretl_errmsg_sprintf(_("The variable %s is read-only"), s);
2053 return E_DATA;
2054 }
2055
2056 if (has_suffix(s, ".gfn")) {
2057 err = maybe_unload_function_package(s, prn);
2058 } else if (gretl_is_user_var(s)) {
2059 err = user_var_delete_by_name(s, prn);
2060 } else {
2061 /* try for a bundle member? */
2062 gchar *genstr = g_strdup_printf("%s=null", s);
2063
2064 err = generate(genstr, NULL, GRETL_TYPE_ANY, OPT_P, prn);
2065 g_free(genstr);
2066 }
2067
2068 return err;
2069 }
2070
2071 /* apparatus to support gretl's stopwatch */
2072
2073 #ifdef WANT_XTIMER
2074
2075 /* Both OpenMP and MPI offer timers that return the
2076 time as a double (seconds and fractions thereof).
2077 If we may be using either of these we'll define
2078 the "xtimers".
2079 */
2080
2081 struct xtimer {
2082 int level;
2083 double t0;
2084 };
2085
2086 static GPtrArray *xtimers;
2087
get_xtime(void)2088 static double get_xtime (void)
2089 {
2090 #if defined(HAVE_MPI)
2091 if (gretl_mpi_initialized()) {
2092 return gretl_mpi_time();
2093 }
2094 #endif
2095 #if defined(WIN32)
2096 return win32_get_time();
2097 #elif defined(_OPENMP)
2098 return omp_get_wtime();
2099 #endif
2100 return 0; /* not reached */
2101 }
2102
new_xtimer(int level)2103 static struct xtimer *new_xtimer (int level)
2104 {
2105 struct xtimer *xt = g_malloc(sizeof *xt);
2106
2107 xt->level = level;
2108 xt->t0 = get_xtime();
2109 g_ptr_array_insert(xtimers, xtimers->len, xt);
2110 return xt;
2111 }
2112
xtimer_init(void)2113 static void xtimer_init (void)
2114 {
2115 if (xtimers == NULL) {
2116 xtimers = g_ptr_array_sized_new(1);
2117 new_xtimer(gretl_function_depth());
2118 }
2119 }
2120
get_xtimer(void)2121 static struct xtimer *get_xtimer (void)
2122 {
2123 struct xtimer *xt;
2124 int i, lev = gretl_function_depth();
2125
2126 if (xtimers == NULL) {
2127 xtimers = g_ptr_array_sized_new(1);
2128 return new_xtimer(lev);
2129 }
2130
2131 for (i=0; i<xtimers->len; i++) {
2132 xt = g_ptr_array_index(xtimers, i);
2133 if (xt->level == lev) {
2134 return xt;
2135 }
2136 }
2137
2138 return new_xtimer(lev);
2139 }
2140
xtimer_stopwatch(void)2141 static double xtimer_stopwatch (void)
2142 {
2143 struct xtimer *xt = get_xtimer();
2144 double t1 = get_xtime();
2145 double ret = t1 - xt->t0;
2146
2147 xt->t0 = t1;
2148
2149 return ret;
2150 }
2151
2152 #endif /* WANT_XTIMER */
2153
2154 #ifdef NEED_ITIMER
2155
2156 /* Alternative timer, giving time in microseconds as a 64-bit int.
2157 If we don't have OpenMP we'll need this when not under MPI.
2158 */
2159
2160 struct itimer {
2161 int level;
2162 gint64 t0;
2163 };
2164
2165 static GPtrArray *itimers;
2166
new_itimer(int level)2167 static struct itimer *new_itimer (int level)
2168 {
2169 struct itimer *it = g_malloc(sizeof *it);
2170
2171 it->level = level;
2172 it->t0 = g_get_monotonic_time();
2173 g_ptr_array_insert(itimers, itimers->len, it);
2174 return it;
2175 }
2176
itimer_init(void)2177 static void itimer_init (void)
2178 {
2179 if (itimers == NULL) {
2180 itimers = g_ptr_array_sized_new(1);
2181 new_itimer(gretl_function_depth());
2182 }
2183 }
2184
get_itimer(void)2185 static struct itimer *get_itimer (void)
2186 {
2187 struct itimer *it;
2188 int i, lev = gretl_function_depth();
2189
2190 if (itimers == NULL) {
2191 itimers = g_ptr_array_sized_new(1);
2192 return new_itimer(lev);
2193 }
2194
2195 for (i=0; i<itimers->len; i++) {
2196 it = g_ptr_array_index(itimers, i);
2197 if (it->level == lev) {
2198 return it;
2199 }
2200 }
2201
2202 return new_itimer(lev);
2203 }
2204
itimer_stopwatch(void)2205 static double itimer_stopwatch (void)
2206 {
2207 struct itimer *it = get_itimer();
2208 gint64 t1 = g_get_monotonic_time();
2209 double ret = (t1 - it->t0) * 1.0e-6;
2210
2211 it->t0 = t1;
2212
2213 return ret;
2214 }
2215
2216 #endif /* NEED_ITIMER */
2217
gretl_stopwatch_cleanup(void)2218 static void gretl_stopwatch_cleanup (void)
2219 {
2220 int i;
2221
2222 #ifdef WANT_XTIMER
2223 if (xtimers != NULL) {
2224 for (i=0; i<xtimers->len; i++) {
2225 g_free(xtimers->pdata[i]);
2226 }
2227 g_ptr_array_free(xtimers, TRUE);
2228 xtimers = NULL;
2229 }
2230 #endif
2231
2232 #ifdef NEED_ITIMER
2233 if (itimers != NULL) {
2234 for (i=0; i<itimers->len; i++) {
2235 g_free(itimers->pdata[i]);
2236 }
2237 g_ptr_array_free(itimers, TRUE);
2238 itimers = NULL;
2239 }
2240 #endif
2241 }
2242
gretl_stopwatch(void)2243 double gretl_stopwatch (void)
2244 {
2245 #if defined(HAVE_MPI)
2246 if (gretl_mpi_initialized()) {
2247 return xtimer_stopwatch();
2248 }
2249 #endif
2250 #if defined(_OPENMP)
2251 return xtimer_stopwatch();
2252 #else
2253 return itimer_stopwatch();
2254 #endif
2255 }
2256
gretl_stopwatch_init(void)2257 static void gretl_stopwatch_init (void)
2258 {
2259 #if defined(HAVE_MPI)
2260 if (gretl_mpi_initialized()) {
2261 xtimer_init();
2262 return;
2263 }
2264 #endif
2265 #if defined(_OPENMP)
2266 xtimer_init();
2267 #else
2268 itimer_init();
2269 #endif
2270 }
2271
2272 /* BLAS detection and analysis */
2273
2274 static char OB_core[32];
2275 static char OB_parallel[12];
2276
2277 static int blas_variant;
2278
parse_ldd_output(const char * s)2279 static int parse_ldd_output (const char *s)
2280 {
2281 char found[4] = {0};
2282 char line[512];
2283 int i = 0;
2284 int ret = BLAS_UNKNOWN;
2285
2286 *line = '\0';
2287
2288 while (*s) {
2289 if (*s == '\n') {
2290 /* got to the end of a line */
2291 line[i] = '\0';
2292 if (strstr(line, "Accelerate.frame")) {
2293 found[3] = 1;
2294 } else if (strstr(line, "libmkl")) {
2295 found[2] = 1;
2296 } else if (strstr(line, "atlas")) {
2297 found[1] = 1;
2298 } else if (strstr(line, "libblas")) {
2299 found[0] = 1;
2300 }
2301 *line = '\0';
2302 i = 0;
2303 } else {
2304 line[i++] = *s;
2305 }
2306 s++;
2307 }
2308
2309 if (found[3]) {
2310 ret = BLAS_VECLIB;
2311 } else if (found[2]) {
2312 ret = BLAS_MKL;
2313 } else if (found[1]) {
2314 ret = BLAS_ATLAS;
2315 } else if (found[0]) {
2316 ret = BLAS_NETLIB;
2317 } else {
2318 fputs("detect blas: found no relevant libs!\n", stderr);
2319 }
2320
2321 return ret;
2322 }
2323
2324 #if !defined(WIN32)
2325
detect_blas_via_ldd(void)2326 static int detect_blas_via_ldd (void)
2327 {
2328 gchar *prog;
2329 gchar *sout = NULL;
2330 gchar *errout = NULL;
2331 gint status = 0;
2332 GError *gerr = NULL;
2333 int variant = 0;
2334
2335 #ifdef OS_OSX
2336 gchar *argv[4];
2337 prog = g_strdup("gretlcli"); /* ?? */
2338 argv[0] = "otool";
2339 argv[1] = "-L";
2340 argv[2] = prog;
2341 argv[3] = NULL;
2342 #else
2343 gchar *argv[3];
2344 prog = g_strdup(GRETL_PREFIX "/bin/gretlcli");
2345 argv[0] = "ldd";
2346 argv[1] = prog;
2347 argv[2] = NULL;
2348 #endif
2349
2350 g_spawn_sync(NULL, argv, NULL, G_SPAWN_SEARCH_PATH,
2351 NULL, NULL, &sout, &errout,
2352 &status, &gerr);
2353
2354 if (gerr != NULL) {
2355 fprintf(stderr, "%s\n", gerr->message);
2356 g_error_free(gerr);
2357 } else if (status != 0) {
2358 fprintf(stderr, "%s exited with status %d\n", argv[0], status);
2359 } else if (sout != NULL) {
2360 variant = parse_ldd_output(sout);
2361 } else {
2362 fprintf(stderr, "%s: %s\n", argv[0], "Got no output");
2363 }
2364
2365 g_free(sout);
2366 g_free(errout);
2367 g_free(prog);
2368
2369 return variant;
2370 }
2371
2372 #endif /* neither Windows nor Mac */
2373
2374 #include <dlfcn.h>
2375
register_openblas_details(void * handle)2376 static void register_openblas_details (void *handle)
2377 {
2378 char *(*OB_get_corename) (void);
2379 int (*OB_get_parallel) (void);
2380
2381 OB_get_corename = dlsym(handle, "openblas_get_corename");
2382 OB_get_parallel = dlsym(handle, "openblas_get_parallel");
2383
2384 if (OB_get_corename != NULL) {
2385 char *s = OB_get_corename();
2386
2387 if (s != NULL) {
2388 *OB_core = '\0';
2389 strncat(OB_core, s, 31);
2390 }
2391 } else {
2392 fprintf(stderr, "Couldn't find openblas_get_corename()\n");
2393 }
2394
2395 if (OB_get_parallel != NULL) {
2396 int p = OB_get_parallel();
2397
2398 if (p == 0) {
2399 strcpy(OB_parallel, "none");
2400 } else if (p == 1) {
2401 strcpy(OB_parallel, "pthreads");
2402 } else if (p == 2) {
2403 strcpy(OB_parallel, "OpenMP");
2404 }
2405 } else {
2406 fprintf(stderr, "Couldn't find openblas_get_parallel()\n");
2407 }
2408 }
2409
2410 /* below: called in creating $sysinfo bundle */
2411
get_openblas_details(char ** s1,char ** s2)2412 int get_openblas_details (char **s1, char **s2)
2413 {
2414 if (*OB_core == '\0' || *OB_parallel == '\0') {
2415 return 0;
2416 } else {
2417 *s1 = OB_core;
2418 *s2 = OB_parallel;
2419 return 1;
2420 }
2421 }
2422
2423 /* for $sysinfo bundle */
2424
blas_variant_string(void)2425 const char *blas_variant_string (void)
2426 {
2427 if (blas_variant == BLAS_NETLIB) {
2428 return "netlib";
2429 } else if (blas_variant == BLAS_ATLAS) {
2430 return "atlas";
2431 } else if (blas_variant == BLAS_OPENBLAS) {
2432 return "openblas";
2433 } else if (blas_variant == BLAS_MKL) {
2434 return "mkl";
2435 } else if (blas_variant == BLAS_VECLIB) {
2436 return "veclib";
2437 } else {
2438 return "unknown";
2439 }
2440 }
2441
2442 /* for gretl_matrix.c, libset.c */
2443
blas_is_openblas(void)2444 int blas_is_openblas (void)
2445 {
2446 return blas_variant == BLAS_OPENBLAS;
2447 }
2448
2449 static void (*OB_set_num_threads) (int);
2450 static int (*OB_get_num_threads) (void);
2451
blas_set_num_threads(int nt)2452 void blas_set_num_threads (int nt)
2453 {
2454 if (OB_set_num_threads != NULL) {
2455 OB_set_num_threads(nt);
2456 }
2457 }
2458
blas_get_num_threads(void)2459 int blas_get_num_threads (void)
2460 {
2461 if (OB_get_num_threads != NULL) {
2462 return OB_get_num_threads();
2463 } else {
2464 return 0;
2465 }
2466 }
2467
blas_init(void)2468 static void blas_init (void)
2469 {
2470 void *ptr = NULL;
2471
2472 #if defined(OS_OSX) && defined(PKGBUILD)
2473 blas_variant = BLAS_VECLIB;
2474 return;
2475 #endif
2476
2477 ptr = dlopen(NULL, RTLD_NOW);
2478
2479 if (ptr != NULL) {
2480 OB_set_num_threads = dlsym(ptr, "openblas_set_num_threads");
2481 OB_get_num_threads = dlsym(ptr, "openblas_get_num_threads");
2482 if (OB_set_num_threads != NULL) {
2483 blas_variant = BLAS_OPENBLAS;
2484 register_openblas_details(ptr);
2485 }
2486 }
2487
2488 if (blas_variant != BLAS_OPENBLAS) {
2489 #ifdef WIN32
2490 blas_variant = BLAS_NETLIB; /* ?? */
2491 #else
2492 blas_variant = detect_blas_via_ldd();
2493 #endif
2494 }
2495 }
2496
2497 /* library init and cleanup functions */
2498
2499 /**
2500 * libgretl_init:
2501 *
2502 * In a program that uses libgretl, this function should be
2503 * called once, before any other libgretl functions are
2504 * used. See also libgretl_cleanup(), and libgretl_mpi_init().
2505 **/
2506
libgretl_init(void)2507 void libgretl_init (void)
2508 {
2509 libset_init();
2510 gretl_rand_init();
2511 gretl_xml_init();
2512 gretl_stopwatch_init();
2513 if (!gretl_in_tool_mode()) {
2514 blas_init();
2515 num_threads_init(blas_variant);
2516 }
2517 #if HAVE_GMP
2518 mpf_set_default_prec(get_mp_bits());
2519 #endif
2520 }
2521
2522 #ifdef HAVE_MPI
2523
2524 /**
2525 * libgretl_mpi_init:
2526 * @self: the MPI rank of the calling process.
2527 * @np: the number of MPI processes.
2528 * @dcmt: if non-zero, set up per-process RNG using DCMT.
2529 *
2530 * This function provides an alternative to libgretl_init()
2531 * which should be used when a libgretl program is to be run in
2532 * MPI mode.
2533 **/
2534
libgretl_mpi_init(int self,int np,int dcmt)2535 int libgretl_mpi_init (int self, int np, int dcmt)
2536 {
2537 int err;
2538
2539 libset_init();
2540
2541 /* let geneval know */
2542 set_mpi_rank_and_size(self, np);
2543
2544 /* load MPI library functions */
2545 err = gretl_MPI_init();
2546 if (err) {
2547 return err;
2548 }
2549
2550 if (dcmt && np > 1) {
2551 /* use DCMT for multiple RNGs */
2552 gretl_dcmt_init(np, self, 0);
2553 } else {
2554 /* just use one RNG */
2555 gretl_rand_init();
2556 }
2557
2558 gretl_xml_init();
2559 gretl_stopwatch_init();
2560 #if HAVE_GMP
2561 mpf_set_default_prec(get_mp_bits());
2562 #endif
2563
2564 /* be relatively quiet by default */
2565 set_gretl_echo(0);
2566 set_gretl_messages(0);
2567
2568 return 0;
2569 }
2570
2571 #else
2572
gretl_mpi_initialized(void)2573 int gretl_mpi_initialized (void)
2574 {
2575 /* stub for non-MPI build */
2576 return 0;
2577 }
2578
2579 #endif
2580
other_gretl_running(const char * fname)2581 static int other_gretl_running (const char *fname)
2582 {
2583 FILE *fp = gretl_fopen(fname, "rb");
2584 int ret = 0;
2585
2586 if (fp != NULL) {
2587 /* Each entry in gretl.pid will contain two
2588 integers: PID and sequence number. So if
2589 we can read more than two values that must
2590 mean there's more than one gretl process
2591 running.
2592 */
2593 long lv[4];
2594 int np;
2595
2596 np = fscanf(fp, "%ld %ld %ld %ld", &lv[0],
2597 &lv[1], &lv[2], &lv[3]);
2598 if (np > 2) {
2599 ret = 1;
2600 }
2601 fclose(fp);
2602 }
2603
2604 return ret;
2605 }
2606
dotdir_cleanup(void)2607 static void dotdir_cleanup (void)
2608 {
2609 const char *dotdir = gretl_dotdir();
2610 const char *workdir = gretl_workdir();
2611 int err;
2612
2613 if (!strcmp(dotdir, workdir)) {
2614 return;
2615 }
2616
2617 err = gretl_chdir(dotdir);
2618
2619 if (!err) {
2620 GDir *dir = gretl_opendir(".");
2621 const gchar *fname;
2622 int skipit = 0;
2623
2624 if (dir != NULL) {
2625 while ((fname = g_dir_read_name(dir)) != NULL) {
2626 if (!strcmp(fname, "gretl.pid")) {
2627 skipit = other_gretl_running(fname);
2628 break;
2629 }
2630 }
2631 if (!skipit) {
2632 g_dir_rewind(dir);
2633 while ((fname = g_dir_read_name(dir)) != NULL) {
2634 if (gretl_isdir(fname)) {
2635 if (*fname == '.' && strlen(fname) > 3) {
2636 /* failed auto-dot directory? */
2637 gretl_deltree(fname);
2638 }
2639 } else if (strcmp(fname, "..") &&
2640 strcmp(fname, ".") &&
2641 strcmp(fname, ".gretl2rc") &&
2642 strcmp(fname, "gretl.pid") &&
2643 strcmp(fname, "addons.idx") &&
2644 strcmp(fname, "mail.dat")) {
2645 gretl_remove(fname);
2646 }
2647 }
2648 }
2649 g_dir_close(dir);
2650 }
2651 }
2652 }
2653
libgretl_session_cleanup(int mode)2654 void libgretl_session_cleanup (int mode)
2655 {
2656 /* "last model" for accessors */
2657 gretl_saved_objects_cleanup();
2658
2659 /* trash dataset-related items */
2660 gretl_transforms_cleanup();
2661 gretl_lists_cleanup();
2662 gretl_tests_cleanup();
2663 gretl_plotx(NULL, OPT_NONE);
2664
2665 /* scrub options set via "setopt" */
2666 setopt_cleanup();
2667
2668 if (mode != SESSION_CLEAR_DATASET) {
2669 destroy_user_vars();
2670 }
2671 }
2672
2673 /**
2674 * libgretl_cleanup:
2675 *
2676 * In a program that uses libgretl, this function may be
2677 * called to free various chunks of memory after the program
2678 * is finished with libgretl. Do not attempt to call any
2679 * other libgretl functions after invoking this cleanup.
2680 *
2681 * See also libgretl_init().
2682 **/
2683
libgretl_cleanup(void)2684 void libgretl_cleanup (void)
2685 {
2686 libgretl_session_cleanup(SESSION_CLEAR_ALL);
2687
2688 gretl_rand_free();
2689 gretl_functions_cleanup();
2690 libset_cleanup();
2691 gretl_command_hash_cleanup();
2692 gretl_function_hash_cleanup();
2693 lapack_mem_free();
2694 forecast_matrix_cleanup();
2695 stored_options_cleanup();
2696 option_printing_cleanup();
2697 gnuplot_cleanup();
2698 bufgets_cleanup();
2699 plugins_cleanup();
2700 gretl_bundle_cleanup();
2701 gretl_typemap_cleanup();
2702 gretl_stopwatch_cleanup();
2703 #ifdef USE_CURL
2704 gretl_www_cleanup();
2705 #endif
2706 builtin_strings_cleanup();
2707 last_result_cleanup();
2708
2709 #ifdef HAVE_MPI
2710 if (!gretl_mpi_initialized()) {
2711 dotdir_cleanup();
2712 }
2713 #else
2714 dotdir_cleanup();
2715 #endif
2716
2717 #ifdef USE_RLIB
2718 gretl_R_cleanup();
2719 #endif
2720
2721 gretl_script_dirs_cleanup();
2722 gretl_xml_cleanup();
2723 }
2724
2725 /* record and retrieve hypothesis test results */
2726
2727 enum {
2728 SET_TEST_STAT,
2729 GET_TEST_STAT,
2730 GET_TEST_PVAL,
2731 GET_TEST_LNL,
2732 GET_TEST_BRK,
2733 TESTS_CLEANUP
2734 };
2735
2736 #define getcode(c) (c == GET_TEST_STAT || c == GET_TEST_PVAL || \
2737 c == GET_TEST_LNL || c == GET_TEST_BRK)
2738
2739 static GretlType last_test_type;
2740
record_or_get_test_result(double teststat,double pval,double lnl,double inbrk,int code)2741 static double record_or_get_test_result (double teststat,
2742 double pval,
2743 double lnl,
2744 double inbrk,
2745 int code)
2746 {
2747 static double val = NADBL;
2748 static double pv = NADBL;
2749 static double ll = NADBL;
2750 static double brk = NADBL;
2751 double ret = NADBL;
2752
2753 if (code == SET_TEST_STAT) {
2754 last_test_type = GRETL_TYPE_DOUBLE;
2755 val = teststat;
2756 pv = pval;
2757 ll = lnl;
2758 brk = inbrk;
2759 } else if (getcode(code) && last_test_type == GRETL_TYPE_DOUBLE) {
2760 if (code == GET_TEST_STAT) {
2761 ret = val;
2762 } else if (code == GET_TEST_PVAL) {
2763 ret = pv;
2764 } else if (code == GET_TEST_LNL) {
2765 ret = ll;
2766 } else if (code == GET_TEST_BRK) {
2767 ret = brk;
2768 }
2769 }
2770
2771 return ret;
2772 }
2773
2774 static gretl_matrix *
record_or_get_test_matrix(gretl_matrix * tests,gretl_matrix * pvals,int code,int * err)2775 record_or_get_test_matrix (gretl_matrix *tests,
2776 gretl_matrix *pvals,
2777 int code, int *err)
2778 {
2779 static gretl_matrix *vals = NULL;
2780 static gretl_matrix *pvs = NULL;
2781 gretl_matrix *ret = NULL;
2782
2783 if (code == TESTS_CLEANUP) {
2784 gretl_matrix_free(vals);
2785 gretl_matrix_free(pvs);
2786 vals = pvs = NULL;
2787 last_test_type = GRETL_TYPE_NONE;
2788 return NULL;
2789 }
2790
2791 if (code == SET_TEST_STAT) {
2792 last_test_type = GRETL_TYPE_MATRIX;
2793 gretl_matrix_free(vals);
2794 vals = tests;
2795 gretl_matrix_free(pvs);
2796 pvs = pvals;
2797 } else if (getcode(code)) {
2798 gretl_matrix *src = (code == GET_TEST_STAT)? vals : pvs;
2799
2800 if (src != NULL) {
2801 ret = gretl_matrix_copy(src);
2802 if (ret == NULL) {
2803 *err = E_ALLOC;
2804 }
2805 } else {
2806 *err = E_BADSTAT;
2807 }
2808 }
2809
2810 return ret;
2811 }
2812
gretl_tests_cleanup(void)2813 static void gretl_tests_cleanup (void)
2814 {
2815 record_or_get_test_matrix(NULL, NULL, TESTS_CLEANUP, NULL);
2816 }
2817
2818 /* Returns the type (scalar, matrix or none) of the last recorded
2819 test statistic/p-value pair.
2820 */
2821
get_last_test_type(void)2822 int get_last_test_type (void)
2823 {
2824 return last_test_type;
2825 }
2826
record_test_result(double teststat,double pval)2827 void record_test_result (double teststat, double pval)
2828 {
2829 record_or_get_test_result(teststat, pval, NADBL, NADBL, SET_TEST_STAT);
2830 }
2831
record_LR_test_result(double teststat,double pval,double lnl)2832 void record_LR_test_result (double teststat, double pval, double lnl)
2833 {
2834 record_or_get_test_result(teststat, pval, lnl, NADBL, SET_TEST_STAT);
2835 }
2836
record_QLR_test_result(double teststat,double pval,double brk)2837 void record_QLR_test_result (double teststat, double pval, double brk)
2838 {
2839 record_or_get_test_result(teststat, pval, NADBL, brk, SET_TEST_STAT);
2840 }
2841
2842 /* Note: the hypothesis-test recorder "takes ownership" of the
2843 two input matrices, @tests and @pvals, which should therefore
2844 NOT be freed by the caller.
2845 */
2846
record_matrix_test_result(gretl_matrix * tests,gretl_matrix * pvals)2847 void record_matrix_test_result (gretl_matrix *tests,
2848 gretl_matrix *pvals)
2849 {
2850 record_or_get_test_matrix(tests, pvals, SET_TEST_STAT, NULL);
2851 }
2852
get_last_test_statistic(void)2853 double get_last_test_statistic (void)
2854 {
2855 return record_or_get_test_result(0, 0, 0, 0, GET_TEST_STAT);
2856 }
2857
get_last_pvalue(void)2858 double get_last_pvalue (void)
2859 {
2860 return record_or_get_test_result(0, 0, 0, 0, GET_TEST_PVAL);
2861 }
2862
get_last_lnl(void)2863 double get_last_lnl (void)
2864 {
2865 return record_or_get_test_result(0, 0, 0, 0, GET_TEST_LNL);
2866 }
2867
get_last_break(void)2868 double get_last_break (void)
2869 {
2870 return record_or_get_test_result(0, 0, 0, 0, GET_TEST_BRK);
2871 }
2872
get_last_test_matrix(int * err)2873 gretl_matrix *get_last_test_matrix (int *err)
2874 {
2875 return record_or_get_test_matrix(NULL, NULL, GET_TEST_STAT, err);
2876 }
2877
get_last_pvals_matrix(int * err)2878 gretl_matrix *get_last_pvals_matrix (int *err)
2879 {
2880 return record_or_get_test_matrix(NULL, NULL, GET_TEST_PVAL, err);
2881 }
2882
2883 #define NEED_ALIGNED_MALLOC 0
2884
2885 #if NEED_ALIGNED_MALLOC
2886
2887 /* We needed this -- in the absence of posix_memalign -- when
2888 experimenting with the array variant of SFMT (random.c). Since that
2889 did not seem to be faster than the "one at a time" variant when we
2890 last experimented, we junked the array-using code, which makes the
2891 following functionality redundant. However, we may want to try
2892 using SFMT arrays again. In that case, see git commit
2893 65f3395bbdb01450496fb20fcf0f8f15092f9eeb for the last revision of
2894 random.c that included the array code.
2895 */
2896
2897 /*
2898 malloc and free for alignments greater than that guaranteed by the C
2899 library, based on Steven G. Johnson's public domand code at
2900 http://ab-initio.mit.edu/~stevenj/align.c
2901 */
2902
2903 #ifdef HAVE_STDINT_H
2904 # include <stdint.h> /* for uintptr_t */
2905 #else
2906 # define uintptr_t size_t
2907 #endif
2908
2909 #define NOT_POWER_OF_TWO(n) (((n) & ((n) - 1)))
2910 #define UI(p) ((uintptr_t) (p))
2911
2912 #define PTR_ALIGN(p0, align) \
2913 ((void *) (((UI(p0) + (align + sizeof(void*))) \
2914 & (~UI(align - 1)))))
2915
2916 /* pointer must sometimes be aligned; assume sizeof(void*) is a power of two */
2917 #define ORIG_PTR(p) (*(((void **) (UI(p) & (~UI(sizeof(void*) - 1)))) - 1))
2918
real_aligned_malloc(size_t size,size_t align)2919 static void *real_aligned_malloc (size_t size, size_t align)
2920 {
2921 void *p0, *p;
2922
2923 if (NOT_POWER_OF_TWO(align)) {
2924 errno = EINVAL;
2925 return NULL;
2926 }
2927
2928 if (align < sizeof(void *)) {
2929 align = sizeof(void *);
2930 }
2931
2932 /* including the extra sizeof(void*) is overkill on a 32-bit
2933 machine, since malloc is already 8-byte aligned, as long
2934 as we enforce alignment >= 8 ...but oh well */
2935
2936 p0 = malloc(size + align + sizeof(void *));
2937 if (!p0) {
2938 return NULL;
2939 }
2940
2941 p = PTR_ALIGN(p0, align);
2942 ORIG_PTR(p) = p0;
2943
2944 return p;
2945 }
2946
gretl_aligned_malloc(size_t size,size_t alignment)2947 void *gretl_aligned_malloc (size_t size, size_t alignment)
2948 {
2949 if (size < 1) {
2950 return NULL;
2951 } else {
2952 return real_aligned_malloc(size, alignment);
2953 }
2954 }
2955
gretl_aligned_free(void * mem)2956 void gretl_aligned_free (void *mem)
2957 {
2958 if (mem != NULL) {
2959 free(ORIG_PTR(mem));
2960 }
2961 }
2962
2963 #endif /* NEED_ALIGNED_MALLOC */
2964
2965 #ifdef WIN32
2966
check_for_program(const char * prog)2967 int check_for_program (const char *prog)
2968 {
2969 int ret = 0;
2970
2971 if (prog == NULL || *prog == '\0') {
2972 return 0;
2973 }
2974
2975 if (has_suffix(prog, ".exe")) {
2976 ret = win32_check_for_program(prog);
2977 } else {
2978 gchar *test = g_strdup_printf("%s.exe", prog);
2979
2980 ret = win32_check_for_program(test);
2981 g_free(test);
2982 }
2983
2984 return ret;
2985 }
2986
2987 #else /* !WIN32 */
2988
2989 #include <sys/stat.h>
2990 #include <unistd.h>
2991
is_executable(const char * s,uid_t myid,gid_t mygrp)2992 static int is_executable (const char *s, uid_t myid, gid_t mygrp)
2993 {
2994 struct stat buf;
2995 int ok = 0;
2996
2997 if (gretl_stat(s, &buf) != 0) {
2998 return 0;
2999 }
3000
3001 if (buf.st_mode & (S_IFREG|S_IFLNK)) {
3002 if (buf.st_uid == myid && (buf.st_mode & S_IXUSR)) {
3003 ok = 1;
3004 } else if (buf.st_gid == mygrp && (buf.st_mode & S_IXGRP)) {
3005 ok = 1;
3006 } else if (buf.st_uid != myid && buf.st_gid != mygrp &&
3007 (buf.st_mode & S_IXOTH)) {
3008 ok = 1;
3009 }
3010 }
3011
3012 return ok;
3013 }
3014
3015 /**
3016 * check_for_program:
3017 * @prog: name of program.
3018 *
3019 * Returns: 1 if @prog is found in the PATH, otherwise 0.
3020 */
3021
check_for_program(const char * prog)3022 int check_for_program (const char *prog)
3023 {
3024 uid_t myid = getuid();
3025 gid_t mygrp = getgid();
3026 char *path;
3027 char *pathcpy;
3028 char **dirs;
3029 char *fullpath;
3030 char *p;
3031 int max_dlen = 0;
3032 int found = 0;
3033 int i, ndirs;
3034
3035 if (prog == NULL || *prog == '\0') {
3036 return 0;
3037 }
3038
3039 if (*prog == '/') {
3040 return is_executable(prog, myid, mygrp);
3041 }
3042
3043 path = getenv("PATH");
3044 if (path == NULL || *path == '\0') {
3045 return 0;
3046 }
3047
3048 pathcpy = gretl_strdup(path);
3049 if (pathcpy == NULL) {
3050 return 0;
3051 }
3052
3053 ndirs = 1;
3054 p = pathcpy;
3055 while (*p) {
3056 if (*p == ':') ndirs++;
3057 p++;
3058 }
3059
3060 dirs = malloc(ndirs * sizeof *dirs);
3061 if (dirs == NULL) {
3062 free(pathcpy);
3063 return 0;
3064 }
3065
3066 if (ndirs == 1) {
3067 dirs[0] = pathcpy;
3068 max_dlen = strlen(pathcpy);
3069 } else {
3070 for (i=0; i<ndirs; i++) {
3071 int dlen;
3072
3073 dirs[i] = strtok((i == 0)? pathcpy : NULL, ":");
3074 if (dirs[i] == NULL) {
3075 ndirs = i;
3076 break;
3077 }
3078 dlen = strlen(dirs[i]);
3079 if (dlen > max_dlen) {
3080 max_dlen = dlen;
3081 }
3082 }
3083 }
3084
3085 if (ndirs == 0 ||
3086 (fullpath = malloc(max_dlen + strlen(prog) + 2)) == NULL) {
3087 free(dirs);
3088 free(pathcpy);
3089 return 0;
3090 }
3091
3092 for (i=0; i<ndirs && !found; i++) {
3093 sprintf(fullpath, "%s/%s", dirs[i], prog);
3094 found = is_executable(fullpath, myid, mygrp);
3095 }
3096
3097 free(dirs);
3098 free(pathcpy);
3099 free(fullpath);
3100
3101 return found;
3102 }
3103
3104 #endif /* WIN32 or not */
3105
gretl_monotonic_time(void)3106 gint64 gretl_monotonic_time (void)
3107 {
3108 return g_get_monotonic_time();
3109 }
3110
3111 /* implementations of gzip and gunzip for unitary files */
3112
3113 #define GRETL_ZBUFSIZ 131072 /* 128k bytes */
3114
gretl_gzip(char * fname,char * zname)3115 int gretl_gzip (char *fname, char *zname)
3116 {
3117 char buf[GRETL_ZBUFSIZ];
3118 size_t len;
3119 FILE *fp;
3120 gzFile fz;
3121
3122 fp = gretl_fopen(fname, "rb");
3123 if (fp == NULL) {
3124 return E_FOPEN;
3125 }
3126
3127 fz = gretl_gzopen(zname, "wb");
3128 if (fz == Z_NULL) {
3129 fclose(fp);
3130 return E_FOPEN;
3131 }
3132
3133 while ((len = fread(buf, 1, GRETL_ZBUFSIZ, fp)) > 0) {
3134 gzwrite(fz, buf, len);
3135 }
3136
3137 fclose(fp);
3138 gzclose(fz);
3139
3140 return 0;
3141 }
3142
gretl_gunzip(char * zname,char * fname)3143 int gretl_gunzip (char *zname, char *fname)
3144 {
3145 char buf[GRETL_ZBUFSIZ];
3146 size_t len;
3147 FILE *fp;
3148 gzFile fz;
3149
3150 fz = gretl_gzopen(zname, "rb");
3151 if (fz == Z_NULL) {
3152 return E_FOPEN;
3153 }
3154
3155 fp = gretl_fopen(fname, "wb");
3156 if (fp == NULL) {
3157 gzclose(fz);
3158 return E_FOPEN;
3159 }
3160
3161 while ((len = gzread(fz, buf, GRETL_ZBUFSIZ)) > 0) {
3162 fwrite(buf, 1, len, fp);
3163 }
3164
3165 gzclose(fz);
3166 fclose(fp);
3167
3168 return 0;
3169 }
3170
3171 /* Check if we're able to launch an automatic MPI routine, on the
3172 local machine, from within some time-consuming function. This
3173 requires that MPI is available but not already running, and
3174 the local machine has at least 2 processors.
3175 */
3176
auto_mpi_ok(void)3177 int auto_mpi_ok (void)
3178 {
3179 int ret = 0;
3180
3181 #ifdef HAVE_MPI
3182 if (gretl_mpi_initialized()) {
3183 ; /* No: can't run MPI under MPI */
3184 } else if (gretl_n_processors() < 2) {
3185 ; /* No: can't do local MPI */
3186 } else {
3187 /* Yes, if mpiexec is installed */
3188 ret = check_for_mpiexec();
3189 }
3190 #endif
3191
3192 return ret;
3193 }
3194