DiBS <- function(targets, receivers, rho_0, longest_edge, verbose = TRUE, obj = 'avg', 
                 sensor = 'fermi', b = 0.25, gap_method = 'none', gap = 0.05) {
  # Main function to run DiBS
  #
  # Args:
  #   targets: Data frame with object, x, y and value.
  #   receivers: Data frame with object, x and y.
  #   rho_0: Range of the day.
  #   longest_edge: Termination condition.
  #   verbose: If TRUE, prints progress. Default is TRUE.
  #   obj: Objective. Possible values: 'avg', 'min'. Default is 'avg'.
  #   sensor: Sensor model. Possible values: 'cookie', 'fermi', 'exp'. Default is 'fermi'.
  #   b: Diffusivity parameter for Fermi model. Default is 0.25.
  #   gap_method: Optimality gap method. Possible values: 'none', 'center', 'corners', 'corners.center'. Default is 'none'.
  #   gap: Accepted optimality gap. Default is 0.05.
  #
  # Returns:
  #   Final set of sectors, where the first row yields the solution.
  terminate <- FALSE; iteration <- 0
  data <- .create.data(targets, receivers)
  sectors <- .create.sectors(data, rho_0, obj, sensor, b, 'initial')
  while(!terminate) {
    iteration <- iteration + 1
    sectors <- sectors[order(-sectors$upper.bound, sectors$longest.edge),]
    gamma <- sectors[1,]
    gamma$lower.bound <- .lower.bound(data, gamma, rho_0, obj, sensor, b, gap_method)
    gamma$optimality <- gamma$lower.bound / gamma$upper.bound
    if (gamma$longest.edge <= longest_edge | 1 - gamma$optimality <= gap) 
      terminate <- TRUE
    else {
      if (verbose) {cat(paste0('Iteration ', iteration, ':\n')); print(gamma); cat('\n')}
      if (gamma$xmax - gamma$xmin <= longest_edge) xs <- c(gamma$xmin, gamma$xmax)
      else xs <- seq(gamma$xmin, gamma$xmax, length.out = 3)
      if (gamma$ymax - gamma$ymin <= longest_edge) xs <- c(gamma$xmin, gamma$xmax)
      else ys <- seq(gamma$ymin, gamma$ymax, length.out = 3)
      sectors.prime <- .create.sectors(data, rho_0, obj, sensor, b, xs = xs, ys = ys,
                                       name = paste0('iter', iteration))
      sectors <- rbind(sectors[-1,], sectors.prime)
    }
  }
  cat(paste0('Finished after ', iteration, ' iterations (evaluated ',
             nrow(sectors) + iteration - 1, ' sectors):\n'))
  print(gamma)
  return (sectors)
}

.create.data <- function(targets, receivers) {
  # Creates the data frame used by DiBS containing target locations,
  # values and distances to receivers
  #
  # Args:
  #   targets: Data frame with object, x, y and value.
  #   receivers: Data frame with object, x and y.
  #
  # Returns:
  #   Created data frame.
  data <- targets[,c('x','y','value')]
  row.names(data) <- targets$object
  data[paste0('dist.', as.factor(receivers$object))] <- 
    sapply(row.names(receivers), function(r)
      sapply(row.names(targets), function(t) 
        sqrt((targets[t,'x']-receivers[r,'x'])^2+(targets[t,'y']-receivers[r,'y'])^2)))                                                 
  return (data)
}

.create.sectors <- function(data, rho_0, obj, sensor, b, name, xs = NULL, ys = NULL) {
  # Creates sectors from the given data and evaluates 
  # each sector's longest edge and upper bound. Uses xs and ys as
  # slices. If xs and ys are NULL then it creates the initial split.
  #
  # Args:
  #   data: Data set created by function .create.data().
  #   rho_0: Range of the day
  #   obj: Objective.
  #   sensor: Sensor Model.
  #   b: Diffusivity parameter for Fermi model.
  #   name: identifier for the set of sectors
  #   xs: Vector of vertical slices. Default is NULL.
  #   ys: Vector of horizontal slices. Default is NULL.
  #
  # Returns:
  #   Data frame with sector coordinates and upper bounds.
  if (is.null(xs)) xs <- seq(min(data$x), max(data$x), length.out=3)
  if (is.null(ys)) ys <- seq(min(data$y), max(data$y), length.out=3)
  sectors <- data.frame(
    sector = paste0(name, '_', 1:((length(xs) - 1) * (length(ys) - 1))),
    xmin   = rep(xs[-length(xs)], length(ys) - 1), 
    xmax   = rep(xs[-1],          length(ys) - 1),
    ymin   = rep(ys[-length(ys)], each = length(xs) - 1),
    ymax   = rep(ys[-1],          each = length(xs) - 1))
  sectors$longest.edge <- pmax(sectors$xmax - sectors$xmin, sectors$ymax - sectors$ymin)
  sectors$upper.bound <- sapply(1:nrow(sectors), function(s) {
    switch(obj, avg = value <- 0, min = value <- Inf,
           stop(paste0('objective "', obj, '" not known')))
    if (distr == 'cookie') value <- 0
    for(t in 1:nrow(data)) {
      x <- min(max(sectors$xmin[s], data$x[t]), sectors$xmax[s])
      y <- min(max(sectors$ymin[s], data$y[t]), sectors$ymax[s])
      Pt <- .Pt(data, t, x, y, rho_0, sensor, b)
      if(sensor == 'cookie') {
        value <- value + Pt*data[t, 'value']
      } else {
        switch(obj, 
               avg = value <- value + data[t, 'value'] * Pt / nrow(data),
               min = value <- min(value, data[t, 'value'] * Pt)
        )
      }
    }
    return (value)
  })
  return (sectors)
}

