## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(echo=TRUE, fig.align="center")
has_ggplot <- requireNamespace("ggplot2", quietly=TRUE)

## ----echo=FALSE, eval=has_ggplot----------------------------------------------
ggplot2::theme_update(plot.background=ggplot2::element_rect(fill=NA))

## -----------------------------------------------------------------------------
library(levelSets)
Ex <- circuitFailure_3dEx
fobj <- Ex$fobj  # 'fnObj' object with the response function
theta_mle <- Ex$theta_mle
theta_vcov <- Ex$theta_vcov
thresh <- Ex$thresh  # threshold value of log-likelihood that defines an 
                     # asymptotic 95% conf region

## -----------------------------------------------------------------------------
fobj <- update(fobj, resptol=c(0, 1e-3))

## -----------------------------------------------------------------------------
rect1 <- inpRect(rbind(theta_mle - 3*sqrt(diag(theta_vcov)), 
                       theta_mle + 3*sqrt(diag(theta_vcov))), 
                 spec=fobj)
scl1 <- inpScale(theta_vcov, spec=fobj)
set.seed(1)
rays1 <- randomRays(100, rect=rect1, scale=scl1)
segs1 <- bdryFromRays(fobj, rays=rays1, lsetthresh=thresh)
summary(segs1)
pairs(segs1, rect=rect1, main="95% LR confidence region (100 rays)")

## -----------------------------------------------------------------------------
srch2 <- bdrySearch(fobj, lsetthresh=thresh, rect=rect1, scale=scl1, 
                    initOrigin=theta_mle, 
                    control=list(axisDegree=2, updateScale=2))
segs2 <- lsetSegs(srch2)
summary(segs2)
pairs(segs2, rect=srch2$rect, main="Adaptive ray origins (81 rays total)")

## -----------------------------------------------------------------------------
srch2b <- bdrySearch(fobj, lsetthresh=thresh, 
                     # Build on previous search results, including updated 
                     # search region and scaling:
                     currentBdry=srch2, rect=srch2$rect, scale=srch2$scale,
                     control=list(axisDegree=2, 
                                  originCrit="maxproj", 
                                  srchDirection=c(0, 0, 1), 
                                  maxSteps=3))
segs2b <- lsetSegs(srch2b)
summary(segs2b)  # Added 3 new origins, 27 new rays
pairs(segs2b, rect=srch2b$rect, main="Add rays to focus on large 'scale'")

## -----------------------------------------------------------------------------
slices <- sliceMat(list(scale=c(10, 30, 50, 70, 90, 110)), spec=fobj)

## -----------------------------------------------------------------------------
rect3 <- boundingRect(rect1, slices, expand=1.1)

## -----------------------------------------------------------------------------
srch3 <- slicedBdryFromRays(fobj, lsetthresh=thresh, slices=slices, 
                            rect=rect3, scale=scl1, 
                            currentBdry=segs1)  # list with one component per slice
slsegs3 <- lapply(srch3, lsetSegs)  # list of 'lsetSegs' objects
# Number of level set segments found per slice:
sapply(slsegs3, length)
if (requireNamespace("ggplot2", quietly=TRUE)) {
  plt <- plotSegsList(slsegs3, dims=c("dfrac", "shape"), rect=rect3)
  print(plt + ggplot2::theme(panel.background=ggplot2::element_blank()) + 
        ggplot2::labs(title="Confidence region sliced by 'scale'"))
} 

## -----------------------------------------------------------------------------
rect4 <- boundingRect(srch2$rect, slices, expand=1.1)
srch4 <- slicedBdrySearch(fobj, lsetthresh=thresh, slices=slices, 
                          currentBdry=srch2, rect=rect4, 
                          scale=srch2$scale, 
                          control=list(axisDegree=2, 
                                       maxSteps=6))
slsegs4 <- lapply(srch4, lsetSegs)  # list of 'lsetSegs' objects
# Number of level set segments found per slice:
sapply(slsegs4, length)

## -----------------------------------------------------------------------------
if (requireNamespace("ggplot2", quietly=TRUE)) {
  plt <- plotSegsList(slsegs4, dims=c("dfrac", "shape"), rect=rect3)
  print(plt + ggplot2::theme(panel.background=ggplot2::element_blank()) + 
        ggplot2::labs(title="Confidence region sliced by 'scale' (adaptive search)"))
}