R语言等高线图绘制代码优化求助：提升绘制效率

Hey there! Let's fix that slow contour plot code of yours. Drawing thousands of individual points to fill a region is really inefficient—we can do way better by leaning into R's vectorized operations and built-in plotting functions instead. Let's break this down step by step.

First, let's address the core issues in your original code #

1. Slow point plotting: Using points() to render over a million individual points (with your 0.001 step size) is a huge bottleneck. R has to handle each point separately, which drags down performance.
2. Constraint logic mismatch: Your h() function marks regions where either c1 or c2 is positive as non-zero, but you specified you only want to fill the region where both c1 and c2 are positive. Let's correct that first.

Optimized Base R Solution #

We'll use vectorized operations (no more slow sapply!) and mask-based filling to render the grey region in one go, not thousands of points.

# Define your grid (keep the 0.001 step if you need ultra-smooth results)
x1 <- seq(0, 1, 0.001)
x2 <- seq(0, 1, 0.001)

# Vectorized function for f (no need for a separate return statement)
f <- function(x1, x2) x1 + x2

# Vectorized function to create a mask for the region we want to fill
# Returns TRUE where both c1 > 0 AND c2 > 0
constraint_mask <- function(x1, x2) {
  c1 <- 1.5 - x1 - 2*x2 - 0.5*sin(2*pi*(x1^2 - 2*x2))
  c2 <- x1^2 + x2^2 - 1.5
  c1 > 0 & c2 > 0
}

# Generate the f values and constraint mask as matrices
z <- outer(x1, x2, f)
mask <- outer(x1, x2, constraint_mask)

# Plot the feasible region's f values first
image(x1, x2, z, xlab = expression(x[1]), ylab = expression(x[2]))
contour(x1, x2, z, add = TRUE)

# Fill the grey region efficiently using image()
# Create a color matrix: lightgrey where mask is TRUE, transparent otherwise
color_mat <- ifelse(mask, "lightgrey", NA)
image(x1, x2, color_mat, add = TRUE, col = c(NA, "lightgrey"))

Why this works better:

• Vectorized operations: We replaced sapply with direct matrix comparisons, which are implemented in C under the hood—way faster than looping through each value.
• Mask-based filling: Instead of plotting individual points, we use image() to overlay the entire grey region in a single render pass.
• Correct constraint logic: The constraint_mask function now accurately targets the region where both constraints are positive, matching your requirement.

Even Cleaner Solution with ggplot2 #

If you're open to using the ggplot2 package, it's built for efficient grid visualization, and the code becomes much more readable and flexible:

library(ggplot2)

# Create a data frame with all grid points
grid_df <- expand.grid(x1 = x1, x2 = x2)

# Calculate f values and the constraint mask in one go
grid_df <- transform(grid_df,
  f = x1 + x2,
  fill_grey = {
    c1 <- 1.5 - x1 - 2*x2 - 0.5*sin(2*pi*(x1^2 - 2*x2))
    c2 <- x1^2 + x2^2 - 1.5
    c1 > 0 & c2 > 0
  }
)

# Build the plot
ggplot(grid_df, aes(x = x1, y = x2)) +
  # Render the f values as a raster (fast, vectorized grid rendering)
  geom_raster(aes(fill = f)) +
  # Add contour lines for f
  geom_contour(aes(z = f), color = "black") +
  # Overlay the grey region (only where fill_grey is TRUE)
  geom_raster(aes(alpha = fill_grey), fill = "lightgrey") +
  # Set alpha to 1 for the grey region, 0 otherwise
  scale_alpha_manual(values = c("FALSE" = 0, "TRUE" = 1)) +
  # Customize labels and theme
  labs(x = expression(x[1]), y = expression(x[2])) +
  theme_minimal() +
  # Hide the alpha legend since we don't need it
  guides(alpha = "none")

Why this is great:

• geom_raster handles grid data in a single optimized operation—no more slow point plotting.
• The code is more readable, with clear separation of data preparation and plotting.
• You get easy access to color scales, themes, and other customization options without writing extra low-level code.

Both solutions will drastically improve the performance of your plot, especially with fine grid steps. The base R version is perfect if you want to stick to core functions, while ggplot2 offers a more modern, flexible approach.

内容的提问来源于stack exchange，提问作者John Smith