cardoid_membership <- function(x) {
  a <- 1/2
  r2 <- apply(x, 1, function(x) x[[1]]^2+(x[[2]]-1)^2)
  theta <- apply(x, 1, function(x) atan2(x[[2]]-1, x[[1]]))
  card_membership <- r2 <= a*(1-sin(theta))
  c_rad <- 0.1
  r2 <- apply(x, 1, function(x) (x[[1]]+0.1)^2+(x[[2]]-0.6)^2)
  circ_membership <- r2 <= c_rad
  return(card_membership & !circ_membership)
}

t_df <- expand.grid(x = seq(-1, 1, length.out = 100), y = seq(0, 1.5, length.out = 100))
mem <- cardoid_membership(x = t_df)
plot(t_df[mem,], pch = 16, cex = 0.5, xlim = c(-1, 1), ylim = c(0, 1.5))

# Scales inputs: important since the emulators should take inputs purely in [-1,1]
scale_input <- function(x, r, forward = TRUE) {
  centers <- purrr::map(r, ~(.x[2]+.x[1])/2)
  scales <- purrr::map(r, ~(.x[2]-.x[1])/2)
  if (is.null(names(x))) {
    centers <- unlist(centers, use.names = F)
    scales <- unlist(scales, use.names = F)
  }
  if (forward)
    output <- (x-centers)/scales
  else
    output <- x * scales + centers
  if (!"data.frame" %in% class(output))
    return(data.frame(output))
  return(output)
}
# Evaluate multiple functions over points
eval_funcs <- function(funcs, points, ...) {
  pointsdim <- (length(dim(points)) != 0)
  manyfuncs <- (typeof(funcs) != "closure")
  if (manyfuncs && pointsdim)
    return(apply(points, 1, function(x) purrr::map_dbl(funcs, purrr::exec, x, ...)))
  if (manyfuncs)
    return(purrr::map_dbl(funcs, purrr::exec, points, ...))
  if (pointsdim) {
    return(tryCatch(apply(points, 1, funcs, ...),
                    error = function(cond1) {
                      tryCatch(purrr::exec(funcs, points, ...),
                               error = function(cond2) {
                                 print(cond1, cond2)
                               })
                    }))
  }
  return(purrr::exec(funcs, points, ...))
}
## LHS Generation
lhs_gen_card <- function(ranges, n_points, points.factor = 10) {
  points <- eval_funcs(scale_input, setNames(data.frame(2 * (lhs::randomLHS(n_points * points.factor, length(ranges)) - 0.5)), names(ranges)), ranges, FALSE)
  point_imps <- cardoid_membership(points)
  return(list(points = points, imps = point_imps))
}
line_samp_card <- function(s_points, ranges, n_lines = 20, ppl = 50, ...) {
  in_range <- function(data, ranges) {
    apply(data, 1, function(x) all(purrr::map_lgl(seq_along(ranges), ~x[.] >= ranges[[.]][1] && x[.] <= ranges[[.]][2])))
  }
  if (nrow(s_points) < 2) return(s_points)
  n_lines <- min(nrow(s_points)*(nrow(s_points)-1)/2, n_lines)
  if (ppl %% 4 == 1) ppl <- ppl + 1
  s_lines <- lapply(1:(10*n_lines), function(x) {
    pts <- s_points[sample(nrow(s_points), 2),]
    pt_dist <- dist(pts)
    return(list(p1 = pts[1,], p2 = pts[2,], d = pt_dist))
  })
  s_lines <- s_lines[!duplicated(purrr::map_dbl(s_lines, ~.$d))]
  best_pts <- s_lines[order(purrr::map_dbl(s_lines, ~.$d), decreasing = TRUE)][1:n_lines]
  get_limits <- function(points) {
    point_dist <- sqrt(sum((points[[1]]-points[[2]])^2))
    range_dist <- sqrt(sum(purrr::map_dbl(ranges, ~(.[[2]]-.[[1]])^2)))
    dist_ratio <- range_dist/point_dist
    l_seq <- seq(1-dist_ratio, dist_ratio, length.out = 500)
    line_points <- setNames(do.call('rbind.data.frame', purrr::map(l_seq, ~points[[1]] + .*(points[[2]]-points[[1]]))), names(ranges))
    is_in_range <- in_range(line_points, ranges)
    valid_ls <- l_seq[is_in_range]
    return(c(valid_ls[1], valid_ls[length(valid_ls)]))
  }
  samp_pts <- lapply(best_pts, function(x) {
    #these_lims <- get_limits(x)
    tryCatch(
      {
        do.call('rbind', lapply(seq(-x[[3]], x[[3]], length.out = ppl), function(y) (x[[1]]+x[[2]])/2 + y*(x[[2]]-x[[1]])))
        #do.call('rbind', lapply(seq(these_lims[1], these_lims[2], length.out = ppl), function(y) (x[[1]]+x[[2]])/2 + y*(x[[2]]-x[[1]]))),
      },
      error = function(e) {
        return(NULL)
      }
    )
  })
  samp_pts <- samp_pts[!purrr::map_lgl(samp_pts, is.null)]
  samp_pts <- purrr::map(samp_pts, ~.[in_range(., ranges),])
  samp_pts <- samp_pts[!purrr::map_lgl(samp_pts, ~is.null(nrow(.)) || is.null(.) || nrow(.) == 0)]
  imps <- purrr::map(samp_pts, cardoid_membership)
  include_pts <- purrr::map(seq_along(samp_pts), function(x) {
    pts <- samp_pts[[x]]
    imp <- imps[[x]]
    included <- purrr::map_lgl(seq_along(imp), function(y) {
      if (!imp[y]) return(FALSE)
      if (y == 1 || y == length(imp)) return(TRUE)
      if (!imp[y+1] || !imp[y-1]) return(TRUE)
      return(FALSE)
    })
    return(pts[included,])
  })
  out_df <- do.call('rbind', include_pts)
  uniqueness <- row.names(unique(signif(out_df, 6)))
  return(list(points = do.call('rbind', samp_pts), final_pts = out_df[uniqueness, ]))
}
runifs <- function(n, d, c = rep(0, d), r = 1) {
  init_points <- matrix(rnorm(n*d, 0, 1), nrow = n, byrow = T)
  exps <- rexp(n, 0.5)
  denoms <- (exps + apply(init_points, 1, function(x) sum(x^2)))^(0.5)
  swept <- sweep(init_points, 1, denoms, "/")
  return(t(apply(swept, 1, function(x) x*r + c)))
}
punifs <- function(x, c = rep(0, length(x)), r = 1) {
  return(ifelse(sum((x-c)^2 / r^2) <= 1, 1, 0))
}
imp_card <- function(n_points, s_points, ranges, cutoff = 3, nth = 1, distro = 'sphere', sd = NULL, to_file = NULL, ...) {
  if (nrow(s_points) >= n_points)
    return(s_points)
  m_points <- n_points - nrow(s_points)
  new_points <- s_points
  in_range <- function(data, ranges) {
    apply(data, 1, function(x) all(purrr::map_lgl(seq_along(ranges), ~x[.] >= ranges[[.]][1] && x[.] <= ranges[[.]][2])))
  }
  s_trafo <- sweep(sweep(s_points, 2, apply(s_points, 2, mean), "-"), 2, apply(s_points, 2, sd), "/")
  s_estruct <- eigen(cov(s_trafo))
  s_estruct$values <- purrr::map_dbl(s_estruct$values, function(x) {if(x < 1e-10) 1e-10 else x})
  pca_transform <- function(x, forward = TRUE) {
    if (forward) x <- sweep(sweep(x, 2, apply(s_points, 2, mean), "-"), 2, apply(s_points, 2, sd), "/")
    if ("data.frame" %in% class(x)) x <- data.matrix(x)
    if (forward) return(x %*% s_estruct$vectors %*% diag(1/sqrt(s_estruct$values)))
    pre_traf <- x %*% diag(sqrt(s_estruct$values)) %*% t(s_estruct$vectors)
    return(sweep(sweep(pre_traf, 2, apply(s_points, 2, sd), "*"), 2, apply(s_points, 2, mean), "+"))
  }
  propose_points <- function(sp, sd, how_many = n_points) {
    sp_trafo <- pca_transform(sp)
    wp <- sp_trafo[sample(nrow(sp_trafo), how_many, replace = TRUE), ]
    if (distro == "normal") pp <- t(apply(wp, 1, function(x) mvtnorm::rmvnorm(1, mean = unlist(x, use.names = F), sigma = sd)))
    else pp <- t(apply(wp, 1, function(x) runifs(1, length(s_points), unlist(x, use.names = F), r = sd)))
    prop_points <- data.frame(pp)
    back_traf <- setNames(data.frame(pca_transform(pp, FALSE)), names(ranges))
    valid <- in_range(back_traf, ranges) & cardoid_membership(back_traf)
    if (distro == "normal") {
      tweights <- apply(prop_points, 1, function(x) 1/nrow(sp_trafo) * sum(apply(sp_trafo, 1, function(y) mvtnorm::dmvnorm(x, mean = y, sigma = sd))))
      min_w <- min(apply(sp_trafo, 1, function(x) 1/nrow(sp_trafo) * sum(apply(sp_trafo, 1, function(y) mvtnorm::dmvnorm(x, mean = y, sigma = sd)))))
      weights <- min_w/tweights
    }
    else
      weights <- apply(prop_points, 1, function(x) sum(apply(sp_trafo, 1, function(y) punifs(x, y, r = sd))))
    allow <- runif(length(weights)) < 1/weights
    accepted <- valid & allow
    return(back_traf[accepted,])
  }
  if (is.null(sd)) {
    if (distro == "normal")
      sd <- diag(2, length(ranges))
    else
      sd <- rep(2, length(ranges))
  }
  accept_rate <- NULL
  upper_accept <- 0.125
  lower_accept <- 0.075
  while ((is.null(accept_rate) || accept_rate > upper_accept || accept_rate < lower_accept) && nrow(new_points) < n_points) {
    if (!is.null(accept_rate)) {
      if (accept_rate > upper_accept) sd <- sd * 1.1
      else sd <- sd * 0.9
    }
    how_many <- max(floor(n_points/4), 250)
    prop_points <- propose_points(s_points, sd, how_many)
    new_points <- rbind(new_points, prop_points)
    uniqueness <- row.names(unique(signif(new_points, 7)))
    new_points <- new_points[uniqueness,]
    if (!is.null(to_file)) write.csv(new_points, file = to_file, row.names = FALSE)
    accept_rate <- nrow(prop_points)/how_many
  }
  while (nrow(new_points) < n_points) {
    prop_points <- propose_points(s_points, sd, ceiling(1.5*m_points/accept_rate))
    new_points <- rbind(new_points, prop_points)
    uniqueness <- row.names(unique(signif(new_points, 7)))
    new_points <- new_points[uniqueness,]
    if (!is.null(to_file)) write.csv(new_points, file = to_file, row.names = FALSE)
  }
  return(new_points)
}

