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))