# R/rational-function.R In weiyaw/blackbox: A SMC-SA Algorithm for Constrained Optimisation

```## ALL CODES FOR FITTING LEAST SQUARES AND TUKEY BIWEIGHT RATIONAL FUNCTIONS

#' Get the numerators of the derivatives of a set of rational functions.
#'
#' Assuming that the rational function is r(x) = f(x) / g(x), this function
#' returns f'(x)g(x) - f(x)g'(x).
#'
#' @param fs a matrix of numerator polynomial in ascending order representing,
#'     each column corresponds to a polynomial. (num mat)
#' @param gs a matrix of denominator coefficients in ascending order, each
#'     column corresponds to a polynomial. (num mat)
#' @param EPS error tolerance. (num)
#'
#' @return a matrix of polynomial coefficients.
get_deriv_nums <- function(fs, gs, EPS = 1e-6) {

lg <- NROW(gs)
lf <- NROW(fs)

if (lg == 1) {
## if the denominator is a constant
stop('This is a polynomial')

} else if (lf == 1) {
## if the numerator is a constant
## derivatives of denominator (no reversing)
gs_deriv <- gs[seq.int(2, lg), , drop = FALSE] * seq_len(lg - 1)
## result
res <- matrix(NA, nrow = lg - 1, ncol = NCOL(fs))
for (i in seq_len(NROW(gs_deriv))) {
res[i, ] <- -fs * gs_deriv[i, ]
}

} else {
## if the numerator is a polynomial
## derivatives of denominator (reversed)
gs_deriv_r <- gs[seq.int(lg, 2), , drop = FALSE] *
(lg - seq_len(lg - 1))
## derivatives of numerator (reversed)
fs_deriv_r <- fs[seq.int(lf, 2), , drop = FALSE] *
(lf - seq_len(lf - 1))
## result
res <- matrix(NA, nrow = lg + lf - 2, ncol = NCOL(fs))
for (i in seq_len(NCOL(fs))) {
res[, i] <- convolve(gs[, i], fs_deriv_r[, i], type = 'o') -
convolve(fs[, i], gs_deriv_r[, i], type = 'o')
}
}

## FAST IMPLEMENTATION, MIGHT NOT WORK IN CERTAIN CASES.
## This operation only remove the zero when differentiating rational
## function with same degrees of numerator and denominator.
if (lg == lf) {
return(res[-NROW(res), , drop = FALSE])
} else {
return(res)
}
}

#' Evaluate polynomials.
#'
#' The polynomials are evaluated at given points using the Horner scheme.
#'
#' @param poly_coefs a matrix of polynomial coefficients in ascending order,
#'     each column corresponds to a polynomial. (num mat)
#' @param x a vector of position where the polynomials are evaluated at. (num vec)
#'
#' @return a matrix, each row and column corresponds to a point in \code{x} and
#'     a polynomial respectively.
eval_pols <- function(poly_coefs, x) {

ncoefs <- NCOL(poly_coefs)
nx <- length(x)
res <- matrix(0, ncoefs, nx)
X <- matrix(rep(x, each = ncoefs), ncoefs, nx)
for (bi in seq.int(NROW(poly_coefs), 1)) {
res <- X * res + poly_coefs[bi, ]
}
t(res)
}

#' Are these polynomials positive?
#'
#' Check if a set of polynomial is positive (or negative) inside a given
#' interval.
#'
#' This can be done by checking wherether the roots of the polynomials have even
#' multiplicities. In addition, one of the points in the interval must have the
#' desired sign.
#'
#' @param poly_coefs a matrix of polynomial coefficients in ascending order,
#'     each column corresponds to a polynomial. (num mat)
#' @param a the lower bound. (num)
#' @param b the upper bound. (num)
#' @param positive if it is \code{FALSE}, check if the polynomial is
#'     negative. (bool)
#' @param EPS error tolerance. (num)
#'
#' @return a boolean vector indicating whether the condition in the description
#'     is met.
is_positives <- function(poly_coefs, a, b, positive = TRUE, EPS = 1e-06) {
if (a >= b) {
stop("a must be strictly larger than b.")
}

## Choose five points between the limits
if (is.infinite(a) && is.infinite(b)) {
chkpoints <- 1:5
} else if (is.infinite(a)) {
chkpoints <- b - 1:5
} else if (is.infinite(b)) {
chkpoints <- a + 1:5
} else {
chkpoints <- seq.int(a, b, length.out = 5)
}

## Check the signs of these points
if (positive) {
signs_mat <- eval_pols(poly_coefs, chkpoints) >= -EPS
} else {
signs_mat <- eval_pols(poly_coefs, chkpoints) <= EPS
}

## All points should have the same sign
satisfy <- colSums(signs_mat) == 5

## Check the multiplicities of the roots
for (i in which(satisfy)) {
roots <- polyroot(poly_coefs[, i])
re.roots <- Re(roots)[abs(Im(roots)) < EPS]
re.roots <- re.roots[a + EPS <= re.roots & re.roots <= b - EPS]
if (length(re.roots) == 0) {
satisfy[i] <- TRUE
} else {
mltplcty <- rowSums(outer(re.roots, re.roots,
function(x, y) abs(x - y) < EPS))
satisfy[i] <- all(mltplcty%%2 == 0L)
}
}
satisfy
}

#' Is there any roots in these polynomials?
#'
#' For a given set of polynomials, find all roots in a given interval and return
#' \code{TRUE} if there are no roots in it.
#'
#' @param poly_coefs a matrix of polynomial coefficients in ascending order,
#'     each column corresponds to a polynomial. (num mat)
#' @param a the lower bound. (num)
#' @param b the upper bound. (num)
#' @param EPS error tolerance. (num)
#'
#' @return a boolean vector indicating whether the condition in the description
#'     is met.
is_no_roots <- function(poly_coefs, a, b, EPS = 1e-06) {
satisfy <- rep(NA, NCOL(poly_coefs))
for (i in seq_len(NCOL(poly_coefs))) {
roots <- polyroot(poly_coefs[, i])
re.roots <- Re(roots)[abs(Im(roots)) < EPS]
re.roots <- re.roots[a - EPS <= re.roots & re.roots <= b + EPS]
satisfy[i] <- length(re.roots) == 0
}
satisfy
}

#' Is there any roots in these polynomials? (factory)
#'
#' Produce a function that checks whether a given set of rational functions is
#' monotone in a given interval.
#'
#' @param dim_f the number of coefficients in the numerator. (num)
#' @param dim_g the number of coefficients in the denominator, excluding the
#'     constant term. (num)
#' @param xmin the lower bound. (num)
#' @param xmax the upper bound. (num)
#' @param increasing set to \code{TRUE} when checking for monotonic increasing,
#'     \code{FALSE} for monotonic decreasing.
#' @param EPS error tolerance. (num)
#'
#' @return a function that takes in a matrix of rational function coefficients
#'     (numerator precede denominator, in ascending order), each column
#'     corresponds to a rational function, and returns a boolean vector
#'     indicating whether the condition in the description is met.
is_monotones_fac <- function(dim_f, dim_g, xmin = 0, xmax = Inf, increasing = TRUE,
EPS = 1e-06) {
force(xmin)
force(xmax)
force(increasing)
force(EPS)
idx_f <- seq.int(1, length.out = dim_f)
idx_g <- seq.int(dim_f + 1, length.out = dim_g)
n_dims <- dim_f + dim_g

function(rat_coefs) {
if (NROW(rat_coefs) != n_dims) {
stop("Input's dimension mismatch with specified dimensions.")
}
## Get the numerators and denominators of the rational functions
rat_f <- rat_coefs[idx_f, , drop = FALSE]
rat_g <- rbind(rep(1, NCOL(rat_coefs)),
rat_coefs[idx_g, , drop = FALSE])

## Continuity condition
satisfy <- is_no_roots(rat_g, xmin, xmax, EPS)

if (any(satisfy)) {
## Monotone condition
derivs <- get_deriv_nums(rat_f[, satisfy, drop = FALSE],
rat_g[, satisfy, drop = FALSE],
EPS)
satisfy[satisfy] <- is_positives(derivs, xmin, xmax, increasing,
EPS)
}
satisfy
}
}

#' Residual sum of squares of rational function models
#'
#' Produce a function that calculates the residual sum of squares of a set of
#' rational function models.
#'
#' The produced function can be minimised using SMC-SA, multiSA or CEPSO.
#'
#' @param y a response vector. (num vec)
#' @param x a predictor vector. (num vec)
#' @param dim_f the number of coefficients in the numerator. (num)
#' @param dim_g the number of coefficients in the denominator, excluding the
#'     constant term. (num)
#'
#' @return a function that takes in a matrix of rational function coefficients
#'     (numerator precede denominator, in ascending order), each column
#'     corresponds to a rational function, and returns a vector of redisual sum
#'     of squares.
rational_rss_fac <- function(y, x, dim_f, dim_g){

## Evaluate function's arguments eagerly
force(x)
force(y)
force(dim_f)
force(dim_g)

## Type-checking
if (!(is.vector(x) & is.vector(y))) {
stop("Non-vector dataset")
}

## Sequence of indices for both numerator and denominator
idx_f <- seq.int(1, dim_f)
idx_g <- seq.int(dim_f + 1, dim_f + dim_g)
total_dim <- dim_f + dim_g

## Check the dimension of the response
if (is.matrix(y) && (NCOL(y) == 1 || NROW(y) == 1)) {
y <- as.vector(y)
} else if (!is.vector(y)) {
stop("Response data must be a vector / 1-dim matrix.")
}

## A function that returns the residuals for each column of "theta" that
## must be supplied as a matrix
function(theta){
n_row <- NROW(theta)
n_col <- NCOL(theta)

if (is.vector(theta)) {
warning("Argument of residual function supplied as vector.")
theta <- as.matrix(theta)
}

if (n_row != total_dim) {
warning("Argument of residual funciton mismatch with dimension supplied.")
}

fit <- eval_pols(theta[idx_f, , drop=FALSE], x) /
eval_pols(rbind(rep(1, len=n_col), theta[idx_g, , drop=FALSE]), x)
return(colSums((y - fit)^2))
}
}

#' Tukey's biweight of rational function models
#'
#' Produce a function that calculates the Tukey's biweight of a set of rational
#' function models.
#'
#' The explicit formula used is rho(x) = (c^2)/6 * (1-(1-(x/c)^2)^3) if x < |c|,
#' or rho(x) = 1 if x > |c|, where c is an arbitrary constant. The sum of
#' rho(y_i - r(x_i)), where r(.) is the rational function, is then
#' minimised. The produced function can be minimised using SMC-SA, multiSA or
#' CEPSO.
#'
#' @param y a response vector. (num vec)
#' @param x a predictor vector. (num vec)
#' @param dim_f the number of coefficients in the numerator. (num)
#' @param dim_g the number of coefficients in the denominator, excluding the
#'     constant term. (num)
#' @param cst the constant in Tukey's biweight. (num)
#'
#' @return a function that takes in a matrix of rational function coefficients
#'     (numerator precede denominator, in ascending order), each column
#'     corresponds to a rational function, and returns a vector Tukey's
#'     biweight.
rational_biweight_fac <- function(y, x, dim_f, dim_g, cst = 4.685){
force(x)
force(y)
force(dim_f)
force(dim_g)
force(cst)

## Sequence of indices for both numerator and denominator
idx_f <- seq.int(1, dim_f)
idx_g <- seq.int(dim_f + 1, dim_f + dim_g)
total_dim <- dim_f + dim_g

## Check the dimension of the response
if (is.matrix(y) && (NCOL(y) == 1 || NROW(y) == 1)) {
y <- as.vector(y)
} else if (!is.vector(y)) {
stop("Response data must be a vector / 1-dim matrix.")
}

## The biweight function (integration of phi)
## x is a matrix
rho <- function(x, cst = 4.685){
res <- matrix((cst^2)/6, NROW(x), NCOL(x))
outside <- abs(x) > cst
res[!outside] <- (cst^2)/6 * (1-(1-(x[!outside]/cst)^2)^3)
return(res)
}

## A function that return biweight for each column
## "theta" must be supplied as matrix
function(theta) {
n_row <- NROW(theta)
n_col <- NCOL(theta)

if (is.vector(theta)) {
warning("Argument of residual function supplied as vector.")
theta <- as.matrix(theta)
}

if (n_row != total_dim) {
stop("Argument of objective mismatch with the dimension supplied.")
}

fit <- eval_pols(theta[idx_f, , drop=FALSE], x) /
eval_pols(rbind(rep(1, len=n_col), theta[idx_g, , drop=FALSE]), x)
return(colSums(rho(y - fit)))
}

}

#' Line plots of rational functions
#'
#' Plot one or multiple rational functions, given a set of coefficients.

#' @param fs a matrix of numerator polynomial in ascending order representing,
#'     each column corresponds to a polynomial. (num mat)
#' @param gs a matrix of denominator coefficients in ascending order, each
#'     column corresponds to a polynomial. (num mat)
#' @param xlim the x-axis range to plot
#' @param sketch if set to \code{TRUE}, this function add a rational function
#'     onto an existing plot, rather than creating a new graphical device.
#' @param fine how many points between xlims?
#' @param ... optional graphical parameters to be passed to plot function.
plot_rat <- function(fs, gs, xlim = c(0, 10), sketch = TRUE, fine = 1000, ...) {

## fs and gs must be matrices
fs <- as.matrix(fs)
gs <- as.matrix(gs)

x <- seq(xlim, xlim, length.out = fine)
y <- eval_pols(fs, x) / eval_pols(gs, x)

if (sketch) {
lines(y[, 1] ~ x, ...)
} else {
plot(y[, 1] ~ x, type = 'l', ...)
}

## If there are multiple rational functions, plot the rest
if (NCOL(fs) > 1) {
for (i in seq_len(NCOL(fs))[-1]) {
lines(y[, i] ~ x, col = i+1, ...)
}
}
}
```
weiyaw/blackbox documentation built on June 7, 2019, 5:12 a.m.