## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 3.2,
  fig.height = 3.2,
  dpi = 96
)

## ----setup--------------------------------------------------------------------
library(rcicr)

# Graphical parameters are changed below to draw the images without margins;
# this records the originals so the last chunk can put them back.
old_par <- par(no.readonly = TRUE)

## ----show-helper--------------------------------------------------------------
# zlim matters more than it looks. image() stretches whatever range it is given
# across the full palette, so without a fixed zlim every linear rescaling of the
# same image renders identically -- which would make the scaling comparison
# below silently meaningless. Pass zlim = c(0, 1) whenever the point is what the
# pixel values actually are; the default just stretches to the data's own range,
# which is what you want when only the structure matters.
show <- function(m, title, zlim = range(m, na.rm = TRUE)) {
  op <- par(mar = c(0, 0, 1.4, 0))
  # image() takes [x, y] with y increasing upwards, so a matrix indexed
  # [row, col] has to be transposed and flipped to display the right way up.
  image(t(m[nrow(m):1, ]), col = gray.colors(256), axes = FALSE, asp = 1, # nolint: seq_linter.
        main = title, zlim = zlim, useRaster = TRUE)
  par(op)
}

## ----install, eval = FALSE----------------------------------------------------
# # install.packages("remotes")
# remotes::install_github("rdotsch/rcicr@*release")   # or @v1.2.3 for one specific release

## ----base-face, fig.width = 2.2, fig.height = 2.2-----------------------------
n <- 64
rows <- matrix(seq(-1, 1, length.out = n), n, n)
cols <- matrix(seq(-1, 1, length.out = n), n, n, byrow = TRUE)
x <- cols
y <- -rows # row indices grow downwards, so flip to make +y point up

face <- exp(-(x^2 / 0.45 + y^2 / 0.75))                            # head
face <- face - 0.55 * exp(-(((x + 0.3)^2 + (y - 0.25)^2) / 0.012)) # left eye
face <- face - 0.55 * exp(-(((x - 0.3)^2 + (y - 0.25)^2) / 0.012)) # right eye
face <- face - 0.35 * exp(-(x^2 / 0.10 + (y + 0.42)^2 / 0.006))    # mouth
face <- (face - min(face)) / (max(face) - min(face))

base_face <- tempfile(fileext = ".png")
png::writePNG(face, base_face)

show(face, "synthetic base face", zlim = c(0, 1))

## ----generate-stimuli, results = "hide"---------------------------------------
stimulus_path <- tempfile("stimuli")
dir.create(stimulus_path)

generateStimuli2IFC(
  base_face_files = list(face = base_face),
  n_trials        = 120,
  img_size        = 64,
  stimulus_path   = stimulus_path,
  seed            = 1,
  nscales         = 3,
  ncores          = 1,
  save_as_png     = FALSE # TRUE in a real study: this writes the actual stimuli
)

rdata_file <- list.files(stimulus_path, pattern = "\\.Rdata$", full.names = TRUE)[1]

## ----simulate-----------------------------------------------------------------
e <- new.env()
load(rdata_file, envir = e)
params <- e$stimuli_params[["face"]]

# The template is itself a noise image, so it is exactly expressible in the same
# basis the stimuli are drawn from.
set.seed(99)
template <- generateNoiseImage(rnorm(max(e$p$patchIdx)), e$p)

# Each trial's noise image, and how strongly it matches the template.
stack <- vapply(seq_len(nrow(params)),
                function(i) generateNoiseImage(params[i, ], e$p),
                matrix(0, 64, 64))
evidence <- apply(stack, 3, function(z) base::sum(z * template))
evidence <- evidence / sd(evidence)

# Three observers with the same template but increasing internal noise: the
# third is much less consistent than the first.
set.seed(7)
simulate <- function(internal_noise) {
  ifelse(evidence + rnorm(length(evidence), 0, internal_noise) > 0, 1, -1)
}

responses <- data.frame(
  participant = rep(c("p01", "p02", "p03"), each = nrow(params)),
  stimulus    = rep(seq_len(nrow(params)), 3),
  response    = c(simulate(0.5), simulate(1.5), simulate(3))
)

head(responses)

## ----one-ci-------------------------------------------------------------------
ci_p01 <- generateCI(
  stimuli     = responses$stimulus[responses$participant == "p01"],
  responses   = responses$response[responses$participant == "p01"],
  baseimage   = "face",
  rdata       = rdata_file,
  save_as_png = FALSE
)

names(ci_p01)

## ----show-recovery, fig.show = "hold", fig.width = 2.6, fig.height = 2.6------
show(template, "true template")
show(ci_p01$ci, "recovered CI")

## ----recovery-cor-------------------------------------------------------------
cor(as.vector(ci_p01$ci), as.vector(template))

## ----scaling-demo, fig.show = "hold", fig.width = 2.1, fig.height = 2.1-------
for (method in c("none", "constant", "matched", "independent")) {
  res <- generateCI(
    stimuli     = responses$stimulus[responses$participant == "p01"],
    responses   = responses$response[responses$participant == "p01"],
    baseimage   = "face", rdata = rdata_file, save_as_png = FALSE,
    scaling     = method, scaling_constant = 0.5
  )
  # $scaled, not $combined, so the effect of scaling is visible rather than
  # hidden under the base image -- and zlim fixed to the displayable range, so
  # that what you see is the actual pixel values.
  show(res$scaled, method, zlim = c(0, 1))
}

## ----batch, results = "hide"--------------------------------------------------
cis <- batchGenerateCI(
  data        = responses,
  by          = "participant",
  stimuli     = "stimulus",
  responses   = "response",
  baseimage   = "face",
  rdata       = rdata_file,
  save_as_png = FALSE
)

## ----batch-names--------------------------------------------------------------
names(cis)

## ----autoscale, fig.show = "hold", fig.width = 2.1, fig.height = 2.1----------
scaled <- autoscale(cis, save_as_pngs = FALSE)

for (nm in names(scaled)) {
  # $scaled, not $combined -- see the note below.
  show(scaled[[nm]]$scaled, sub(".*_", "", nm), zlim = c(0, 1))
}

## ----autoscale-combined, fig.width = 2.4, fig.height = 2.4--------------------
p01 <- scaled[["face_participant_p01"]]
show((p01$scaled + p01$base) / 2, "p01 over base", zlim = c(0, 1))

## ----infoval, eval = FALSE----------------------------------------------------
# # Slow: simulates `iter` classification images from random responses. Do this once
# # per stimulus set; the result is cached back into the .Rdata file.
# generateReferenceDistribution2IFC(rdata_file, iter = 10000)
# 
# computeInfoVal2IFC(target_ci = ci_p01, rdata = rdata_file)

## ----zmap, results = "hide", fig.width = 3.2, fig.height = 3.2----------------
zmap_dir <- tempfile("zmaps")

ci_z <- generateCI(
  stimuli     = responses$stimulus[responses$participant == "p01"],
  responses   = responses$response[responses$participant == "p01"],
  baseimage   = "face", rdata = rdata_file, save_as_png = FALSE,
  zmap = TRUE, zmapmethod = "quick", threshold = 1.5,
  zmaptargetpath = zmap_dir, zmapdecoration = FALSE
)

## ----zmap-show----------------------------------------------------------------
# Pixels that did not clear the threshold are set to NA.
range(ci_z$zmap, na.rm = TRUE)
mean(!is.na(ci_z$zmap)) # fraction of the image flagged

## ----citation, eval = FALSE---------------------------------------------------
# citation("rcicr")

## ----reset-par, include = FALSE-----------------------------------------------
par(old_par)