# Josua Kugler, Christian Merten

# a) ----------------------------------------------------------------------- 
midpoint <- function(f, a, b) (b-a) * f((a+b)/2)
trapezoid <- function(f, a, b) (b-a) * (f(a) + f(b))/2

# b) ----------------------------------------------------------------------
nc_integrate <- function(f, lower, upper, n, rule) {
    xs <- seq(lower, upper, length.out=n+1)
    sum(rule(f, xs[-(n+1)], xs[-1]))
}

# c) ------------------------------------------------------------------------
closeable_seq <- function(from, to, len, closed=TRUE) {
    if (closed) seq(from, to, length.out = len)
    else seq(from, to, length.out = len+2)[-c(1, len+2)]
}

newton_cotes <- function(coef, closed=TRUE) {
    m <- length(coef)
    w <- coef / sum(coef)
    function(f, a, b) {
        t <- mapply(closeable_seq, a, b, len=m, closed=closed)
        (b-a) * colSums(w*matrix(f(t), ncol=length(a)))
    }
}

# d) ------------------------------------------------------------------------

# some objects for param_fun-list
sin1x <- function(x) {
  y <- suppressWarnings(sin(1/x))
  y[is.na(y)] <- 0
  y
}
set.seed(0)
x <- c(0:1, runif(5))
y <- runif(7)

# list of functions with integral interval
param_fun <- list(
  poly = list(f = function(x) x^4 - x^3 - 3*x^2 + x + 2, lower = 0, upper = 2),
  sin1x = list(f = sin1x, lower = 0, upper = 1),
  lin = list(f = approxfun(x, y), lower = 0, upper = 1),
  spline = list(f = splinefun(x, y), lower = 0, upper = 1)
)

# true values of integrals
true <- c(
  poly = 0.4, 
  sin1x = 0.504067061906928, 
  lin = 0.472878602164825, 
  spline = 0.97236924451286)

# options for creating integral-rules
param_rule <- list(
  midpoint = list(coef = 1, closed=FALSE),
  trapezoid = list(coef = c(1,1), closed=TRUE),
  simpson = list(coef = c(1,4,1), closed=TRUE),
  boole = list(coef = c(7,32,12,32,7), closed=TRUE),
  open5 = list(coef = c(611, -453, 562, 562, -453, 611), closed=FALSE)
)

# how many function evaluations are allowed
param_n <- 5*(2:20)

library(tidyverse)

param <- as_tibble(
  expand.grid(n = param_n,
              fun_name = names(param_fun),
              rule_name = names(param_rule),
              stringsAsFactors=FALSE))

nc_rules <- sapply(param_rule, do.call, what=newton_cotes)

# run nc_integrate with given options
call_nc <- function(n, fun_name, rule_name) {
  opts <- param_fun[[fun_name]]
  opts$rule = nc_rules[[rule_name]]
  # to make things fair:
  # a rule which evaluates f at m points may only be called n/m times
  opts$n = round(n / length(param_rule[[rule_name]]$coef))
  do.call(nc_integrate, opts)
}

param %>% rowwise() %>%
    mutate(value=nc_integrate(param_fun[[fun_name]]$f,
                              lower = param_fun[[fun_name]]$lower,
                              upper = param_fun[[fun_name]]$upper,
                              n = n,
                              rule = nc_rules[[rule_name]]),
           true=true[[fun_name]],
           error=abs(true-value)) -> res

# plot results
plots <- lapply(names(param_fun), function(nm) 
  res %>%
    filter(fun_name == nm) %>%
    ggplot(aes(x = n, y = error, color = rule_name)) +
    scale_y_log10() +
    geom_line() + geom_point() + labs(title = nm)
)
gridExtra::grid.arrange(grobs = plots, nrow=2)