# Example:  Response function is a combination of two quadratic surfaces in 
# a 2-dimensional input space.  The response surface has two modes.
#
# The level set of interest is the set of input points for which response 
# is greater than or equal to -40.  It consists of two separate elliptical 
# regions.
#
# Two versions of the example are created.  The first has no constraint on 
# the feasible region of the response function.  In the second version the 
# feasible region is artificially set to be points with X1 >= 0.

# This script creates and saves (as 'dbl_ellipse_2dEx.rda' and 
# 'dbl_ellipse_X1bnd_2dEx.rda' in directory '../data') lists with components:
#  fobj : 'fnObj' object with the response function, response and input  
#      space specifications, and feasible region function.
#  thresh : Scalar with the contour value defining the level set 
#      of interest (-40).
#  origin0 : 1-row matrix with the coordinates of a point (2, 1) 
#      in the level set to use as a ray origin.
#  rect0 : 'inpRect' object with an initial hyper-rectangle in input 
#      space, to guide the search for the level set boundary.
#  rays0 : 'inpRays' object containing 50 random rays with origin and 
#      rectangle given by 'origin0' and 'rect0'.
#  segs, segs_checked : 'lsetSegs' objects with level set segments for 
#      'fobj'.  'segs' was obtained by applying 'bdryFromRays()' to 'rays0'. 
#      Some of the segments are invalid, and 'segs_checked' contains 
#      the segments after being checked and fixed by 'lsetSegsCheck()' 
#      (all valid).

# 12 Jan 2026 / Richard Raubertas

local({
  respf <- function(x) {
     c1 <- x[, 1]^2 + 2*x[, 1]*x[, 2] + 2*(x[, 2]^2)
     # Same ellipse, but shifted to be centered at (8, 8):
     c2 <- (x[, 1]-8)^2 + 2*(x[, 1]-8)*(x[, 2]-8) + 2*((x[, 2]-8)^2)
     #resp <- -1 * pmin(c1, c2)
     resp <- ifelse(x[, 1] + 2*x[, 2] < 15, -1*c1, -1*c2)
     resp
  }
  fobj1 <- fnObj(c("X1", "X2"), respfn=respf)
  # Feasible region:  X1 >= 0.
  feasf <- function(x) { (x[, 1] >= 0) }
  fobj2 <- fnObj(c("X1", "X2"), respfn=respf, feasfn=feasf)
  
  thresh <- -40
  origin0 <- matrix(c(2, 1), nrow=1)
  rect0 <- inpRect(rbind(c(-5, -10), c(15, 12)), spec=fobj1)
  set.seed(1)
  rays0 <- randomRays(50, origin=origin0, rect=rect0)
  
  segs <- bdryFromRays(fobj1, rays=rays0, lsetthresh=thresh)
  segs_checked <- lsetSegsCheck(segs, fnobj=fobj1, posns=(1:9)/10)$segs
  dbl_ellipse_2dEx <- list(fobj=fobj1, thresh=thresh, origin0=origin0, 
                           rect0=rect0, rays0=rays0, segs=segs, 
                           segs_checked=segs_checked)


  segs <- bdryFromRays(fobj2, rays=rays0, lsetthresh=thresh)
  segs_checked <- lsetSegsCheck(segs, fnobj=fobj2, posns=(1:4)/5)$segs
  dbl_ellipse_X1bnd_2dEx <- list(fobj=fobj2, thresh=thresh, origin0=origin0, 
                                 rect0=rect0, rays0=rays0, segs=segs, 
                                 segs_checked=segs_checked)

  objs <- c("dbl_ellipse_2dEx", "dbl_ellipse_X1bnd_2dEx")
  for (obj in objs) {
    save(list=obj, file=paste0("../../data/", obj, ".rda"), compress=TRUE)
  }

  # Clean up.  (This may not be necessary due to 'local()' wrap.)
  rm(respf, fobj1, feasf, fobj2, thresh, origin0, rect0, rays0, segs, 
     segs_checked, dbl_ellipse_2dEx, dbl_ellipse_X1bnd_2dEx, obj, objs)
})

if (FALSE) {
  # Code like this can be inserted at the start of example sections of 
  # help pages for functions that work with 'lsetSegs' objects:
  Ex <- dbl_ellipse_X1bnd_2dEx  # see '?dbl_ellipse_2dEx'
  segs <- Ex$segs  # 8 invalid segs
  segs_checked <- Ex$segs_checked
  
  rect0 <- Ex$rect0
  plot(segs, rect=rect0, main="Segments found along 50 random rays")
  plot(segs_checked, rect=rect0, 
       main="Segments along 50 random rays, invalid segs fixed")

}
