## ----setup, include = FALSE---------------------------------------------------
Sys.setenv(OMP_NUM_THREADS = "2")

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 5
)

## ----load-package-------------------------------------------------------------
library(OptOR)

## ----discrete-array-----------------------------------------------------------
X <- matrix(
  c(
    1, 1, 0, 0,
    1, 1, 0, 0,
    0, 0, 1, 0
  ),
  nrow = 3,
  byrow = TRUE
)

X

## ----discrete-optimization----------------------------------------------------
grid_result <- optimal_grid_hr(
  X,
  verbose = FALSE
)

grid_result

## ----discrete-volume----------------------------------------------------------
grid_widths <- grid_result$ul - grid_result$ll + 1L
grid_widths
prod(grid_widths)

## ----define-fct---------------------------------------------------------------
disk_fct <- list(
  c = 0,
  b = c(0, 0),
  Q = diag(2),
  lim_ul = 1
)

fcts <- list(disk_fct)

## ----classify-grid------------------------------------------------------------
grid_point <- calc_X(
  fcts = fcts,
  n = 8,
  rg_ll = c(-1, -1),
  rg_ul = c(1, 1),
  gmode = "point"
)

grid_point$X

## ----inspect-grid-result------------------------------------------------------
str(grid_point, max.level = 1)

## ----conservative-grid--------------------------------------------------------
grid_conservative <- calc_X(
  fcts = fcts,
  n = 8,
  rg_ll = c(-1, -1),
  rg_ul = c( 1,  1),
  gmode = "conservative"
)

grid_conservative$X

## ----continuous-optimization--------------------------------------------------
cont_result <- optimal_cont_hr(
  fcts = fcts,
  n = 30,
  rg_ll = c(-1, -1),
  rg_ul = c( 1,  1),
  gmode = "conservative",
  verbose = FALSE
)

cont_result

## ----continuous-volume--------------------------------------------------------
cont_widths <- cont_result$ul - cont_result$ll
V_cons <- prod(cont_widths)

cont_widths
V_cons

## ----verify-region------------------------------------------------------------
verification <- find_extrema(
  fct = disk_fct,
  hr_ll = cont_result$ll,
  hr_ul = cont_result$ul
)

verification

## ----verify-upper-limit-------------------------------------------------------
verification$max <= disk_fct$lim_ul + 1e-10

## ----optimistic-bound---------------------------------------------------------
optimistic_result <- optimal_cont_hr(
  fcts = fcts,
  n = 30,
  rg_ll = c(-1, -1),
  rg_ul = c( 1,  1),
  gmode = "optimistic",
  verbose = FALSE
)

optimistic_widths <- optimistic_result$ul - optimistic_result$ll
V_optimistic <- prod(optimistic_widths)

optimistic_widths
V_optimistic

## ----approximation-metric-----------------------------------------------------
d <- 2L

R <- (V_cons / V_optimistic)^(1 / d)
R

## ----prepare-circle-figure, include=FALSE-------------------------------------
n_plot <- 8L
plot_rg_ll <- c(-1, -1)
plot_rg_ul <- c( 1,  1)

plot_fct <- list(
  c = 0,
  b = c(0, 0),
  Q = diag(2),
  lim_ul = 1
)
plot_fcts <- list(plot_fct)

plot_grid_conservative <- calc_X(
  fcts = plot_fcts,
  n = n_plot,
  rg_ll = plot_rg_ll,
  rg_ul = plot_rg_ul,
  gmode = "conservative"
)

plot_grid_optimistic <- calc_X(
  fcts = plot_fcts,
  n = n_plot,
  rg_ll = plot_rg_ll,
  rg_ul = plot_rg_ul,
  gmode = "optimistic"
)

plot_hr_conservative <- optimal_cont_hr(
  fcts = plot_fcts,
  n = n_plot,
  rg_ll = plot_rg_ll,
  rg_ul = plot_rg_ul,
  gmode = "conservative",
  verbose = FALSE
)

# The conservative classification may contain boundary cells coded as 2.
# Only cells coded as 1 are guaranteed feasible and may enter the discrete
# conservative optimization.
plot_X_conservative_binary <- array(
  as.integer(plot_grid_conservative$X == 1L),
  dim = dim(plot_grid_conservative$X)
)

stopifnot(all(plot_X_conservative_binary %in% c(0L, 1L)))

plot_grid_hr_conservative <- optimal_grid_hr(
  plot_X_conservative_binary,
  verbose = FALSE
)

plot_hr_optimistic <- optimal_cont_hr(
  fcts = plot_fcts,
  n = n_plot,
  rg_ll = plot_rg_ll,
  rg_ul = plot_rg_ul,
  gmode = "optimistic",
  verbose = FALSE
)

true_half_width <- 1 / sqrt(2)
plot_hr_true <- list(
  ll = rep(-true_half_width, 2),
  ul = rep(true_half_width, 2)
)

plot_cell_width <- (plot_rg_ul - plot_rg_ll) / n_plot
plot_x_centres <- plot_rg_ll[1] +
  (seq_len(n_plot) - 0.5) * plot_cell_width[1]
plot_y_centres <- plot_rg_ll[2] +
  (seq_len(n_plot) - 0.5) * plot_cell_width[2]

plot_cells <- expand.grid(
  i = seq_len(n_plot),
  j = seq_len(n_plot)
)
plot_cells$x <- plot_x_centres[plot_cells$i]
plot_cells$y <- plot_y_centres[plot_cells$j]

