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

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

## ----eval=TRUE, echo=FALSE----------------------------------------------------
if (!has_rf) {
  message("Can't create vignette 'levelSets Example:  Classifiers' because ", 
          "package 'randomForest' is not available")
  knitr::knit_exit()
} else if (!has_peng) {
  message("Can't create vignette 'levelSets Example:  Classifiers' because ", 
          "package 'palmerpenguins' is not available")
  knitr::knit_exit()
}

## -----------------------------------------------------------------------------
penguin_data <- stats::na.omit(as.data.frame(palmerpenguins::penguins))
# Shorter names for feature (input) space dimensions:
names(penguin_data) <- sub("_length_mm", "_len", names(penguin_data))
names(penguin_data) <- sub("_depth_mm", "_dep", names(penguin_data))
names(penguin_data) <- sub("body_mass_g", "weight", names(penguin_data))
penguin_data$sex <- ifelse(penguin_data$sex == "female", 1, 2)  # now numeric
# Features to use for modeling:
pvars <- c("bill_len", "bill_dep", "flipper_len", "weight", "sex")

set.seed(1)
rf1 <- randomForest::randomForest(penguin_data[, pvars], 
                                  y=penguin_data[, "species"], 
                                  importance=TRUE)
print(rf1)
randomForest::varImpPlot(rf1)

## -----------------------------------------------------------------------------
respfn <- function(x, model) {
  if (nrow(x) > 0) {
    predprobs <- predict(model, newdata=x, type="prob")  # 3-column matrix
                                              # of class probabilities
    unname(predprobs[, "Adelie"])
  } else  numeric(0)  # 'predict.randomForest' does not like empty 'newdata'
}

## -----------------------------------------------------------------------------
feasbnds <- apply(data.matrix(penguin_data[, pvars]), 2, range)
library(levelSets)
fobj <- fnObj(pvars, respfn=respfn, feasbnds=feasbnds, model=rf1, 
              derivmethod="none")
d <- inpDim(fobj)  # 5

## -----------------------------------------------------------------------------
rect0 <- boundingRect(feasbnds, expand=1.1, spec=fobj)
scl0 <- inpScale(feasbnds[2, ] - feasbnds[1, ], spec=fobj)

## -----------------------------------------------------------------------------
thresh <- 0.75

## -----------------------------------------------------------------------------
rays0 <- axisRays(fobj, rect=rect0, scale=scl0)
profs <- respProfiles(fobj, rays=rays0, posns=c(-10:10)/10, 
                      lsetthresh=thresh)
summary(profs)

## ----eval=requireNamespace("ggplot2", quietly=TRUE)---------------------------
plot(profs, groupBy=inpNames(fobj))  # 'groupBy' puts each ray in a separate panel

## ----eval=!requireNamespace("ggplot2", quietly=TRUE), echo=FALSE--------------
# message("'ggplot2' package for plotting 'respProfiles' is not available")

## -----------------------------------------------------------------------------
profinfo <- respInfo(profs)  # data frame, including column 'lset'
print(subset(profinfo, lset))

## -----------------------------------------------------------------------------
srch1 <- bdrySearch(fobj, lsetthresh=thresh, rect=rect0, scale=scl0, 
                    initOrigin=ptCoord(profs)[8, ])
segs1 <- lsetSegs(srch1)

## -----------------------------------------------------------------------------
chk <- lsetSegsCheck(segs1, fnobj=fobj, posns=(1:4)/5)
table(chk$validIn)  # all TRUE

## -----------------------------------------------------------------------------
summary(segs1)
pairs(segs1, rect=rect0, main=paste0("Segments in feature space with ", 
                                     "Pr(Adelie) >= ", thresh))

## -----------------------------------------------------------------------------
slices2 <- sliceMat(list(sex=c(1, 2), weight=c(3000, 4500, 6000)), 
                    spec=fobj)  # 2 * 3 = 6 slices
srch2 <- slicedBdrySearch(fobj, lsetthresh=thresh, rect=rect0, scale=scl0, 
                          slices=slices2, currentBdry=srch1, 
                          control=list(axisDegree=2))

## -----------------------------------------------------------------------------
if (requireNamespace("ggplot2", quietly=TRUE)) {
  plt <- plotSegsList(srch2, dims=c("bill_len", "bill_dep"), rect=feasbnds, 
                      wrap=list(ncol=2))
  print(plt + ggplot2::labs(title=paste0("Pr(Adelie) >= ", thresh, 
                                         " (bill_dep vs bill_len)")))
  plt <- plotSegsList(srch2, dims=c("bill_len", "flipper_len"), rect=feasbnds, 
                      wrap=list(ncol=2))
  print(plt + ggplot2::labs(title=paste0("Pr(Adelie) >= ", thresh, 
                                         " (flipper_len vs bill_len)")))
  plt <- plotSegsList(srch2, dims=c("bill_dep", "flipper_len"), rect=feasbnds, 
                      wrap=list(ncol=2))
  print(plt + ggplot2::labs(title=paste0("Pr(Adelie) >= ", thresh, 
                                         " (flipper_len vs bill_dep)")))
}