1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1997--2019 The R Core Team
4 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
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_ext/Arith.h>
26 #include <R_ext/Error.h>
27 #include <R_ext/Applic.h>
28 #include <Rinternals.h> // for R_xlen_t
29 #ifdef DEBUG_approx
30 # include <R_ext/Print.h>
31 #endif
32
33 #ifdef ENABLE_NLS
34 #include <libintl.h>
35 #define _(String) dgettext ("stats", String)
36 #else
37 #define _(String) (String)
38 #endif
39
40 /* Linear and Step Function Interpolation */
41
42 /* Assumes that ordinates are in ascending order
43 * The right interval is found by bisection
44 * Linear/constant interpolation then takes place on that interval
45 */
46
47 /* NB: interv(.) in ../../../appl/interv.c
48 ------ is conceptually a special case of this, where y = 1:n */
49
50 typedef struct {
51 double ylow;
52 double yhigh;
53 double f1;
54 double f2;
55 int kind;
56 int na_rm;
57 } appr_meth;
58
approx1(double v,double * x,double * y,R_xlen_t n,appr_meth * Meth)59 static double approx1(double v, double *x, double *y, R_xlen_t n,
60 appr_meth *Meth)
61 {
62 /* Approximate y(v), given (x,y)[i], i = 0,..,n-1 */
63
64 #ifdef DEBUG_approx
65 REprintf("approx1(*, n = %.0f, Meth:=(f1=%g,f2=%g, kind=%d)):\n",
66 (double)n, Meth->f1, Meth->f2, Meth->kind);
67 #endif
68
69 if(!n) return R_NaN;
70
71 R_xlen_t
72 i = 0,
73 j = n - 1;
74 /* handle out-of-domain points */
75 if(v < x[i]) return Meth->ylow;
76 if(v > x[j]) return Meth->yhigh;
77
78
79 /* find the correct interval by bisection */
80 while(i < j - 1) { /* x[i] <= v <= x[j] */
81 R_xlen_t ij = (i+j) / 2;
82 /* i+1 <= ij <= j-1 */
83 if(v < x[ij]) j = ij; else i = ij;
84 /* still i < j */
85 #ifdef DEBUG_approx
86 REprintf(" (i,j) = (%.0f,%.0f)\n", (double)i, (double)j);
87 #endif
88 }
89 /* provably have i == j-1 */
90
91 /* interpolation */
92
93 if(v == x[j]) return y[j];
94 if(v == x[i]) return y[i];
95 /* impossible: if(x[j] == x[i]) return y[i]; */
96
97 if(Meth->kind == 1) /* linear */
98 return y[i] + (y[j] - y[i]) * ((v - x[i])/(x[j] - x[i]));
99 else /* 2 : constant */
100 return (Meth->f1 != 0.0 ? y[i] * Meth->f1 : 0.0)
101 + (Meth->f2 != 0.0 ? y[j] * Meth->f2 : 0.0);
102 }/* approx1() */
103
104
105 /* Testing done only once - in a separate function */
106 static void
R_approxtest(double * x,double * y,R_xlen_t nxy,int method,double f,int na_rm)107 R_approxtest(double *x, double *y, R_xlen_t nxy, int method, double f, int na_rm)
108 {
109 switch(method) {
110 case 1: /* linear */
111 break;
112 case 2: /* constant */
113 if(!R_FINITE(f) || f < 0 || f > 1)
114 error(_("approx(): invalid f value"));
115 break;
116 default:
117 error(_("approx(): invalid interpolation method"));
118 break;
119 }
120 /* check interpolation method */
121 if(na_rm) { // (x,y) should not have any NA's anymore
122 for(R_xlen_t i = 0; i < nxy; i++)
123 if(ISNAN(x[i]) || ISNAN(y[i]))
124 error(_("approx(): attempted to interpolate NA values"));
125 } else { // na.rm = FALSE ==> at least y may contain NA's
126 for(R_xlen_t i = 0; i < nxy; i++)
127 if(ISNAN(x[i]))
128 error(_("approx(x,y, .., na.rm=FALSE): NA values in x are not allowed"));
129 }
130 }
131
132 /* R Frontend for Linear and Constant Interpolation, no testing */
133
134 static void
R_approxfun(double * x,double * y,R_xlen_t nxy,double * xout,double * yout,R_xlen_t nout,int method,double yleft,double yright,double f,int na_rm)135 R_approxfun(double *x, double *y, R_xlen_t nxy, double *xout, double *yout,
136 R_xlen_t nout, int method, double yleft, double yright, double f, int na_rm)
137 {
138 appr_meth M = {0.0, 0.0, 0.0, 0.0, 0}; /* -Wall */
139
140 M.f2 = f;
141 M.f1 = 1 - f;
142 M.kind = method;
143 M.ylow = yleft;
144 M.yhigh = yright;
145 M.na_rm = na_rm;
146 #ifdef DEBUG_approx
147 REprintf("R_approxfun(x,y, nxy = %.0f, .., nout = %.0f, method = %d, ...)",
148 (double)nxy, (double)nout, Meth->kind);
149 #endif
150 for(R_xlen_t i = 0; i < nout; i++)
151 yout[i] = ISNAN(xout[i]) ? xout[i] : approx1(xout[i], x, y, nxy, &M);
152 }
153
154 #include <Rinternals.h>
155 #include "statsR.h"
ApproxTest(SEXP x,SEXP y,SEXP method,SEXP f,SEXP na_rm)156 SEXP ApproxTest(SEXP x, SEXP y, SEXP method, SEXP f, SEXP na_rm)
157 {
158 R_xlen_t nx = XLENGTH(x);
159 R_approxtest(REAL(x), REAL(y), nx,
160 asInteger(method), asReal(f), asLogical(na_rm));
161 // if no error was signalled,
162 return R_NilValue;
163 }
164
Approx(SEXP x,SEXP y,SEXP v,SEXP method,SEXP yleft,SEXP yright,SEXP f,SEXP na_rm)165 SEXP Approx(SEXP x, SEXP y, SEXP v, SEXP method,
166 SEXP yleft, SEXP yright, SEXP f, SEXP na_rm)
167 {
168 SEXP xout = PROTECT(coerceVector(v, REALSXP));
169 R_xlen_t nx = XLENGTH(x), nout = XLENGTH(xout);
170 SEXP yout = PROTECT(allocVector(REALSXP, nout));
171 R_approxfun(REAL(x), REAL(y), nx, REAL(xout), REAL(yout), nout,
172 asInteger(method), asReal(yleft), asReal(yright), asReal(f), asLogical(na_rm));
173 UNPROTECT(2);
174 return yout;
175 }
176