How to generate an elliptical cylinder, populate it with randomly distributed points, and measure instances of overlap between those points in R

Question

I would like to determine the probability that a randomly distributed object of Type A occupies or touches (overlaps) the same space as any randomly distributed object of Type B when populated inside an elliptical cylinder. I would then like to loop this simulation many times to generate a more reliable probability value.

I am able to draw the elliptical cylinder using the shape package:

library(shape)
emptyplot(c(-5, 5), c(-15, 15), main = "filled elliptic cylinder")
filledcylinder(rx = 9, ry = 5, len= 2, angle = 00, col = "white",  
               lcol = "black", lcolint = "grey") 

I do not know how to add points (i.e. objects A and B) to this graph. However, I suspect graphical expression is not the way to go with this task (though I find visualising helpful). I suspect a better approach will be to create a function to describe the elliptical cylinder, similar to the cone in the following example, and run the simulation without graphical output:

# Create a function to describe a cone
cone <- function(x, y){
  sqrt(x ^ 2 + y ^ 2)
}

# prepare variables.
x <- y <- seq(-1, 1, length = 30)
z <- outer(x, y, cone)

# plot as a 3D surface for visual reference (even though I actually want a volume)
persp(x, y, z,
      main="Perspective Plot of a Cone",
      zlab = "Height",
      theta = 30, phi = 15,
      col = "orange", shade = 0.4)

Sadly I do not know how to do this for my elliptical cylinder. I am aware of the paramaters for describing an elliptical cylinder from the following source:

https://mathworld.wolfram.com/EllipticCylinder.html

Unfortunately, I do not understand much of it. I hope the dimensions given in my filledcylinder can act as a guide. Ultimately the dimension values do not matter, what matters is the code structure into which values can be entered.

As for the objects:

Let there be 50 Type A objects and 50 Type B objects of size x=0.4, y=0.4, z=0.4 (same units as in my graphical elliptical cylinder example).

All objects are to be distributed at random within the volume of the elliptical cylinder, with the exception that objects of Type A cannot overlap with another object of Type A, and objects of Type B cannot overlap with other objects of Type B. Type A objects may overlap with Type B objects.

I would like to output the number of Type A objects that overlap with any Type B object in the given volume, this number as a percentage of total Type A objects, and as a percentage of total all objects for each run of the simulation.

I do not know how to even start to do this.

If you can help, I'm afraid statistics, geometry and non-basic R expressions will need to be explained as if to a (not particularly bright) child.

Thank you very very much for your time!

Answer 1

An implementation with heavily commented code for explanations. This assumes the A- and B-type objects must be entirely within the elliptical cylinder.

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
}

Time 10K simulation replications:

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

Answer 2

One approach to get an approximate answer to this problem is to discretize things. Set up a volume as a 3 dimensional array of zeros, then randomly generate the parameters of your shapes one at a time.

For each generated shape, find all the elements of the array that would be inside the shape. If any locations would be outside the cylinder or overlap a shape of the same type, try again. Once you have a legal shape, mark those array entries (e.g. 1 for type A, 2 for type B). Do all type A first, then all type B, and keep count of the times when shape B occupies a space that was previously marked for shape A.