plot_cells$conservative_code <- mapply(
  function(i, j) plot_grid_conservative$X[i, j],
  plot_cells$i,
  plot_cells$j
)
plot_cells$optimistic_code <- mapply(
  function(i, j) plot_grid_optimistic$X[i, j],
  plot_cells$i,
  plot_cells$j
)

plot_cells$classification <- "Outside optimistic approximation"
plot_cells$classification[plot_cells$optimistic_code != 0] <-
  "Optimistic approximation"
plot_cells$classification[plot_cells$conservative_code == 1] <-
  "Conservative"
plot_cells$classification <- factor(
  plot_cells$classification,
  levels = c(
    "Outside optimistic approximation",
    "Optimistic approximation",
    "Conservative"
  )
)

plot_theta <- seq(0, 2 * pi, length.out = 1000)
plot_circle <- data.frame(
  x = cos(plot_theta),
  y = sin(plot_theta)
)

plot_rectangles <- data.frame(
  region = factor(
    c("Conservative", "True optimal", "Optimistic"),
    levels = c("Conservative", "True optimal", "Optimistic")
  ),
  xmin = c(
    plot_hr_conservative$ll[1],
    plot_hr_true$ll[1],
    plot_hr_optimistic$ll[1]
  ),
  xmax = c(
    plot_hr_conservative$ul[1],
    plot_hr_true$ul[1],
    plot_hr_optimistic$ul[1]
  ),
  ymin = c(
    plot_hr_conservative$ll[2],
    plot_hr_true$ll[2],
    plot_hr_optimistic$ll[2]
  ),
  ymax = c(
    plot_hr_conservative$ul[2],
    plot_hr_true$ul[2],
    plot_hr_optimistic$ul[2]
  )
)

circle_plot <- ggplot2::ggplot() +
  ggplot2::geom_tile(
    data = plot_cells,
    ggplot2::aes(x = x, y = y, fill = classification),
    width = plot_cell_width[1],
    height = plot_cell_width[2],
    colour = "black",
    linewidth = 0.25
  ) +
  ggplot2::geom_path(
    data = plot_circle,
    ggplot2::aes(x = x, y = y),
    linewidth = 1.2,
    colour = "black"
  ) +
  ggplot2::geom_rect(
    data = plot_rectangles,
    ggplot2::aes(
      xmin = xmin,
      xmax = xmax,
      ymin = ymin,
      ymax = ymax,
      colour = region
    ),
    fill = NA,
    linewidth = 1.5
  ) +
  ggplot2::scale_fill_manual(
    values = c(
      "Outside optimistic approximation" = "white",
      "Optimistic approximation" = "grey75",
      "Conservative" = "grey45"
    ),
    name = "Grid classification"
  ) +
  ggplot2::scale_colour_manual(
    values = c(
      "Conservative" = "blue",
      "True optimal" = "red",
      "Optimistic" = "orange"
    ),
    name = "Operating region",
    guide = ggplot2::guide_legend(
      order = 1,
      override.aes = list(linewidth = 1.5, fill = NA)
    )
  ) +
  ggplot2::guides(
    fill = ggplot2::guide_legend(order = 2)
  ) +
  ggplot2::coord_equal(
    xlim = c(plot_rg_ll[1], plot_rg_ul[1]),
    ylim = c(plot_rg_ll[2], plot_rg_ul[2]),
    expand = FALSE
  ) +
  ggplot2::labs(
    x = expression(x[1]),
    y = expression(x[2]),
    title = "Grid approximations and optimal rectangles",
    subtitle = expression(x[1]^2 + x[2]^2 <= 1)
  ) +
  ggplot2::theme_bw(base_size = 13) +
  ggplot2::theme(
    panel.border = ggplot2::element_blank(),
    legend.position = "right",
    plot.title.position = "plot",
    legend.key.width = grid::unit(1, "cm"),
    legend.key.height = grid::unit(1, "cm")
  )

V_cons_plot <- prod(plot_hr_conservative$ul - plot_hr_conservative$ll)
V_true_plot <- prod(plot_hr_true$ul - plot_hr_true$ll)
V_optimistic_plot <- prod(plot_hr_optimistic$ul - plot_hr_optimistic$ll)

## ----circle-grid-figure, echo=FALSE, fig.width=10, fig.height=7.5, fig.align='center', out.width='100%', fig.cap='Conservative and optimistic grid approximations of the unit disk. The conservative operating region is guaranteed feasible, the optimistic region provides a global upper volume bound, and the analytically known maximum volume lies between them.'----
circle_plot

## ----circle-approximation-metrics---------------------------------------------
d_plot <- 2L

R_true <- (V_cons_plot / V_true_plot)^(1 / d_plot)
R_outer <- (V_cons_plot / V_optimistic_plot)^(1 / d_plot)

c(
  conservative_to_true = R_true,
  conservative_to_outer_bound = R_outer
)

## ----width-constraint---------------------------------------------------------
width_result <- optimal_cont_hr(
  fcts = fcts,
  n = 30,
  rg_ll = c(-1, -1),
  rg_ul = c( 1,  1),
  ctype = c("width", "none"),
  cwidth = c(1.8, 0),
  gmode = "conservative",
  verbose = FALSE
)

width_result$ul - width_result$ll