# R/fourierin.R In gbasulto/fourierin: Computes Numeric Fourier Integrals

#### Documented in fourierinfourierin_1dfourierin_2d

## -------------------------------------------------------------------
#' Univariate Fourier integrals
#'
#' It computes Fourier integrals of functions of one and two
#' variables on a regular grid.
#'
#' See vignette for more detailed examples.
#'
#' @param f function or a vector of size m. If a function is
#'     provided, it must be able to be evaluated at vectors. If a
#'     vector of values is provided, such evaluations must have been
#'     obtained on a regular grid and the Fourier integral is faster
#'     is m is a power of 2.
#' @param lower_int Lower integration limit(s).
#' @param upper_int Upper integration limit(s).
#' @param lower_eval Lower evaluation limit(s). It can be NULL if an
#'     evaluation grid is provided.
#' @param upper_eval Upper evaluation limit(s). It can be NULL if an
#'     evaluation grid is provided.
#'     transform. It is usually equal to 0, -1 or 1.
#'     of the Fourier transform. It is usually equal to 1, -1, 2pi or
#'     -2pi.
#' @param resolution A vector of integers (faster if powers of two)
#'     determining the resolution of the evaluation grid. Not
#'     required if f is a vector.
#' @param eval_grid Optional matrix with d columns with the points
#'     where the Fourier integral will be evaluated. If it is
#'     provided, the FFT will not be used.
#' @param use_fft Logical value specifying whether the FFT will be
#'     used.
#' @return If w is given, only the values of the Fourier integral are
#'     returned, otherwise, a list with the elements \item{w}{A vector
#'     of size m where the integral was computed.}  \item{values}{A
#'     complex vector of size m with the values of the integral}
#'
#' @examples
#' ##--- Example 1 ---------------------------------------------------
#' ##--- Recovering std. normal from its characteristic function -----
#' library(fourierin)
#'
#' #' Function to to be used in integrand
#' myfun <- function(t) exp(-t^2/2)
#'
#'                                         # Compute Foueien integral
#' out <- fourierin_1d(f = myfun,
#'                     lower_int = -5, upper_int = 5,
#'                     lower_eval = -3, upper_eval = 3,
#'                     resolution = 64)
#'
#' ## Extract grid and values
#' grid <- out$w #' values <- Re(out$values)
#'
#' plot(grid, values, type = "l", col = 3)
#' lines(grid, dnorm(grid), col = 4)
#'
#' ##--- Example 2 -----------------------------------------------
#' ##--- Computing characteristic function of a gamma r. v. ------
#'
#' library(fourierin)
#'
#' ## Function to to be used in integrand
#' myfun <- function(t) dgamma(t, shape, rate)
#'
#' ## Compute integral
#' shape <- 5
#' rate <- 3
#' out <- fourierin_1d(f = myfun, lower_int = 0, upper_int = 6,
#'                     lower_eval = -4, upper_eval = 4,
#'
#' grid <- out$w # Extract grid #' re_values <- Re(out$values)             # Real values
#' im_values <- Im(out$values) # Imag values #' #' # Now compute the real and #' # imaginary true values of the #' # characteric function. #' true_cf <- function(t, shape, rate) (1 - 1i*t/rate)^-shape #' true_re <- Re(true_cf(grid, shape, rate)) #' true_im <- Im(true_cf(grid, shape, rate)) #' #' # Compare them. We can see a #' # slight discrepancy on the #' # tails, but that is fixed #' # when resulution is #' # increased. #' plot(grid, re_values, type = "l", col = 3) #' lines(grid, true_re, col = 4) #' #' # Same here #' plot(grid, im_values, type = "l", col = 3) #' lines(grid, true_im, col = 4) #' @export fourierin_1d <- function(f, lower_int, upper_int, lower_eval = NULL, upper_eval = NULL, const_adj, freq_adj, resolution = NULL, eval_grid = NULL, use_fft = TRUE) { ## Condensed notation a <- lower_int b <- upper_int c <- lower_eval d <- upper_eval r <- const_adj s <- freq_adj resol <- resolution w <- eval_grid ## fourierin_1d <- function(f, a, b, c = NULL, d = NULL, ## r, s, resol = NULL, w = NULL, ## use_fft = TRUE) { ## Flag to determine wheter a list or a vector will be returned w_given <- !is.null(w) ## If function values are provided, then the resolution ## is the length of the vector of values. if (!is.function(f)) resol <- length(f) ## Increment in the frequency domain. gam <- (d - c)/resol ## Freq. dom. vector. If w is provided, FFT will NOT be used. if (is.null(w)) { if (is.null(c) | is.null(d)) { stop("c and d must be provided.") } w <- seq(c, d - gam, length.out = resol) } else { use_fft <- FALSE } ## If f is the function, it needs to be evaluated in ## the time domain values. if (is.function(f)) { del <- (b - a)/resol # Increment in the time domain. t <- seq(a + del/2, b - del/2, length.out = resol) # Freq. dom. vector. f_t <- f(t) # Function values ## Rutinary check if(is.null(f_t)) stop("Function f is null.") } else { f_t <- f } if (!use_fft) { out <- switch(is.complex(f_t) + 1, fourierin_1d_nonregular_cpp(f_t, a, b, w, resol, r, s), fourierin_cx_1d_nonregular_cpp(f_t, a, b, w, resol, r, s)) } else { out <- switch(is.complex(f_t) + 1, fourierin_1d_cpp(f_t, a, b, c, d, r, s), fourierin_cx_1d_cpp(f_t, a, b, c, d, r, s)) } ## If w is given, return only the values of the integral, ## otherwise alse return w. if(w_given) return(out) return(list(w = w, # Return list. values = drop(out))) } #' Bivariate Fourier integrals #' #' It computes Fourier integrals for functions of one and two #' variables. #' #' @inheritParams fourierin_1d #' @return If w is given, only the values of the Fourier integral are #' returned, otherwise, a list with three elements #' \item{w1}{Evaluation grid for first entry} \item{w2}{Evaluation #' grid for second entry} \item{values}{m1 x m2 matrix of complex #' numbers, corresponding to the evaluations of the integral} #' #' @examples #' ##--- Recovering std. normal from its characteristic function ----- #' library(fourierin) #' #' ##-Parameters of bivariate normal distribution #' mu <- c(-1, 1) #' sig <- matrix(c(3, -1, -1, 2), 2, 2) #' #' ##-Multivariate normal density #' ##-x is n x d #' f <- function(x) { #' ##-Auxiliar values #' d <- ncol(x) #' z <- sweep(x, 2, mu, "-") #' #' ##-Get numerator and denominator of normal density #' num <- exp(-0.5*rowSums(z * (z %*% solve(sig)))) #' denom <- sqrt((2*pi)^d*det(sig)) #' #' return(num/denom) #' } #' #' ##-Characteristic function #' ##-s is n x d #' phi <- function(s) { #' complex(modulus = exp(- 0.5*rowSums(s*(s %*% sig))), #' argument = s %*% mu) #' } #' #' ##-Approximate cf using Fourier integrals #' eval <- fourierin_2d(f, lower_int = c(-8, -6), upper_int = c(6, 8), #' lower_eval = c(-4, -4), upper_eval = c(4, 4), #' const_adj = 1, freq_adj = 1, #' resolution = c(128, 128)) #' #' ## Extract values #' t1 <- eval$w1
#' t2 <- eval$w2 #' t <- as.matrix(expand.grid(t1 = t1, t2 = t2)) #' approx <- eval$values
#' true <- matrix(phi(t), 128, 128)        # Compute true values
#'
#' ##-This is a section of the characteristic functions
#' i <- 65
#' plot(t2, Re(approx[i, ]), type = "l", col = 2,
#'      ylab = "",
#'      xlab = expression(t[2]),
#'      main = expression(paste("Real part section at ",
#'                              t[1], "= 0")))
#' lines(t2, Re(true[i, ]), col = 3)
#' legend("topleft", legend = c("true", "approximation"),
#'        col = 3:2, lwd = 1)
#'
#' ##-Another section, now of the imaginary part
#' plot(t1, Im(approx[, i]), type = "l", col = 2,
#'      ylab = "",
#'      xlab = expression(t[1]),
#'      main = expression(paste("Imaginary part section at ",
#'                              t[2], "= 0")))
#' lines(t1, Im(true[, i]), col = 3)
#' legend("topleft", legend = c("true", "approximation"),
#'        col = 3:2, lwd = 1)
#' @export
fourierin_2d <- function(f, lower_int, upper_int,
lower_eval = NULL, upper_eval = NULL,
resolution = NULL, eval_grid = NULL,
use_fft = TRUE){
## Condensed notation
a <- lower_int
b <- upper_int
c <- lower_eval
d <- upper_eval
resol <- resolution
w <- eval_grid

## If function values are provided, then the resolution is the
## length of the vector of values.
if(!is.function(f)) resol <- dim(f)

## Increment in the frequency domain.
gam <- (d - c)/resol

## Freq. dom. vectors.
if (is.null(w)){
if (is.null(c) | is.null(d)) {
stop("c and d must be provided.")
}
w1 <- seq(c[1], d[1] - gam[1], length.out = resol[1])
w2 <- seq(c[2], d[2] - gam[2], length.out = resol[2])
}

## If f is the function, it needs to be evaluated in the time
## domain values.
if(is.function(f)){
del <- (b - a)/resol # Increment in the time domain.
t1 <- seq(a[1] + del[1]/2, b[1] - del[1]/2,
length.out = resol[1]) # Freq. dom. vector.
t2 <- seq(a[2] + del[2]/2, b[2] - del[2]/2,
length.out = resol[1]) # Freq. dom. vector.
t <- as.matrix(expand.grid(t1, t2))
f_vals <- matrix(f(t), resol[1], resol[2])

# Rutinary check
if(is.null(f_vals)) stop("Function f is null.")

} else{
f_vals <- f
}

if (!use_fft) {
w_temp <- switch(is.null(w) + 1,
w,
as.matrix(expand.grid(w1, w2))
)

out <- switch(is.complex(f_vals) + 1,
fourierin_2d_nonregular_cpp(f_vals, a, b,
w_temp,
resol, r, s),
fourierin_cx_2d_nonregular_cpp(f_vals, a, b,
w_temp, resol,
r, s))
## If no grid was provided, the values of the Fourier integral
## are put into a matrix.
if(is.null(w)) out <- matrix(out, resol[1], resol[2])
} else {
out <- switch(is.complex(f_vals) + 1,
fourierin_2d_cpp(f_vals, a, b, c, d, r, s),
fourierin_cx_2d_cpp(f_vals, a, b, c, d, r, s))
}

## Return grid if it was not provided.
if (is.null(w)) {
return(list(w1 = w1,
w2 = w2,
values = out))

} else {
return(out)
}
}

