--- title: "Overview of the itp package" author: "Paul Northrop" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Overview of the itp package} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: itp.bib csl: taylor-and-francis-chicago-author-date.csl --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 5, fig.height = 3, fig.align='center', global.par = TRUE ) ``` The Interpolate, Truncate, Project (ITP) root-finding algorithm was developed in @OliveiraTakahashi2021. 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). The `itp` function implements the ITP method to find a root $x^*$ of the function $f: \mathbb{R} \rightarrow \mathbb{R}$ in the interval $[a, b]$, where $f(a)f(b) < 0$. If $f$ is continuous over $[a, b]$ then $f(x^*) = 0$. If $f$ is discontinuous over $[a, b]$ then $x^*$ may be an estimate of a point of discontinuity at which the sign of $f$ changes, that is, $f(x^* - \delta)f(x^* + \delta) \leq 0$, where $0 \leq \delta \leq \epsilon$ for some tolerance value $\epsilon$. We use some of the examples presented in Table 1 of @OliveiraTakahashi2021 to illustrate the use of this function and run the examples using the `uniroot` function in the `stats` package as a means of comparison, using a convergence tolerance of $10^{-10}$ in both cases. The `itp` function uses the following default values of the tuning parameters: $\kappa_1 = 0.2 / (b - a)$, $\kappa_2 = 2$ and $n_0 = 1$, but these may be changed by the user. See Sections 2 and 3 of @OliveiraTakahashi2021 for information. The following function prints output from `uniroot` in the same style as the output from `itp`. ```{r} # Method to print part of uniroot output print.list <- function(x, digits = max(3L, getOption("digits") - 3L)) { names(x)[1:3] <- c("root", "f(root)", "iterations") print.default(format(x[1:3], digits = digits), print.gap = 2L, quote = FALSE) } ``` ```{r setup} library(itp) ``` First, we consider supplying the argument `f` to `itp` using an R function. Then we show how to supply an external pointer to a C++ function. For the simple examples given in the `itp` package only a modest improvement in speed is observed (and expected). However, being able to call `itp()` on the C++ side may have benefits in more challenging problems. ## Supplying to `itp` an R function ### Well-behaved functions These functions are infinitely differentiable and contain only one simple root over $[−1, 1]$. #### Lambert: $f(x) = x e^x - 1$ ```{r} # Lambert lambert <- function(x) x * exp(x) - 1 itp(lambert, c(-1, 1)) uniroot(lambert, c(-1, 1), tol = 1e-10) ``` ```{r, echo = FALSE} oldpar <- par(mar = c(4, 4, 1, 1)) curve(lambert, -1, 1, main = "Lambert") abline(h = 0, lty = 2) abline(v = itp(lambert, c(-1, 1))$root, lty = 2) par(oldpar) ``` #### Trigonometric 1: $f(x) = \tan(x - 1 / 10)$ ```{r} # Trigonometric 1 trig1 <- function(x) tan(x - 1 /10) itp(trig1, c(-1, 1)) uniroot(trig1, c(-1, 1), tol = 1e-10) ``` ```{r, echo = FALSE} oldpar <- par(mar = c(4, 4, 1, 1)) curve(trig1, -1, 1, main = "Trigonometric 1") abline(h = 0, lty = 2) abline(v = itp(trig1, c(-1, 1))$root, lty = 2) par(oldpar) ``` ### Ill-behaved functions #### Polynomial 3 (non-simple root): $f(x) = (10^6 x - 1) ^ 3$ This function has a non-simple root at $10^{-6}$, with a multiplicity of 3. ```{r} # Polynomial 3 poly3 <- function(x) (x * 1e6 - 1) ^ 3 itp(poly3, c(-1, 1)) # Using n0 = 0 leads to (slightly) fewer iterations, in this example poly3 <- function(x) (x * 1e6 - 1) ^ 3 itp(poly3, c(-1, 1), n0 = 0) uniroot(poly3, c(-1, 1), tol = 1e-10) ``` ```{r, echo = FALSE} oldpar <- par(mar = c(4, 4, 1, 1)) curve(poly3, -1, 1, main = "Polynomial 3") abline(h = 0, lty = 2) abline(v = itp(poly3, c(-1, 1))$root, lty = 2) par(oldpar) ``` #### Staircase (discontinuous): $f(x) = \lceil 10 x - 1 \rceil + 1/2$ This function has discontinuities, including one at the location of the root. ```{r} # Staircase staircase <- function(x) ceiling(10 * x - 1) + 1 / 2 itp(staircase, c(-1, 1)) uniroot(staircase, c(-1, 1), tol = 1e-10) ``` ```{r, echo = FALSE} oldpar <- par(mar = c(4, 4, 1, 1)) 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) ``` #### Warsaw (multiple roots): $f(x) = I(x > -1)\left(1 + \sin\left(\frac{1}{1+x}\right)\right)-1$ This function has multiple roots: we find two of them. ```{r} # Warsaw warsaw <- function(x) ifelse(x > -1, sin(1 / (x + 1)), -1) # Function increasing over the interval itp(warsaw, c(-1, 1)) uniroot(warsaw, c(-1, 1), tol = 1e-10) # Function decreasing over the interval itp(warsaw, c(-0.85, -0.8)) uniroot(warsaw, c(-0.85, -0.8), tol = 1e-10) ``` ```{r, echo = FALSE} oldpar <- par(mar = c(4, 4, 1, 1)) curve(warsaw, -1, 1, main = "Warsaw", n = 1000) abline(h = 0, lty = 2) abline(v = itp(warsaw, c(-1, 1))$root, lty = 2) abline(v = itp(warsaw, c(-0.85, -0.8))$root, lty = 2) par(oldpar) ``` In terms of a naive comparison based on the number of iterations `itp` and `uniroot` perform similarly, except in the repeated-root "Polynomial 3" example, where `itp` requires fewer iterations. ## Supplying to `itp` an external pointer to a C++ function The general approach follows the article [Passing user-supplied C++ functions](https://gallery.rcpp.org/articles/passing-cpp-function-pointers/) in the [Rcpp Gallery](https://gallery.rcpp.org/). The user writes a C++ function to calculate $f$. This function must have a particular structure. As an example consider the following function, which implements the Lambert function considered above. double lambert_cpp(const double& x, const List& pars) { return x * exp(x) - 1.0 ; } The function returns a value of double type and has two arguments: the first is the main argument and is of double type and the second is a list containing the values of additional parameters whose values are not specified inside the function. This list must be present, even if an empty list will be passed to the function. This allows the user to change the values of any parameters in $f$ without editing the function. One way to provide C++ functions is to create them in a file, say `user_fns.cpp`. Example content is provided below. #include using namespace Rcpp; // [[Rcpp::interfaces(r, cpp)]] // User-supplied C++ functions for f. // The only interface is double fun(const double& x, const List& pars). // The second (List) argument must be included even if the function has no // additional arguments. // Each function must be prefaced by the line: // [[Rcpp::export]] // [[Rcpp::export]] double lambert_cpp(const double& x, const List& pars) { return x * exp(x) - 1.0 ; } // [[Rcpp::export]] SEXP xptr_create(std::string fstr) { typedef double (*funcPtr)(const double& x, const List& pars) ; if (fstr == "lambert") return(XPtr(new funcPtr(&lambert_cpp))) ; else return(XPtr(R_NilValue)) ; } The full file is available on the [itp Github page](https://raw.githubusercontent.com/paulnorthrop/itp/main/src/user_fns.cpp). The functions in this file are compiled and made available to R, either using the `Rcpp::sourceCpp` function (e.g. `Rcpp::sourceCpp("user_fns.cpp")`) or using RStudio's Source button on the editor toolbar. The example content below also includes the function `xptr_create`, which creates an external pointer to a C++ function. It is this external pointer that is passed to `itp`. If the user has written a C++ function, say `new_name`, then they need to add to `xptr_create` (or their own version of this function) two lines of code: else if (fstr == "new_name") return(Rcpp::XPtr(new funcPtr(&new_name))) ; to create an external pointer for `new_name` using `create_xptr`. ### Lambert: $f(x) = x e^x - 1$ We repeat the Lambert example, obtaining the same results as above when we used an R function to supply the function. ```{r lambert_root_Cpp} # Lambert, using an external pointer to a C++ function lambert_ptr <- xptr_create("lambert") res <- itp(lambert_ptr, c(-1, 1)) res ``` ## C++ 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`. This may be useful if you wish to call an ITP function on the C++ side. ```{r lambert_root_itp_c} # Calling itp_c() res <- itp_c(lambert_ptr, pars = list(), a = -1, b = 1) res ``` ## Plot method Objects of class `itp` have a plot method that, by default, produces a plot of the function over the interval over which the root was sought and indicates the location of the root using dashed lines. ```{r plot_itp, fig.show='hold'} oldpar <- par(mar = c(4, 4, 1, 1)) plot(res, main = "Lambert") par(oldpar) ``` ## References