一个带有大量注解代码的实现。它假设A型和B型对象必须完全位于椭圆柱体中。

library(data.table)

rObj <- function(rx, ry, h, n, dims, eps = 2) {
  # Function to create a random sample (by rejection) of non-overlapping
  # rectangular prism objects inside an elliptical cylinder whose ellipse is
  # centered at x = 0, y = 0 and whose height ranges from -dims[3]/2 to h -
  # dims[3]/2. The objects have dimensions (x, y, z) = dims, and all edges are
  # parallel or orthogonal to each of the x, y, or z axes.
  # INPUTS:
  #   rx:   length of the ellipse
  #   ry:   width of the ellipse
  #   h:    height of the elliptical cylinder
  #   n:    number of non-overlapping objects to return
  #   dims: dimensions of the rectangular prism objects (vector of length 3)
  #   eps:  oversampling factor
  # OUTPUT: a data.table with 3 columns and n rows. Each row gives the
  #         coordinates of the centroid of a sampled object
  dt <- data.table()
  while(nrow(dt) < n) {
    # increase oversampling if it is not the first pass
    if (nrow(dt)) eps <- eps*2
    rho <- sqrt(runif(eps*n))
    phi <- runif(eps*n, 0, 2*pi)
    dt <- data.table(
      # sample object centroids
      # see https://stackoverflow.com/questions/5529148/algorithm-calculate-pseudo-random-point-inside-an-ellipse
      # First, uniformly sample on an ellipse centered on x = 0, y = 0,
      # with xlength = rx - dims[1] and ylength = ry - dims[2]
      # (any object with a centroid outside of this ellipse will stick out of
      # the elliptical cylinder, although some with a centroid within the
      # smaller ellipse will still stick out of the elliptical cylinder).
      x = (rx - dims[1])/2*rho*cos(phi),
      y = (ry - dims[2])/2*rho*sin(phi),
      # uniformly sample centroid heights
      z = runif(eps*n, 0, h - dims[3])
    )[
      # remove objects that stick out of bounds
      # The ellipse satisfies (x/(rx/2))^2 + (y/(ry/2))^2 = 1, which is the
      # same as (x/rx)^2 + (y/ry)^2 = 0.25. Taking advantage of symmetry, add
      # half of the x and y dimensions of the objects to the absolute value of
      # x and y (the object corner furthest from the foci of the ellipse) and
      # check if the result satisfies the standard equation.
      ((abs(x) + dims[1]/2)/rx)^2 + ((abs(y) + dims[2]/2)/ry)^2 < 0.25
    ][
      # remove objects that overlap a previously placed object
      # Since each rectangular prism object is oriented with the x, y, z axes,
      # two objects overlap if they are closer than their lengths in each
      # dimension.
      tabulate(
        sequence((.N - 1L):1, 2:.N)[ # row numbers (always keep the first row)
          (dist(x) < dims[1]) & (dist(y) < dims[2]) & (dist(z) < dims[3])
        ],
        .N
      ) == 0L
    ]
  }
  dt[1:n] # keep the first n objects
}

# function to get pairwise distances between two vectors
dist2 <- function(x, y) abs(outer(x, y, "-"))

fsim <- function(rx, ry, h, nA, nB, dimA, dimB, nreps, eps = 2) {
  # function to simulate placement of A and B rectangular prism objects inside
  # an elliptical cylinder and count the number of A-type objects that
  # intersect at least one B-type object. All object edges are parallel or
  # orthogonal to each of the x, y, or z axes.
  
  # INPUTS:
  #   rx:     length of the ellipse
  #   ry:     width of the ellipses
  #   h:      height of the elliptical cylinder
  #   nA:     number of non-overlapping A-type objects to return
  #   nB:     number of non-overlapping B-type objects to return
  #   dimX:   dimensions of the rectangular prism objects (vector of length 3)
  #   nreps:  the number of replications to simulate
  #   eps:    oversampling factor when randomly sampling non-overlapping objects
  #           by rejection
  # OUTPUT: vector of length "nreps" giving the number of A-type objects that
  # intersect at least one B-type object for each replication
  dims <- rowMeans(cbind(dimA, dimB)) # average dimensions of the A and B objects
  out <- integer(nreps) # initialize the output vector
  # repeat the simulation "nreps" times
  for (i in 1:nreps) {
    # get the coordinates of the A- and B-type objects' centroids
    A <- rObj(rx, ry, h, nA, dimA, eps)
    B <- rObj(rx, ry, h, nB, dimB, eps)
    # count the number of A-type objects that intersect at least one B-type
    # object
    out[i] <- sum(rowSums((dist2(A$x, B$x) < dims[1])*(dist2(A$y, B$y) < dims[2])*(dist2(A$z, B$z) < dims[3])) != 0L)
  }
  
  out
}

时间10K模拟复制：

system.time(overlaps <- fsim(9, 5, 2, 50L, 50L, rep(0.4, 3), rep(0.4, 3), 1e4L))
#>    user  system elapsed 
#>   27.19    0.25   27.67
mean(overlaps)
#> [1] 18.7408

一种近似求解这个问题的方法是将物体离散化。将一个体积设置为一个三维的零数组，然后随机地一次生成一个形状的参数。
对于每个生成的形状，查找将位于该形状内部的数组的所有元素。如果任何位置将位于圆柱体外部或与同一类型的形状重叠，请重试。获得法律的形状后，标记这些数组条目（例如，类型A为1，类型B为2）。首先执行所有类型A，然后执行所有类型B，并记录形状B占用先前为形状A标记的空间的次数。