1 /*
2 * R : A Computer Language for Statistical Data Analysis
3
4 * Copyright (C) 1999-2016 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/.
19 */
20
21 #ifdef HAVE_CONFIG_H
22 # include <config.h>
23 #endif
24
25 #include <R.h>
26 #include "ts.h"
27
28 #ifndef min
29 #define min(a, b) ((a < b)?(a):(b))
30 #define max(a, b) ((a < b)?(b):(a))
31 #endif
32
33 // currently ISNAN includes NAs
34 #define my_isok(x) (!ISNA(x) & !ISNAN(x))
35
cfilter(SEXP sx,SEXP sfilter,SEXP ssides,SEXP scircular)36 SEXP cfilter(SEXP sx, SEXP sfilter, SEXP ssides, SEXP scircular)
37 {
38 if (TYPEOF(sx) != REALSXP || TYPEOF(sfilter) != REALSXP)
39 error("invalid input");
40 R_xlen_t nx = XLENGTH(sx), nf = XLENGTH(sfilter);
41 int sides = asInteger(ssides), circular = asLogical(scircular);
42 if(sides == NA_INTEGER || circular == NA_LOGICAL) error("invalid input");
43
44 SEXP ans = allocVector(REALSXP, nx);
45
46 R_xlen_t i, j, nshift;
47 double z, tmp, *x = REAL(sx), *filter = REAL(sfilter), *out = REAL(ans);
48
49 if(sides == 2) nshift = nf /2; else nshift = 0;
50 if(!circular) {
51 for(i = 0; i < nx; i++) {
52 z = 0;
53 if(i + nshift - (nf - 1) < 0 || i + nshift >= nx) {
54 out[i] = NA_REAL;
55 continue;
56 }
57 for(j = max(0, nshift + i - nx); j < min(nf, i + nshift + 1) ; j++) {
58 tmp = x[i + nshift - j];
59 if(my_isok(tmp)) z += filter[j] * tmp;
60 else { out[i] = NA_REAL; goto bad; }
61 }
62 out[i] = z;
63 bad:
64 continue;
65 }
66 } else { /* circular */
67 for(i = 0; i < nx; i++)
68 {
69 z = 0;
70 for(j = 0; j < nf; j++) {
71 R_xlen_t ii = i + nshift - j;
72 if(ii < 0) ii += nx;
73 if(ii >= nx) ii -= nx;
74 tmp = x[ii];
75 if(my_isok(tmp)) z += filter[j] * tmp;
76 else { out[i] = NA_REAL; goto bad2; }
77 }
78 out[i] = z;
79 bad2:
80 continue;
81 }
82 }
83 return ans;
84 }
85
86 /* recursive filtering */
rfilter(SEXP x,SEXP filter,SEXP out)87 SEXP rfilter(SEXP x, SEXP filter, SEXP out)
88 {
89 if (TYPEOF(x) != REALSXP || TYPEOF(filter) != REALSXP
90 || TYPEOF(out) != REALSXP) error("invalid input");
91 R_xlen_t nx = XLENGTH(x), nf = XLENGTH(filter);
92 double sum, tmp, *r = REAL(out), *rx = REAL(x), *rf = REAL(filter);
93
94 for(R_xlen_t i = 0; i < nx; i++) {
95 sum = rx[i];
96 for (R_xlen_t j = 0; j < nf; j++) {
97 tmp = r[nf + i - j - 1];
98 if(my_isok(tmp)) sum += tmp * rf[j];
99 else { r[nf + i] = NA_REAL; goto bad3; }
100 }
101 r[nf + i] = sum;
102 bad3:
103 continue;
104 }
105 return out;
106 }
107
108 /* now allows missing values */
109 static void
acf0(double * x,int n,int ns,int nl,Rboolean correlation,double * acf)110 acf0(double *x, int n, int ns, int nl, Rboolean correlation, double *acf)
111 {
112 int d1 = nl+1, d2 = ns*d1;
113
114 for(int u = 0; u < ns; u++)
115 for(int v = 0; v < ns; v++)
116 for(int lag = 0; lag <= nl; lag++) {
117 double sum = 0.0; int nu = 0;
118 for(int i = 0; i < n-lag; i++)
119 if(!ISNAN(x[i + lag + n*u]) && !ISNAN(x[i + n*v])) {
120 nu++;
121 sum += x[i + lag + n*u] * x[i + n*v];
122 }
123 acf[lag + d1*u + d2*v] = (nu > 0) ? sum/(nu + lag) : NA_REAL;
124 }
125 if(correlation) {
126 if(n == 1) {
127 for(int u = 0; u < ns; u++)
128 acf[0 + d1*u + d2*u] = 1.0;
129 } else {
130 double *se = (double *) R_alloc(ns, sizeof(double));
131 for(int u = 0; u < ns; u++)
132 se[u] = sqrt(acf[0 + d1*u + d2*u]);
133 for(int u = 0; u < ns; u++)
134 for(int v = 0; v < ns; v++)
135 for(int lag = 0; lag <= nl; lag++) { // ensure correlations remain in [-1,1] :
136 double a = acf[lag + d1*u + d2*v] / (se[u]*se[v]);
137 acf[lag + d1*u + d2*v] = (a > 1.) ? 1. : ((a < -1.) ? -1. : a);
138 }
139 }
140 }
141 }
142
acf(SEXP x,SEXP lmax,SEXP sCor)143 SEXP acf(SEXP x, SEXP lmax, SEXP sCor)
144 {
145 int nx = nrows(x), ns = ncols(x), lagmax = asInteger(lmax),
146 cor = asLogical(sCor);
147 x = PROTECT(coerceVector(x, REALSXP));
148 SEXP ans = PROTECT(allocVector(REALSXP, (lagmax + 1)*ns*ns));
149 acf0(REAL(x), nx, ns, lagmax, cor, REAL(ans));
150 SEXP d = PROTECT(allocVector(INTSXP, 3));
151 INTEGER(d)[0] = lagmax + 1;
152 INTEGER(d)[1] = INTEGER(d)[2] = ns;
153 setAttrib(ans, R_DimSymbol, d);
154 UNPROTECT(3);
155 return ans;
156 }
157