t_rang <- list(x = c(-1,1), y = c(0, 1.5))
t_lhs <- lhs_gen_card(list(x = c(-1, 1), y = c(0, 1.5)), 75)
plot(t_lhs$points, pch = 16, cex = ifelse(t_lhs$imps, 1, 0.5), col = ifelse(t_lhs$imps, 'black', 'grey50'), xlim = c(-1, 1), ylim = c(0, 1.5))
lhs_valid <- t_lhs$points[t_lhs$imps,]
t_line <- line_samp_card(lhs_valid, t_rang, n_lines = 50, ppl = 50)
plot(t_line$points, pch = 16, cex = 0.5, col = 'grey50', xlim = c(-1, 1), ylim = c(0, 1.5))
points(t_line$final_pts, pch = 16, cex = 1, col = 'black')
lhs_samp <- sample(1:nrow(lhs_valid), 300)
current_set <- rbind(lhs_valid[lhs_samp,], t_line$final_pts)
t_imp <- imp_card(2000, current_set, t_rang)
plot(current_set, pch = 16, cex = 0.5, col = 'grey50', ylim = c(0, 1.5), xlim = c(-1, 1))
points(t_imp, pch = 16, cex = 0.8, col = 'black', xlim = c(-1, 1), ylim = c(0, 1.5))


c_measure <- op <- NULL
for (i in 1:1000) {
  tp <- t_imp[sample(nrow(t_imp), min(nrow(t_imp), 1000)),]
  if (!"data.frame" %in% class(tp)) tp <- setNames(data.frame(tp), names(t_rang))
  measure <- min(dist(tp))
  if (is.null(c_measure) || measure > c_measure) {
    op <- tp
    c_measure <- measure
  }
}
points <- op

plot(t_imp, pch = 16, cex = 0.5, col = 'grey50', xlim = c(-1, 1), ylim = c(0, 1.5))
points(points, pch = 16, cex = 0.75, col = 'black', xlim = c(-1, 1), ylim = c(0, 1.5))

plot(points, pch = 16, cex = 0.75, col = 'black', xlim = c(-1, 1), ylim = c(0, 1.5))