?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
GnmFunc: make this a GObject.
[gnumeric.git] / plugins / fn-tsa / functions.c
blobdf1c74f0c8772f70b8a9dd5d660f15b7e1831779
1 /*
2 * fn-tsa plugin
3 * functions.c
4 *
5 * Copyright (C) 2006 Laurency Franck
6 * Copyright (C) 2007 Jean Bréfort <jean.brefort@normalesup.org>
7 * Copyright (C) 2009 Andreas Guelzow <aguelzow@pyrshep.ca>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License as
11 * published by the Free Software Foundation; either version 2 of the
12 * License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, see <https://www.gnu.org/licenses/>
21 */
22
23
24 /***************************************************************
25 * This file contains all declarations of time series analysis
26 * functions and their help functions
27 ****************************************************************/
28
29 #include <gnumeric-config.h>
30 #include <gnumeric.h>
31 #include <func.h>
32 #include <parse-util.h>
33 #include <mathfunc.h>
34 #include <rangefunc.h>
35 #include <regression.h>
36 #include <workbook.h>
37 #include <sheet.h>
38 #include <cell.h>
39 #include <collect.h>
40 #include <number-match.h>
41 #include <value.h>
42 #include <expr.h>
43 #include <func-builtin.h>
44 #include <gnm-i18n.h>
45 #include <goffice/goffice.h>
46 #include <gnm-plugin.h>
47 #include <tools/analysis-tools.h>
48 #include <complex.h>
49
50 #include <math.h>
51 #include <stdlib.h>
52 #include <string.h>
53
54 GNM_PLUGIN_MODULE_HEADER;
55
56 enum {
57 FILTER_NONE,
58 FILTER_BARTLETT,
59 FILTER_HANN,
60 FILTER_WELCH,
61 FILTER_GAUSSIAN,
62 };
63
64 enum {
65 INTERPOLATION_LINEAR = 0,
66 INTERPOLATION_LINEAR_AVG,
67 INTERPOLATION_STAIRCASE,
68 INTERPOLATION_STAIRCASE_AVG,
69 INTERPOLATION_SPLINE,
70 INTERPOLATION_SPLINE_AVG,
71 };
72
73 #ifdef GNM_WITH_LONG_DOUBLE
74 # define GnmCSpline GOCSplinel
75 # define gnm_cspline_init go_cspline_initl
76 # define gnm_cspline_destroy go_cspline_destroyl
77 # define gnm_cspline_get_value go_cspline_get_valuel
78 # define gnm_cspline_get_values go_cspline_get_valuesl
79 # define gnm_cspline_get_integrals go_cspline_get_integralsl
80 #else
81 # define GnmCSpline GOCSpline
82 # define gnm_cspline_init go_cspline_init
83 # define gnm_cspline_destroy go_cspline_destroy
84 # define gnm_cspline_get_value go_cspline_get_value
85 # define gnm_cspline_get_values go_cspline_get_values
86 # define gnm_cspline_get_integrals go_cspline_get_integrals
87 #endif
88
89 #define INTERPOLATIONMETHODS { GNM_FUNC_HELP_DESCRIPTION, F_("Possible interpolation methods are:\n" \
90 "0: linear;\n" \
91 "1: linear with averaging;\n" \
92 "2: staircase;\n" \
93 "3: staircase with averaging;\n" \
94 "4: natural cubic spline;\n" \
95 "5: natural cubic spline with averaging.") }
96
97 /**************************************************************************/
98 /**************************************************************************/
99 /* CALCULATION FUNCTIONS FOR TOOLS */
100 /**************************************************************************/
101 /**************************************************************************/
102
103 static void
104 gnm_fourier_fft (gnm_complex const *in, int n, int skip, gnm_complex **fourier, gboolean inverse)
105 {
106 gnm_complex *fourier_1, *fourier_2;
107 int i;
108 int nhalf = n / 2;
109 gnm_float argstep;
110
111 *fourier = g_new (gnm_complex, n);
112
113 if (n == 1) {
114 (*fourier)[0] = in[0];
115 return;
116 }
117
118 gnm_fourier_fft (in, nhalf, skip * 2, &fourier_1, inverse);
119 gnm_fourier_fft (in + skip, nhalf, skip * 2, &fourier_2, inverse);
120
121 argstep = (inverse ? M_PIgnum : -M_PIgnum) / nhalf;
122 for (i = 0; i < nhalf; i++) {
123 gnm_complex dir, tmp;
124
125 dir = GNM_CPOLAR (1, argstep * i);
126 tmp = GNM_CMUL (fourier_2[i], dir);
127
128 (*fourier)[i] = GNM_CSCALE (GNM_CADD (fourier_1[i], tmp), 0.5);
129
130 (*fourier)[i + nhalf] = GNM_CSCALE (GNM_CSUB (fourier_1[i], tmp), 0.5);
131 }
132
133 g_free (fourier_1);
134 g_free (fourier_2);
135 }
136
137
138 /*******LINEAR INTERPOLATION*******/
139
140 static gnm_float*
141 linear_interpolation (const gnm_float *absc, const gnm_float *ord, int nb_knots,
142 const gnm_float *targets, int nb_targets)
143 {
144 int i, j, k, jmax = nb_knots - 1;
145 gnm_float slope, *res;
146 if (nb_knots < 2)
147 return NULL;
148 res = g_new (gnm_float, nb_targets);
149 if (gnm_range_increasing (targets, nb_targets)) {
150 j = 1;
151 k = 0;
152 slope = (ord[1] - ord[0]) / (absc[1] - absc[0]);
153 for (i = 0; i < nb_targets; i++) {
154 while (j < jmax && targets[i] > absc[j])
155 j++;
156 if (k < j - 1) {
157 k = j - 1;
158 slope = (ord[j] - ord[k]) / (absc[j] - absc[k]);
159 }
160 res[i] = (targets[i] - absc[k]) * slope + ord[k];
161 }
162 } else {
163 int l;
164 for (i = 0; i < nb_targets; i++) {
165 k = jmax - 1;
166 if (targets[i] >= absc[k])
167 res[i] = (targets[i] - absc[k]) *
168 (ord[jmax] - ord[k]) / (absc[jmax] - absc[k])
169 + ord[k];
170 else if (targets[i] <= absc[1])
171 res[i] = (targets[i] - absc[0]) *
172 (ord[1] - ord[0]) / (absc[1] - absc[0])
173 + ord[0];
174 else {
175 j = 1;
176 while (k > j + 1) {
177 l = (k + j) / 2;
178 if (targets[i] > absc[l])
179 j = l;
180 else
181 k = l;
182 }
183 res[i] = (targets[i] - absc[j]) *
184 (ord[k] - ord[j]) / (absc[k] - absc[j])
185 + ord[j];
186 }
187 }
188 }
189 return res;
190 }
191
192 /*******LINEAR AVERAGING*******/
193
194 static gnm_float*
195 linear_averaging (const gnm_float *absc, const gnm_float *ord, int nb_knots,
196 const gnm_float *targets, int nb_targets)
197 {
198 int i, j, k, jmax = nb_knots - 1;
199 gnm_float slope, *res, x0, x1;
200 if (nb_knots < 2 || !gnm_range_increasing (targets, nb_targets + 1))
201 return NULL;
202 res = g_new (gnm_float, nb_targets);
203 j = 1;
204 while (j < jmax && targets[0] > absc[j])
205 j++;
206 k = j - 1;
207 slope = (ord[j] - ord[k]) / (absc[j] - absc[k]) / 2.;
208 for (i = 1; i <= nb_targets; i++) {
209 if (targets[i] < absc[j] || j == jmax) {
210 x0 = targets[i - 1] - absc[k];
211 x1 = targets[i] - absc[k];
212 res[i - 1] = (x1 * (slope * x1 + ord[k])
213 - x0 * (slope * x0 + ord[k]))
214 / (x1 - x0);
215 continue;
216 }
217 x0 = targets[i - 1] - absc[k];
218 x1 = absc[j] - absc[k];
219 res[i - 1] = (x1 * (slope * x1 + ord[k])
220 - x0 * (slope * x0 + ord[k]));
221 while (j < jmax && targets[i] > absc[++j]) {
222 k++;
223 x0 = absc[j] - absc[k];
224 slope = (ord[j] - ord[k]) / x0 / 2.;
225 res[i - 1] += x0 * (slope * x0 + ord[k]);
226 }
227 if (j > k + 1) {
228 k = j - 1;
229 slope = (ord[j] - ord[k]) / (absc[j] - absc[k]) / 2.;
230 } else
231 k = j;
232 x0 = targets[i] - absc[k];
233 res[i - 1] += x0 * (slope * x0 + ord[k]);
234 res[i - 1] /= (targets[i] - targets[i - 1]);
235 }
236 return res;
237 }
238
239 /*******STAIRCASE INTERPOLATION*******/
240
241 static gnm_float*
242 staircase_interpolation (const gnm_float *absc, const gnm_float *ord, int nb_knots,
243 const gnm_float *targets, int nb_targets)
244 {
245 int i, j, jmax = nb_knots - 1;
246 gnm_float *res;
247
248 if (nb_knots <= 0)
249 return NULL;
250
251 res = g_new (gnm_float, nb_targets);
252 if (gnm_range_increasing (targets, nb_targets)) {
253 j = 1;
254 for (i = 0; i < nb_targets; i++) {
255 while (j <= jmax && targets[i] >= absc[j])
256 j++;
257 res[i] = ord[j - 1];
258 }
259 } else {
260 int k, l;
261 for (i = 0; i < nb_targets; i++) {
262 k = jmax;
263 if (targets[i] >= absc[jmax])
264 res[i] = ord[jmax];
265 else {
266 j = 0;
267 while (k > j + 1) {
268 l = (k + j) / 2;
269 if (targets[i] >= absc[l])
270 j = l;
271 else
272 k = l;
273 }
274 if (k != j && targets[i] >= absc[k])
275 j = k;
276 res[i] = ord[j];
277 }
278 }
279 }
280 return res;
281 }
282
283 /*******STAIRCASE AVERAGING*******/
284
285 static gnm_float*
286 staircase_averaging (const gnm_float *absc, const gnm_float *ord, int nb_knots,
287 const gnm_float *targets, int nb_targets)
288 {
289 int i, j, jmax = nb_knots - 1;
290 gnm_float *res;
291
292
293 if (nb_knots <= 0)
294 return NULL;
295
296 if (!gnm_range_increasing (targets, nb_targets + 1))
297 return NULL;
298 res = g_new (gnm_float, nb_targets);
299 j = 1;
300 while (j <= jmax && targets[0] >= absc[j])
301 j++;
302 for (i = 1; i <= nb_targets; i++) {
303 if (j > jmax || targets[i] < absc[j]) {
304 res[i - 1] = ord[j - 1];
305 continue;
306 }
307 res[i - 1] = (absc[j] - targets[i - 1]) * ord[j - 1];
308 while (j < jmax && targets[i] >= absc[++j])
309 res[i - 1] += (absc[j] - absc[j - 1]) * ord[j - 1];
310 if (targets[i] >= absc[j])
311 j++;
312 res[i - 1] += (targets[i] - absc[j - 1]) * ord[j - 1];
313 res[i - 1] /= (targets[i] - targets[i - 1]);
314 }
315 return res;
316 }
317
318 /*******SPLINE INTERPOLATION*******/
319
320 static gnm_float*
321 spline_interpolation (const gnm_float *absc, const gnm_float *ord, int nb_knots,
322 const gnm_float *targets, int nb_targets)
323 {
324 gnm_float *res;
325 int i;
326 GnmCSpline *sp = gnm_cspline_init (absc, ord, nb_knots,
327 GO_CSPLINE_NATURAL, 0., 0.);
328 if (!sp)
329 return NULL;
330 if (gnm_range_increasing (targets, nb_targets))
331 res = gnm_cspline_get_values (sp, targets, nb_targets);
332 else {
333 res = g_new (gnm_float, nb_targets);
334 for (i = 0; i < nb_targets; i++)
335 res[i] = gnm_cspline_get_value (sp, targets[i]);
336 }
337 gnm_cspline_destroy (sp);
338 return res;
339 }
340
341 /*******SPLINE AVERAGING*******/
342
343 static gnm_float *
344 spline_averaging (const gnm_float *absc, const gnm_float *ord, int nb_knots,
345 const gnm_float *targets, int nb_targets)
346 {
347 gnm_float *res;
348 int i, imax;
349 GnmCSpline *sp;
350 if (!gnm_range_increasing (targets, nb_targets + 1))
351 return NULL;
352 sp = gnm_cspline_init (absc, ord, nb_knots,
353 GO_CSPLINE_NATURAL, 0., 0.);
354 if (!sp)
355 return NULL;
356 res = gnm_cspline_get_integrals (sp, targets, nb_targets + 1);
357 imax = nb_targets;
358 for (i = 0; i < imax; i++)
359 res[i] /= targets[i + 1] - targets[i];
360 gnm_cspline_destroy (sp);
361 return res;
362 }
363
364 /*******Interpolation procedure********/
365
366 typedef gnm_float* (*INTERPPROC) (const gnm_float*, const gnm_float*,
367 int, const gnm_float*, int);
368
369 /******************************************************************************/
370 /* INTERPOLATION FUNCTION */
371 /******************************************************************************/
372
373 static GnmFuncHelp const help_interpolation[] = {
374 { GNM_FUNC_HELP_NAME, F_("INTERPOLATION:interpolated values corresponding to the given abscissa targets") },
375 { GNM_FUNC_HELP_ARG, F_("abscissae:abscissae of the given data points") },
376 { GNM_FUNC_HELP_ARG, F_("ordinates:ordinates of the given data points") },
377 { GNM_FUNC_HELP_ARG, F_("targets:abscissae of the interpolated data") },
378 { GNM_FUNC_HELP_ARG, F_("interpolation:method of interpolation, defaults to 0 (\'linear\')") },
379 { GNM_FUNC_HELP_DESCRIPTION, F_("The output consists always of one column of numbers.") },
380 INTERPOLATIONMETHODS,
381 { GNM_FUNC_HELP_NOTE, F_("The @{abscissae} should be given in increasing order. If the @{abscissae} is not in "
382 "increasing order the INTERPOLATION function is significantly slower.") },
383 { GNM_FUNC_HELP_NOTE, F_("If any two @{abscissae} values are equal an error is returned.") },
384 { GNM_FUNC_HELP_NOTE, F_("If any of interpolation methods 1 ('linear with averaging'), 3 "
385 "('staircase with averaging'), and 5 ('natural cubic spline with "
386 "averaging') is used, the number "
387 "of returned values is one less than the number of targets and the target "
388 "values must be given in increasing order. The values returned "
389 "are the average heights of the interpolation function on the intervals "
390 "determined by consecutive target values.") },
391 { GNM_FUNC_HELP_NOTE, F_("Strings and empty cells in @{abscissae} and @{ordinates} are ignored.") },
392 { GNM_FUNC_HELP_NOTE, F_("If several target data are provided they must be in the same column in consecutive cells.") },
393 { GNM_FUNC_HELP_EXAMPLES, "=interpolation(array(1,2,3),array(10,20,20),1.5,0)" },
394 { GNM_FUNC_HELP_EXAMPLES, "=interpolation(array(1,2,3),array(10,20,20),array(1.5,4),1)" },
395 { GNM_FUNC_HELP_SEEALSO, "PERIODOGRAM" },
396 { GNM_FUNC_HELP_END, NULL }
397 };
398
399 static GnmValue *
400 gnumeric_interpolation (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
401 {
402 gnm_float *vals0, *vals1, *vals2, *fres;
403 int n0, n2;
404 int interp;
405 GnmValue *error = NULL;
406 GnmValue *res;
407 CollectFlags flags;
408 GnmEvalPos const * const ep = ei->pos;
409 GnmValue const * const PtInterpolation = argv[2];
410 int r, i;
411 GSList *missing2 = NULL, *missing;
412 INTERPPROC interpproc = NULL;
413 gboolean constp = FALSE;
414
415 /* argv[2] */
416
417 int const cols = value_area_get_width (PtInterpolation, ep);
418 int const rows = value_area_get_height (PtInterpolation, ep);
419
420 if (rows == 0 || cols != 1) {
421 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
422 return res;
423 }
424
425 flags = COLLECT_IGNORE_BLANKS | COLLECT_IGNORE_STRINGS | COLLECT_IGNORE_BOOLS | COLLECT_IGNORE_ERRORS;
426 vals2 = collect_floats_value_with_info (PtInterpolation, ei->pos, flags,
427 &n2, &missing2, &error);
428 if (error) {
429 g_slist_free (missing2);
430 return error;
431 }
432
433 /* argv[3] */
434
435 if (argv[3]) {
436 interp = (int) gnm_floor (value_get_as_float (argv[3]));
437 if (interp < 0 || interp > INTERPOLATION_SPLINE_AVG) {
438 g_slist_free (missing2);
439 g_free (vals2);
440 return value_new_error_VALUE (ei->pos);
441 }
442 } else
443 interp = INTERPOLATION_LINEAR;
444
445 switch (interp) {
446 case INTERPOLATION_LINEAR:
447 interpproc = linear_interpolation;
448 break;
449 case INTERPOLATION_LINEAR_AVG:
450 interpproc = linear_averaging;
451 n2--;
452 break;
453 case INTERPOLATION_STAIRCASE:
454 interpproc = staircase_interpolation;
455 break;
456 case INTERPOLATION_STAIRCASE_AVG:
457 interpproc = staircase_averaging;
458 n2--;
459 break;
460 case INTERPOLATION_SPLINE:
461 interpproc = spline_interpolation;
462 break;
463 case INTERPOLATION_SPLINE_AVG:
464 interpproc = spline_averaging;
465 n2--;
466 break;
467 }
468
469 if (n2 <= 0) {
470 g_slist_free (missing2);
471 g_free (vals2);
472 return value_new_error_std (ei->pos, GNM_ERROR_VALUE);
473 }
474
475 /* argv[0] & argv[1] */
476
477 flags = COLLECT_IGNORE_BLANKS | COLLECT_IGNORE_STRINGS | COLLECT_IGNORE_BOOLS;
478 error = collect_float_pairs (argv[0], argv[1], ei->pos, flags, &vals0, &vals1,
479 &n0, &constp);
480
481 if (error) {
482 g_slist_free (missing2);
483 g_free (vals2);
484 return error;
485 }
486
487 /* Check whether the abscissae are increasing, if not order them */
488 if (!gnm_range_increasing (vals0, n0)) {
489 gboolean switched = TRUE;
490 if (constp) {
491 vals0 = g_memdup (vals0, sizeof(gnm_float) * n0);
492 vals1 = g_memdup (vals1, sizeof(gnm_float) * n0);
493 constp = FALSE;
494 }
495 while (switched) {
496 gnm_float *val;
497 switched = FALSE;
498 for (i = 1, val = vals0; i < n0; i++, val++) {
499 if (*val == *(val + 1)) {
500 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE) ;
501 goto done;
502 }
503 if (*val > *(val + 1)) {
504 gnm_float v = *val;
505 *val = *(val + 1);
506 *(val + 1) = v;
507 v = *(vals1 + i);
508 *(vals1 + i) = *(vals1 + i - 1);
509 *(vals1 + i - 1) = v;
510 switched = TRUE;
511 }
512 }
513 }
514 }
515
516 {
517 int n = n2;
518
519 if (missing2)
520 gnm_strip_missing (vals2, &n, missing2);
521 res = value_new_array_non_init (1 , n2);
522 i = 0;
523
524 res->v_array.vals[0] = g_new (GnmValue *, n2);
525
526 fres = interpproc (vals0, vals1, n0, vals2, n);
527 missing = missing2;
528 if (fres) {
529 i = 0;
530 for (r = 0 ; r < n2; ++r)
531 if (missing && r == GPOINTER_TO_INT (missing->data)) {
532 missing = missing->next;
533 res->v_array.vals[0][r] = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
534 } else
535 res->v_array.vals[0][r] = value_new_float (fres[i++]);
536 g_free (fres);
537 } else {
538 for( r = 0 ; r < n2; ++r)
539 res->v_array.vals[0][r] = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
540 }
541
542 }
543
544 done:
545 g_slist_free (missing2);
546 if (!constp) {
547 g_free (vals0);
548 g_free (vals1);
549 }
550 g_free (vals2);
551 return res;
552 }
553
554 /*********************************************************************************************************/
555
556 /******************************************************************************/
557 /* PERIODOGRAM WITH INTERPOLATION FUNCTION */
558 /******************************************************************************/
559
560 static GnmFuncHelp const help_periodogram[] = {
561 { GNM_FUNC_HELP_NAME, F_("PERIODOGRAM:periodogram of the given data") },
562 { GNM_FUNC_HELP_ARG, F_("ordinates:ordinates of the given data") },
563 { GNM_FUNC_HELP_ARG, F_("filter:windowing function to be used, defaults to no filter") },
564 { GNM_FUNC_HELP_ARG, F_("abscissae:abscissae of the given data, defaults to regularly spaced abscissae") },
565 { GNM_FUNC_HELP_ARG, F_("interpolation:method of interpolation, defaults to none") },
566 { GNM_FUNC_HELP_ARG, F_("number:number of interpolated data points") },
567 { GNM_FUNC_HELP_DESCRIPTION, F_("If an interpolation method is used, the number of returned values is one less than the number of targets and the targets values must be given in increasing order.") },
568 { GNM_FUNC_HELP_DESCRIPTION, F_("The output consists always of one column of numbers.") },
569 INTERPOLATIONMETHODS,
570 { GNM_FUNC_HELP_DESCRIPTION, F_("Possible window functions are:\n"
571 "0: no filter (rectangular window)\n"
572 "1: Bartlett (triangular window)\n"
573 "2: Hahn (cosine window)\n"
574 "3: Welch (parabolic window)") },
575 { GNM_FUNC_HELP_NOTE, F_("Strings and empty cells in @{abscissae} and @{ordinates} are ignored.") },
576 { GNM_FUNC_HELP_NOTE, F_("If several target data are provided they must be in the same column in consecutive cells.") },
577 { GNM_FUNC_HELP_SEEALSO, "INTERPOLATION" },
578 { GNM_FUNC_HELP_END, NULL }
579 };
580
581 static GnmValue *
582 gnumeric_periodogram (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
583 {
584 gnm_float *ord, *absc;
585 int filter, interp;
586 int n0, n1, nb;
587 GnmValue *error = NULL;
588 GnmValue *res;
589 CollectFlags flags;
590 GnmEvalPos const * const ep = ei->pos;
591 GnmValue const * const Pt = argv[0];
592 int i;
593 GSList *missing0 = NULL, *missing1 = NULL;
594 gnm_complex *in, *out = NULL;
595
596 int const cols = value_area_get_width (Pt, ep);
597 int const rows = value_area_get_height (Pt, ep);
598
599 if (cols == 1)
600 nb=rows;
601 else {
602 if (rows == 1)
603 nb=cols;
604 else
605 nb=0;
606 }
607
608 if (nb == 0) {
609 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
610 return res;
611 }
612
613 flags=COLLECT_IGNORE_BLANKS | COLLECT_IGNORE_STRINGS | COLLECT_IGNORE_BOOLS;
614
615 ord = collect_floats_value_with_info (argv[0], ei->pos, flags,
616 &n0, &missing0, &error);
617 if (error) {
618 g_slist_free (missing0);
619 return error;
620 }
621
622 if (n0 == 0) {
623 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
624 return res;
625 }
626
627 if (argv[1]) {
628 filter = (int) gnm_floor (value_get_as_float (argv[1]));
629 if (filter < 0 || filter > FILTER_WELCH) {
630 g_slist_free (missing0);
631 g_free (ord);
632 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
633 return res;
634 }
635 } else
636 filter = FILTER_NONE;
637
638 if (argv[2]) {
639 gnm_float *interpolated, *new_ord, start, incr;
640 int n2;
641 INTERPPROC(interpproc) = NULL;
642 absc = collect_floats_value_with_info (argv[2], ei->pos, flags,
643 &n1, &missing1, &error);
644 if (n1 == 1) {
645 g_slist_free (missing1);
646 g_free (absc);
647 goto no_absc;
648 }
649 if (error) {
650 g_slist_free (missing0);
651 g_slist_free (missing1);
652 g_free (absc);
653 return error;
654 }
655 if (n1 == 0) {
656 g_slist_free (missing0);
657 g_slist_free (missing1);
658 g_free (absc);
659 g_free (ord);
660 return value_new_error_std (ei->pos, GNM_ERROR_VALUE);
661 }
662 if (argv[3]) {
663 interp = (int) gnm_floor (value_get_as_float (argv[3]));
664 if (interp < 0 || interp > INTERPOLATION_SPLINE_AVG) {
665 g_slist_free (missing0);
666 g_slist_free (missing1);
667 g_free (absc);
668 g_free (ord);
669 return error;
670 }
671 } else
672 interp = INTERPOLATION_LINEAR;
673
674 if (missing0 || missing1) {
675 GSList *missing = gnm_slist_sort_merge (missing0, missing1);
676 gnm_strip_missing (ord, &n0, missing);
677 gnm_strip_missing (absc, &n1, missing);
678 g_slist_free (missing);
679
680 if (n0 != n1)
681 g_warning ("This should not happen. n0=%d n1=%d\n",
682 n0, n1);
683 }
684 n0 = n1 = MIN (n0, n1);
685 /* here we test if there is abscissas are always increasing, if not,
686 an error is returned */
687 if (n0 < 2 || !gnm_range_increasing (absc, n0) || n0 == 0) {
688 g_free (absc);
689 g_free (ord);
690 return value_new_error_std (ei->pos, GNM_ERROR_VALUE);
691 }
692 if (argv[4]) {
693 n1 = (int) gnm_floor (value_get_as_float (argv[4]));
694 if (n1 < n0) {
695 g_free (absc);
696 g_free (ord);
697 return value_new_error_std (ei->pos, GNM_ERROR_VALUE);
698 }
699 nb = 1;
700 while (nb < n1)
701 nb *= 2;
702 } else {
703 n1 = 1;
704 while (n1 < n0)
705 n1 *= 2;
706 nb = n1;
707 }
708 incr = (absc[n0 - 1] - absc[0]) / n1;
709 switch (interp) {
710 case INTERPOLATION_LINEAR:
711 interpproc = linear_interpolation;
712 start = absc[0];
713 n2 = n1;
714 break;
715 case INTERPOLATION_LINEAR_AVG:
716 interpproc = linear_averaging;
717 start = absc[0] - incr / 2.;
718 n2 = n1 + 1;
719 break;
720 case INTERPOLATION_STAIRCASE:
721 interpproc = staircase_interpolation;
722 start = absc[0];
723 n2 = n1;
724 break;
725 case INTERPOLATION_STAIRCASE_AVG:
726 interpproc = staircase_averaging;
727 start = absc[0] - incr / 2.;
728 n2 = n1 + 1;
729 break;
730 case INTERPOLATION_SPLINE:
731 interpproc = spline_interpolation;
732 start = absc[0];
733 n2 = n1;
734 break;
735 case INTERPOLATION_SPLINE_AVG:
736 interpproc = spline_averaging;
737 start = absc[0] - incr / 2.;
738 n2 = n1 + 1;
739 break;
740 default:
741 g_free (absc);
742 g_free (ord);
743 return value_new_error_std (ei->pos, GNM_ERROR_NA);
744 }
745 interpolated = g_new (gnm_float, n1);
746 for (i = 0; i < n2; i++)
747 interpolated[i] = start + i * incr;
748 new_ord = interpproc (absc, ord, n0, interpolated, n1);
749 g_free (ord);
750 ord = new_ord;
751 if (ord == NULL) {
752 g_free (absc);
753 g_free (interpolated);
754 return value_new_error_std (ei->pos, GNM_ERROR_NA);
755 }
756 n0 = n1;
757 } else {
758 no_absc:
759 /* we have no interpolation to apply, so just take the values */
760 if (missing0) {
761 gnm_strip_missing (ord, &n0, missing0);
762 g_slist_free (missing0);
763 }
764 nb = 1;
765 while (nb < n0)
766 nb *= 2;
767 }
768
769 /* Now apply the filter if any */
770 if (filter != FILTER_NONE) {
771 gnm_float factor;
772 switch (filter) {
773 case FILTER_BARTLETT:
774 factor = n0 / 2.;
775 for (i = 0; i < n0; i++)
776 ord[i] *= 1. - gnm_abs ((i / factor - 1));
777 break;
778 case FILTER_HANN:
779 factor = 2. * M_PIgnum / n0;
780 for (i = 0; i < n0; i++)
781 ord[i] *= 0.5 * (1 - gnm_cos (factor * i));
782 break;
783 case FILTER_WELCH:
784 factor = n0 / 2.;
785 for (i = 0; i < n0; i++)
786 ord[i] *= 1. - (i / factor - 1.) * (i / factor - 1.);
787 break;
788 }
789 }
790
791 /* Transform and return the result */
792 in = g_new0 (gnm_complex, nb);
793 for (i = 0; i < n0; i++){
794 in[i].re = ord[i];
795 }
796 g_free (ord);
797 gnm_fourier_fft (in, nb, 1, &out, FALSE);
798 g_free (in);
799 nb /= 2;
800 if (out && nb > 0) {
801 res = value_new_array_non_init (1 , nb);
802 res->v_array.vals[0] = g_new (GnmValue *, nb);
803 for (i = 0; i < nb; i++)
804 res->v_array.vals[0][i] =
805 value_new_float (gnm_sqrt (
806 out[i].re * out[i].re +
807 out[i].im * out[i].im));
808 g_free (out);
809 } else
810 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
811
812 return res;
813 }
814
815 /******************************************************************************/
816 /* Fourier Transform */
817 /******************************************************************************/
818
819 static GnmFuncHelp const help_fourier[] = {
820 { GNM_FUNC_HELP_NAME, F_("FOURIER:Fourier or inverse Fourier transform") },
821 { GNM_FUNC_HELP_ARG, F_("Sequence:the data sequence to be transformed") },
822 { GNM_FUNC_HELP_ARG, F_("Inverse:if true, the inverse Fourier transform is calculated, defaults to false") },
823 { GNM_FUNC_HELP_ARG, F_("Separate:if true, the real and imaginary parts are given separately, defaults to false") },
824 { GNM_FUNC_HELP_DESCRIPTION, F_("This array function returns the Fourier or inverse Fourier transform of the given data sequence.") },
825 { GNM_FUNC_HELP_DESCRIPTION, F_("The output consists of one column of complex numbers if @{Separate} is false and of two columns of real numbers if @{Separate} is true.") },
826 { GNM_FUNC_HELP_DESCRIPTION, F_("If @{Separate} is true the first output column contains the real parts and the second column the imaginary parts.") },
827 { GNM_FUNC_HELP_NOTE, F_("If @{Sequence} is neither an n by 1 nor 1 by n array, this function returns #VALUE!") },
828 { GNM_FUNC_HELP_END, NULL }
829 };
830
831 static GnmValue *
832 gnumeric_fourier (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
833 {
834 gnm_float *ord;
835 gboolean inverse = FALSE;
836 gboolean sep_columns = FALSE;
837 int n0, nb;
838 GnmValue *error = NULL;
839 GnmValue *res;
840 CollectFlags flags;
841 GnmEvalPos const * const ep = ei->pos;
842 GnmValue const * const Pt = argv[0];
843 int i;
844 GSList *missing0 = NULL;
845 gnm_complex *in, *out = NULL;
846
847 int const cols = value_area_get_width (Pt, ep);
848 int const rows = value_area_get_height (Pt, ep);
849
850 if (cols != 1 && rows != 1) {
851 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
852 return res;
853 }
854
855 flags=COLLECT_IGNORE_BLANKS | COLLECT_IGNORE_STRINGS | COLLECT_IGNORE_BOOLS;
856
857 ord = collect_floats_value_with_info (argv[0], ei->pos, flags,
858 &n0, &missing0, &error);
859 if (error) {
860 g_slist_free (missing0);
861 return error;
862 }
863
864 if (n0 == 0) {
865 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
866 return res;
867 }
868
869 if (argv[1]) {
870 inverse = 0 != (int) gnm_floor (value_get_as_float (argv[1]));
871 if (argv[2]) {
872 sep_columns = (0 != (int)
873 gnm_floor (value_get_as_float (argv[2])));
874 }
875 }
876
877 if (missing0) {
878 gnm_strip_missing (ord, &n0, missing0);
879 g_slist_free (missing0);
880 }
881
882 nb = 1;
883 while (nb < n0)
884 nb *= 2;
885
886 /* Transform and return the result */
887 in = g_new0 (gnm_complex, nb);
888 for (i = 0; i < n0; i++)
889 in[i].re = ord[i];
890 g_free (ord);
891 gnm_fourier_fft (in, nb, 1, &out, inverse);
892 g_free (in);
893
894 if (out && !sep_columns) {
895 res = value_new_array_empty (1 , nb);
896 for (i = 0; i < nb; i++)
897 res->v_array.vals[0][i] = value_new_string_nocopy
898 (gnm_complex_to_string (&(out[i]), 'i'));
899 g_free (out);
900 } else if (out && sep_columns) {
901 res = value_new_array_empty (2 , nb);
902 for (i = 0; i < nb; i++) {
903 res->v_array.vals[0][i] = value_new_float (out[i].re);
904 res->v_array.vals[1][i] = value_new_float (out[i].im);
905 }
906 g_free (out);
907 } else
908 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
909
910 return res;
911 }
912
913 /******************************************************************************/
914
915 static GnmFuncHelp const help_hpfilter[] = {
916 { GNM_FUNC_HELP_NAME, F_("HPFILTER:Hodrick Prescott Filter") },
917 { GNM_FUNC_HELP_ARG, F_("Sequence:the data sequence to be transformed") },
918 { GNM_FUNC_HELP_ARG, F_("\316\273:filter parameter \316\273, defaults to 1600") },
919 { GNM_FUNC_HELP_DESCRIPTION, F_("This array function returns the trend and cyclical components obtained by applying the Hodrick Prescott Filter with parameter @{\316\273} to the given data sequence.") },
920 { GNM_FUNC_HELP_DESCRIPTION, F_("The output consists of two columns of numbers, the first containing the trend component, the second the cyclical component.") },
921 { GNM_FUNC_HELP_NOTE, F_("If @{Sequence} is neither an n by 1 nor 1 by n array, this function returns #VALUE!") },
922 { GNM_FUNC_HELP_NOTE, F_("If @{Sequence} contains less than 6 numerical values, this function returns #VALUE!") },
923 { GNM_FUNC_HELP_END, NULL }
924 };
925
926
927 static void
928 gnm_hpfilter (gnm_float *data, int n, gnm_float lambda, int *err)
929 {
930 gnm_float *a, *b, *c;
931 int i;
932 gnm_float lambda6 = 6 * lambda + 1;
933 gnm_float lambda4 = -4 * lambda;
934 gnm_float h[5] = {0,0,0,0,0};
935 gnm_float g[5] = {0,0,0,0,0};
936 gnm_float j[2] = {0,0};
937 gnm_float h_b, h_c, denom;
938
939 g_return_if_fail (n > 5);
940 g_return_if_fail (data != NULL);
941 g_return_if_fail (err != NULL);
942
943 /* Initializing arrays a, b, and c */
944
945 a = g_new (gnm_float, n);
946 b = g_new (gnm_float, n);
947 c = g_new (gnm_float, n);
948
949 a[0] = lambda + 1;
950 b[0] = -2 * lambda;
951 c[0] = lambda;
952
953 for (i = 1; i < n - 2; i++) {
954 a[i] = lambda6;
955 b[i] = lambda4;
956 c[i] = lambda;
957 }
958
959 a[n - 2] = a[1] = lambda6 - lambda;
960 a[n - 1] = a[0];
961 b[n - 2] = b[0];
962 b[n - 1] = 0;
963 c[n - 2] = 0;
964 c[n - 1] = 0;
965
966 /* Forward */
967 for (i = 0; i < n; i++) {
968 denom = a[i]- h[3]*h[0] - g[4]*g[1];
969 if (denom == 0) {
970 *err = GNM_ERROR_DIV0;
971 goto done;
972 }
973
974 h_b = b[i];
975 g[0] = h[0];
976 b[i] = h[0] = (h_b - h[3] * h[1])/denom;
977
978 h_c = c[i];
979 g[1] = h[1];
980 c[i] = h[1] = h_c/denom;
981
982 a[i] = (data[i] - g[2]*g[4] - h[2]*h[3])/denom;
983
984 g[2] = h[2];
985 h[2] = a[i];
986 h[3] = h_b - h[4] * g[0];
987 g[4] = h[4];
988 h[4] = h_c;
989 }
990
991 data[n - 1] = j[0] = a[n - 1];
992
993 /* Backwards */
994 for (i = n - 1; i > 0; i--) {
995 data[i - 1] = a[i - 1] - b[i - 1] * j[0] - c[i - 1] * j[1];
996 j[1] = j[0];
997 j[0] = data[i - 1];
998 }
999
1000 done:
1001 g_free (a);
1002 g_free (b);
1003 g_free (c);
1004 return;
1005 }
1006
1007 static GnmValue *
1008 gnumeric_hpfilter (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
1009 {
1010 gnm_float *raw, *filtered;
1011 gnm_float lambda;
1012 int n = 0;
1013 GnmValue *error = NULL;
1014 GnmValue *res;
1015 CollectFlags flags;
1016 GnmEvalPos const * const ep = ei->pos;
1017 GnmValue const * const Pt = argv[0];
1018 int i, err = -1;
1019
1020 int const cols = value_area_get_width (Pt, ep);
1021 int const rows = value_area_get_height (Pt, ep);
1022
1023 if (cols != 1 && rows != 1) {
1024 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
1025 return res;
1026 }
1027
1028 flags=COLLECT_IGNORE_BLANKS | COLLECT_IGNORE_STRINGS | COLLECT_IGNORE_BOOLS;
1029
1030 raw = collect_floats_value (argv[0], ei->pos, flags,
1031 &n, &error);
1032 if (error)
1033 return error;
1034
1035 if (n < 6) {
1036 g_free (raw);
1037 res = value_new_error_std (ei->pos, GNM_ERROR_VALUE);
1038 return res;
1039 }
1040
1041 if (argv[1])
1042 lambda = value_get_as_float (argv[1]);
1043 else
1044 lambda = 1600.;
1045
1046
1047 /* Filter and return the result */
1048 filtered = g_new0 (gnm_float, n);
1049 for (i = 0; i < n; i++)
1050 filtered[i] = raw[i];
1051 gnm_hpfilter (filtered, n, lambda, &err);
1052 if (err > -1) {
1053 g_free (raw);
1054 g_free (filtered);
1055 res = value_new_error_std (ei->pos, err);
1056 return res;
1057 }
1058
1059 res = value_new_array_empty (2 , n);
1060 for (i = 0; i < n; i++) {
1061 res->v_array.vals[0][i] = value_new_float (filtered[i]);
1062 res->v_array.vals[1][i] = value_new_float (raw[i] - filtered[i]);
1063 }
1064 g_free (raw);
1065 g_free (filtered);
1066
1067 return res;
1068 }
1069
1070 const GnmFuncDescriptor TimeSeriesAnalysis_functions[] = {
1071
1072 { "interpolation", "AAA|f",
1073 help_interpolation, gnumeric_interpolation, NULL,
1074 GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
1075
1076 { "periodogram", "A|fAff",
1077 help_periodogram, gnumeric_periodogram, NULL,
1078 GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
1079
1080 { "fourier", "A|bb",
1081 help_fourier, gnumeric_fourier, NULL,
1082 GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
1083
1084 { "hpfilter", "A|f",
1085 help_hpfilter, gnumeric_hpfilter, NULL,
1086 GNM_FUNC_RETURNS_NON_SCALAR, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
1087
1088 {NULL}
1089 };
gnumeric mirror