-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathREADME.Rmd
133 lines (98 loc) · 5.55 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 5,
fig.height = 3,
fig.align='center',
fig.path = "man/figures/README-"
)
```
# itp
[](https://ci.appveyor.com/project/paulnorthrop/itp)
[](https://github.com/paulnorthrop/itp/actions/workflows/R-CMD-check.yaml)
[](https://app.codecov.io/github/paulnorthrop/itp?branch=main)
[](https://cran.r-project.org/package=itp)
[](https://cran.r-project.org/package=itp)
[](https://cran.r-project.org/package=itp)
## The Interpolate, Truncate, Project (ITP) Root-Finding Algorithm
The **itp** package implements the Interpolate, Truncate, Project (ITP) root-finding algorithm of [Oliveira and Takahashi (2021)](https://doi.org/10.1145/3423597). Each iteration of the algorithm results in a bracketing interval for the root that is narrower than the previous interval. It's performance compares favourably with existing methods on both well-behaved functions and ill-behaved functions while retaining the worst-case reliability of the bisection method. For details see the authors' [Kudos summary](https://www.growkudos.com/publications/10.1145%25252F3423597/reader) and the Wikipedia article [ITP method](https://en.wikipedia.org/wiki/ITP_method).
### Examples
We use three examples from Section 3 of [Oliveira and Takahashi (2021)](https://doi.org/10.1145/3423597) to illustrate the use of the `itp` function. Each of these functions has a root in the interval $(-1, 1)$. The function can be supplied either as an R function or as an external pointer to a C++ function.
```{r}
library(itp)
```
#### A continuous function
The Lambert function $l(x) = xe^x - 1$ is continuous.
```{r lambert, echo = FALSE}
oldpar <- par(mar = c(4.5, 4, 1, 1))
lambert <- function(x) x * exp(x) - 1
curve(lambert, -1, 1, main = "Lambert")
abline(h = 0, lty = 2)
abline(v = itp(lambert, c(-1, 1))$root, lty = 2)
par(oldpar)
```
The `itp` function finds an estimate of the root, that is, $x^{\ast}$ for which $f(x^{\ast})$ is (approximately) equal to 0. The algorithm continues until the length of the interval that brackets the root is smaller than $2 \epsilon$, where $\epsilon$ is a user-supplied tolerance. The default is $\epsilon = 10^{-10}$.
First, we supply an R function that evaluates the Lambert function.
```{r lambert_root_r}
# Lambert, using an R function
lambert <- function(x) x * exp(x) - 1
itp(lambert, c(-1, 1))
```
Now, we create an external pointer to a C++ function that has been provided in the `itp` package and pass this pointer to the function `itp()`. For more information see the [Overview of the itp package](https://paulnorthrop.github.io/itp/articles/itp-vignette.html) vignette.
```{r lambert_root_Cpp}
# Lambert, using an external pointer to a C++ function
lambert_ptr <- xptr_create("lambert")
itp(lambert_ptr, c(-1, 1))
```
#### The function `itp_c`
Also provided is the function `itp_c`, which is equivalent to `itp`, but the calculations are performed entirely using C++, and the arguments differ slightly: `itp_c` has a named required argument `pars` rather than `...` and it does not have the arguments `interval`, `f.a` or `f.b`.
```{r lambert_root_itp_c}
# Calling itp_c()
res <- itp_c(lambert_ptr, pars = list(), a = -1, b = 1)
res
```
#### A discontinuous function
The staircase function $s(x) = \lceil 10 x - 1 \rceil + 1/2$ is discontinuous.
```{r staircase, echo = FALSE}
oldpar <- par(mar = c(4.5, 4, 1, 1))
staircase <- function(x) ceiling(10 * x - 1) + 1 / 2
curve(staircase, -1, 1, main = "Staircase", n = 10000)
abline(h = 0, lty = 2)
abline(v = itp(staircase, c(-1, 1))$root, lty = 2)
par(oldpar)
```
The `itp` function finds the discontinuity at $x = 0$ at which the sign of the function changes. The value of 0.5 returned for the root `res$root` is the midpoint of the bracketing interval `[res$a, res$b]` at convergence.
```{r staircase_root}
# Staircase
staircase <- function(x) ceiling(10 * x - 1) + 1 / 2
res <- itp(staircase, c(-1, 1))
print(res, all = TRUE)
```
#### A function with multiple roots
The Warsaw function $w(x) = I(x > -1)\left(1 + \sin\left(\frac{1}{1+x}\right)\right)-1$ has multiple roots.
```{r warsaw, echo = FALSE}
oldpar <- par(mar = c(4.5, 4, 1, 1))
warsaw <- function(x) ifelse(x > -1, sin(1 / (x + 1)), -1)
curve(warsaw, -0.9999999, 1, main = "Warsaw", n = 1000)
abline(h = 0, lty = 2)
abline(v = itp(warsaw, c(-1, 1))$root, lty = 2)
par(oldpar)
```
When the initial interval is $[-1, 1]$ the `itp` function finds the root $x \approx -0.6817$. There are other roots that could be found from a different initial interval.
```{r warsaw_root}
# Warsaw
warsaw <- function(x) ifelse(x > -1, sin(1 / (x + 1)), -1)
itp(warsaw, c(-1, 1))
```
### Installation
To get the current released version from CRAN:
```{r installation, eval = FALSE}
install.packages("itp")
```
### Vignette
See the [Overview of the itp package](https://paulnorthrop.github.io/itp/articles/itp-vignette.html) vignette, which can also be accessed using `vignette("itp-vignette", package = "itp")`.