.Pt <- function(data, t, x, y, rho_0, sensor, b) {
  # Calculates detection probability Pt for target t with source at position (x, y)
  #
  # Args:
  #   data: Data set created by function .create.data().
  #   t: target
  #   x: source's x coordinate
  #   y: source's y coordinate
  #   rho_0: Range of the day
  #   sensor: Sensor Model.
  #   b: Diffusivity parameter for Fermi model.
  #
  # Returns:
  #   Detection Probability Pt
  d.ts <- sqrt((x - data$x[t])^2 + (y - data$y[t])^2)
  Pt <- 1
  for(r in 4:ncol(data)) {
    rho.tsr <- sqrt(data[t,r] * d.ts)
    switch(sensor,
           cookie = if (rho.tsr <= rho.0) {Pt <- 0; break},
           fermi  = Pt <- Pt * (1 - 1 / (1 + 10^(((rho.tsr / rho.0) - 1) / b))),
           exp    = Pt <- Pt * (1 - 10^(-0.30103 * rho.tsr / rho.0)),
           stop(paste0('distribution "', distr, '" not known'))
    )
  }
  return (1-Pt)
}

.lower.bound <- function(data, gamma, rho_0, obj, sensor, b, gap_method) {
  # calculates a lower bound for sector gamma
  #
  # Args:
  #   data: Data set created by function .create.data().
  #   gamma: Sector data created in function DiBS()
  #   rho_0: Range of the day.
  #   longest_edge: Termination condition.
  #   obj: Objective. Possible values: 'avg', 'min'.
  #   sensor: Sensor model. Possible values: 'cookie', 'fermi', 'exp'.
  #   b: Diffusivity parameter for Fermi model.
  #   gap_method: Optimality gap method. Possible values: 'none', 'center', 'corners', 'corners.center'.
  #
  # Returns:
  #   Lower bound for sector gamma.
  switch(gap_method,
         none           = return (0),
         center         = {
           potentials <- data.frame(x=mean(c(gamma$xmin, gamma$xmax)),
                                    y=mean(c(gamma$ymin, gamma$ymax)))
         },
         corners        = {
           potentials <- data.frame(x=c(gamma$xmin, gamma$xmax, gamma$xmin, gamma$xmax),
                                    y=c(gamma$ymin, gamma$ymin, gamma$ymax, gamma$ymax))
         },
         center.corners = {
           potentials <- data.frame(x=c(gamma$xmin, gamma$xmax, gamma$xmin, gamma$xmax,
                                        mean(c(gamma$xmin, gamma$xmax))),
                                    y=c(gamma$ymin, gamma$ymin, gamma$ymax, gamma$ymax,
                                        mean(c(gamma$ymin, gamma$ymax))))
         },
         stop(paste0('gap method "', gap_method, '" not known'))
  ) 
  potentials$value <- sapply(1:nrow(potentials), function(p) {
    switch(obj, avg = value <- 0, min = value <- Inf,
           stop(paste0('objective "', obj, '" not known')))
    if (distr == 'cookie') value <- 0
    for(t in 1:nrow(data)) {
      Pt <- .Pt(data, t, potentials$x[p], potentials$y[p], rho_0, sensor, b)
      if(sensor == 'cookie') {
        value <- value + Pt*data[t, 'value']
      } else {
        switch(obj, 
               avg = value <- value + data[t, 'value'] * Pt / nrow(data),
               min = value <- min(value, data[t, 'value'] * Pt)
        )
      }
    }
    return (value)
  })
  potentials <- potentials[order(-potentials$value),]
  return (potentials$value[1])
}

# Example
# The next examples demonstrates the use of DiBS.
targets <- data.frame(
  object = c(  'T1',  'T2',  'T3',  'T4',  'T5',  'T6',  'T7',  'T8',  'T9', 'T10'),
  x      = c(-15.00,  6.25,- 9.25,  0.75,- 4.00,  9.00,- 9.50,  2.75,  3.50,- 4.00),
  y      = c(- 4.75,- 9.50,-11.00,  8.75,  5.50,-11.25,-15.00, 11.50,-12.25,  7.50),
  value  = 1)
receivers <- data.frame(
  object = c(  'R1',  'R2',  'R3'),
  x      = c(  4.25, 14.25,- 6.75),
  y      = c(-14.25,  6.25, 12.50))

# 1. fermi model, maximize average detection probability, longest edge <= 0.25
sectors <- DiBS(targets, receivers, rho_0 = 3, longest_edge = 0.25)

# 2. cookie cutter model, longest edge <= 0.1
sectors <- DiBS(targets, receivers, rho_0 = 3, longest_edge = 0.1, sensor='cookie')

# 3. like 1. but longest edge <= 0.001, use center for lower bound
sectors <- sectors <- DiBS(targets, receivers, rho_0 = 3, longest_edge = 0.001,
                           gap_method = 'center', gap = 0.01)