1#  File src/library/stats/R/approx.R
2#  Part of the R package, https://www.R-project.org
3#
4#  Copyright (C) 1995-2019 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#  A copy of the GNU General Public License is available at
17#  https://www.R-project.org/Licenses/
18
19### approx() and approxfun() are *very similar* -- keep in sync!
20
21## This function is used in approx, approxfun, spline, and splinefun
22## to massage the input (x,y) pairs into standard form:
23## x values unique and increasing, y values collapsed to match
24## (except if ties=="ordered", then not unique)
25regularize.values <- function(x, y, ties, warn.collapsing = TRUE, na.rm = TRUE) {
26    x <- xy.coords(x, y, setLab = FALSE) # -> (x,y) numeric of same length
27    y <- x$y
28    x <- x$x
29    keptNA <- FALSE
30    nx <-
31    if(any(na <- is.na(x) | is.na(y))) {
32	ok <- !na
33      if(na.rm) {
34	x <- x[ok]
35	y <- y[ok]
36        length(x)
37      } else { ## na.rm is FALSE
38          keptNA <- TRUE
39          sum(ok)
40      }
41    } else {
42        length(x)
43    }
44    if (!identical(ties, "ordered")) {
45	ordered <-
46	    if(is.function(ties) || is.character(ties))# fn or name of one
47		FALSE
48	    else if(is.list(ties) && length(T <- ties) == 2L && is.function(T[[2]])) {
49		## e.g. ties ==  list("ordered", mean)
50		ties <- T[[2]]
51		identical(T[[1]], "ordered")
52	    } else
53		stop("'ties' is not \"ordered\", a function, or list(<string>, <function>)")
54	if(!ordered && is.unsorted(if(keptNA) x[ok] else x)) {
55	    o <- order(x)
56	    x <- x[o]
57	    y <- y[o]
58	}
59	if (length(ux <- unique(x)) < nx) {
60	    if (warn.collapsing)
61		warning("collapsing to unique 'x' values")
62	    # tapply bases its uniqueness judgement on character representations;
63	    # we want to use values (PR#14377)
64	    y <- as.vector(tapply(y, match(x,x), ties))# as.v: drop dim & dimn.
65	    x <- ux
66	    stopifnot(length(y) == length(x))# (did happen in 2.9.0-2.11.x)
67            if(keptNA) ok <- !is.na(x)
68	}
69    }
70    list(x=x, y=y, keptNA=keptNA, notNA = if(keptNA) ok)
71}
72
73approx <- function(x, y = NULL, xout, method = "linear", n = 50,
74		   yleft, yright, rule = 1, f = 0, ties = mean, na.rm = TRUE)
75{
76    method <- pmatch(method, c("linear", "constant"))
77    if (is.na(method)) stop("invalid interpolation method")
78    stopifnot(is.numeric(rule), (lenR <- length(rule)) >= 1L, lenR <= 2L)
79    if(lenR == 1) rule <- rule[c(1,1)]
80    r <- regularize.values(x, y, ties, missing(ties), na.rm=na.rm)
81                                        # -> (x,y) numeric of same length
82    y <- r$y
83    x <- r$x
84    noNA <- na.rm || !r$keptNA
85    nx <- if(noNA)
86              length(x) # large vectors ==> non-integer
87          else
88              sum(r$notNA)
89    if (is.na(nx)) stop("invalid length(x)")
90    if (nx <= 1) {
91	if(method == 1)# linear
92	    stop("need at least two non-NA values to interpolate")
93	if(nx == 0) stop("zero non-NA points")
94    }
95
96    if (missing(yleft))
97	yleft <- if (rule[1L] == 1) NA else y[1L]
98    if (missing(yright))
99	yright <- if (rule[2L] == 1) NA else y[length(y)]
100    stopifnot(length(yleft) == 1L, length(yright) == 1L, length(f) == 1L)
101    if (missing(xout)) {
102	if (n <= 0) stop("'approx' requires n >= 1")
103        xout <-
104            if(noNA)
105                seq.int(x[1L], x[nx], length.out = n)
106            else {
107                xout <- x[r$notNA]
108                seq.int(xout[1L], xout[length(xout)], length.out = n)
109            }
110    }
111    x <- as.double(x); y <- as.double(y)
112    .Call(C_ApproxTest, x, y, method, f, na.rm)
113    yout <- .Call(C_Approx, x, y, xout, method, yleft, yright, f, na.rm)
114    list(x = xout, y = yout)
115}
116
117approxfun <- function(x, y = NULL, method = "linear",
118		   yleft, yright, rule = 1, f = 0, ties = mean, na.rm = TRUE)
119{
120    method <- pmatch(method, c("linear", "constant"))
121    if (is.na(method)) stop("invalid interpolation method")
122    stopifnot(is.numeric(rule), (lenR <- length(rule)) >= 1L, lenR <= 2L)
123    if(lenR == 1) rule <- rule[c(1,1)]
124    x <- regularize.values(x, y, ties, missing(ties), na.rm=na.rm)
125                                        # -> (x,y) numeric of same length
126    nx <- if(na.rm || !x$keptNA)
127              length(x$x) # large vectors ==> non-integer
128          else
129              sum(x$notNA)
130    if (is.na(nx)) stop("invalid length(x)")
131    if (nx <= 1) {
132	if(method == 1)# linear
133	    stop("need at least two non-NA values to interpolate")
134	if(nx == 0) stop("zero non-NA points")
135    }
136    y <- x$y
137    if (missing(yleft))
138	yleft <- if (rule[1L] == 1) NA else y[1L]
139    if (missing(yright))
140	yright <- if (rule[2L] == 1) NA else y[length(y)]
141    stopifnot(length(yleft) == 1L, length(yright) == 1L, length(f) == 1L)
142    rm(rule, ties, lenR, nx) # we do not need nx, but summary.stepfun did.
143
144    ## 1. Test input consistency once
145    x <- as.double(x$x); y <- as.double(y)
146    .Call(C_ApproxTest, x, y, method, f, na.rm)
147
148    ## 2. Create and return function that does not test input validity...
149    function(v) .approxfun(x, y, v, method, yleft, yright, f, na.rm)
150}
151
152## avoid capturing internal calls
153## default for 'na.rm': for old saved approxfun() {incl ecdf()} results
154.approxfun <- function(x, y, v,  method, yleft, yright, f, na.rm=TRUE)
155    .Call(C_Approx, x, y, v, method, yleft, yright, f, na.rm)
156