#' Compute Fourier integrals
#'
#' It computes Fourier integrals for functions of one and two
#' variables.
#'
#' See plenty of detailed examples in the vignette.
#'
#' @inheritParams fourierin_1d
#' @return A list with the elements n-dimensional array and n vectors
#'     with their corresponding resolution. Specifically,
#'     \item{values}{A n-dimensional (resol_1 x resol_2 x ... x
#'     resol_n) complex array with the values.}  \item{w1}{A vector of
#'     size resol_1} \item{...}{ } \item{wn}{A vector of size resol_n}
#'
#' @examples
#' ##--- Example 1 ---------------------------------------------------
#' ##--- Recovering std. normal from its characteristic function -----
#' library(fourierin)
#'
#' ## Function to be used in the integrand
#' myfnc <- function(t) exp(-t^2/2)
#'
#' ## Compute integral
#' out <- fourierin(f = myfnc, lower_int = -5, upper_int = 5,
#'                  lower_eval= -3, upper_eval = 3, const_adj = -1,
#'                  freq_adj = -1, resolution = 64)
#'
#' ## Extract grid and values
#' grid <- out$w #' values <- Re(out$values)
#'
#' ## Compare with true values of Fourier transform
#' plot(grid, values, type = "l", col = 3)
#' lines(grid, dnorm(grid), col = 4)
#'
#'
#' ##--- Example 2 ---------------------------------------------------
#' ##--- Computing characteristic function of a gamma r. v. ----------
#'
#' library(fourierin)
#'
#' ## Function to be used in integrand
#' myfnc <- function(t) dgamma(t, shape, rate)
#'
#' ## Compute integral
#' shape <- 5
#' rate <- 3
#' out <- fourierin(f = myfnc, lower_int = 0, upper_int = 6,
#'                  lower_eval = -4, upper_eval = 4,
#'
#' ## Extract values
#' grid <- out$w # Extract grid #' re_values <- Re(out$values)             # Real values
#' im_values <- Im(out$values) # Imag values #' #' ## Now compute the real and imaginary true values of the #' ## characteric function. #' true_cf <- function(t, shape, rate) (1 - 1i*t/rate)^-shape #' true_re <- Re(true_cf(grid, shape, rate)) #' true_im <- Im(true_cf(grid, shape, rate)) #' #' ## Compare them. We can see a slight discrepancy on the tails, #' ## but that is fixed when resulution is increased. #' plot(grid, re_values, type = "l", col = 3) #' lines(grid, true_re, col = 4) #' #' # Same here #' plot(grid, im_values, type = "l", col = 3) #' lines(grid, true_im, col = 4) #' #' ##--- Example 3 ------------------------------------------------- #' ##--- Recovering std. normal from its characteristic function --- #' library(fourierin) #' #' ##-Parameters of bivariate normal distribution #' mu <- c(-1, 1) #' sig <- matrix(c(3, -1, -1, 2), 2, 2) #' #' ##-Multivariate normal density #' ##-x is n x d #' f <- function(x) { #' ##-Auxiliar values #' d <- ncol(x) #' z <- sweep(x, 2, mu, "-") #' ##-Get numerator and denominator of normal density #' num <- exp(-0.5*rowSums(z * (z %*% solve(sig)))) #' denom <- sqrt((2*pi)^d*det(sig)) #' return(num/denom) #' } #' #' ## Characteristic function #' ## s is n x d #' phi <- function(s) { #' complex(modulus = exp(- 0.5*rowSums(s*(s %*% sig))), #' argument = s %*% mu) #' } #' #' ##-Approximate cf using Fourier integrals #' eval <- fourierin(f, lower_int = c(-8, -6), upper_int = c(6, 8), #' lower_eval = c(-4, -4), upper_eval = c(4, 4), #' const_adj = 1, freq_adj = 1, #' resolution = c(128, 128)) #' #' ## Extract values #' t1 <- eval$w1
#' t2 <- eval$w2 #' t <- as.matrix(expand.grid(t1 = t1, t2 = t2)) #' approx <- eval$values
#' true <- matrix(phi(t), 128, 128)        # Compute true values
#'
#'
#' ## This is a section of the characteristic function
#' i <- 65
#' plot(t2, Re(approx[i, ]), type = "l", col = 2,
#'      ylab = "",
#'      xlab = expression(t[2]),
#'      main = expression(paste("Real part section at ",
#'                              t[1], "= 0")))
#' lines(t2, Re(true[i, ]), col = 3)
#' legend("topleft", legend = c("true", "approximation"),
#'        col = 3:2, lwd = 1)
#'
#' ##-Another section, now of the imaginary part
#' plot(t1, Im(approx[, i]), type = "l", col = 2,
#'      ylab = "",
#'      xlab = expression(t[1]),
#'      main = expression(paste("Imaginary part section at ",
#'                              t[2], "= 0")))
#' lines(t1, Im(true[, i]), col = 3)
#' legend("topleft", legend = c("true", "approximation"),
#'        col = 3:2, lwd = 1)
#' @export
fourierin <- function(f, lower_int, upper_int,
lower_eval = NULL, upper_eval = NULL,
resolution = NULL, eval_grid = NULL,
use_fft = TRUE){

## Get dimension of function from lower integration limit.
n <- length(lower_int)

switch(n,
return(fourierin_1d(f, lower_int, upper_int,
lower_eval, upper_eval,
eval_grid, use_fft)),
return(fourierin_2d(f, lower_int, upper_int,
lower_eval, upper_eval,
resolution, eval_grid, use_fft))
) # End switch
}

#' @useDynLib fourierin, .registration=TRUE
#' @importFrom Rcpp sourceCpp
NULL