---
title: "Comparing `CIPerm` with 'naive' approach"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{naive}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

We run simple timing comparisons to show how our package's approach with Nguyen (2009) compares against a "naive" grid-based search approach to confidence intervals from permutation methods.

# Defining streamlined functions

We can use `CIPerm`'s `cint(dset(x, y))` directly,
and it already wins against the naive for-loop or array-based approaches.
But since `dset` and `cint` both have extra cruft built into them,
for an even "fairer" comparison we recreate the core of `cint(dset())` here
without all the extra variables and without passing copies of dataframes.

First we re-define the internal function `roundOrCeiling()`, then we define a streamlined `cint.nguyen()`:

```{r function-nguyen}
#### Define a streamlined function for Nguyen approach ####

# 1. use round(siglevel*num), if siglevel*num is basically an integer,
#    to within the default numerical tolerance of all.equal();
# 2. or otherwise use ceiling(siglevel*num) instead.
roundOrCeiling <- function(x) {
  ifelse(isTRUE(all.equal(round(x), x)), # is x==round(x) to numerical tolerance?
         round(x),
         ceiling(x))
}

cint.nguyen <- function(x, y, nmc = 10000, conf.level = 0.95) {
  # Two-tailed CIs only, for now

  sig <- 1 - conf.level

  # Code copied/modified from within CIPerm::dset()
  n <- length(x)
  m <- length(y)
  N <- n + m
  num <- choose(N, n) # number of possible combinations
  # Form a matrix where each column contains indices in new "group1" for that comb or perm
  if(nmc == 0 | num <= nmc) { # take all possible combinations
    dcombn1 <- utils::combn(1:N, n)
  } else { # use Monte Carlo sample of permutations, not all possible combinations
    dcombn1 <- replicate(nmc, sample(N, n))
    dcombn1[,1] <- 1:n # force the 1st "combination" to be original data order
    num <- nmc
  }
  # Form the equivalent matrix for indices in new "group2"
  dcombn2 <- apply(dcombn1, 2, function(x) setdiff(1:N, x))

  # Form the corresponding matrices of data values, not data indices
  combined <- c(x, y)
  group1_perm <- matrix(combined[dcombn1], nrow = n)
  group2_perm <- matrix(combined[dcombn2], nrow = m)
  # For each comb or perm, compute difference in group means, k, and w_{k,d}
  diffmean <- colMeans(group1_perm) - colMeans(group2_perm)
  k <- colSums(matrix(dcombn1 %in% ((n+1):N),
                      nrow = n))
  wkd <- (diffmean[1] - diffmean) / (k * (1/n + 1/m))

  # Code copied/modified from within CIPerm::cint()
  # Sort wkd values and find desired quantiles
  w.i <- sort(wkd, decreasing = FALSE, na.last = FALSE)
  siglevel <- (1 - conf.level)/2
  index <- roundOrCeiling(siglevel*num) - 1
  #   When dset's nmc leads us to use Monte Carlo sims,
  #   we may get some permutations equivalent to orig data
  #   i.e. we may get SEVERAL k=0 and therefore several w.i=NaN.
  nk0 <- sum(k == 0)
  # Start counting from (1+nk0)'th element of w.i
  # (not the 1st, which will always be 'NaN' since k[1] is 0)
  LB <- w.i[1 + nk0 + index]
  UB <- w.i[(num - index)]

  CI <- c(LB, UB)
  conf.achieved <- 1 - (2*(index+1) / num)
    message(paste0("Achieved conf. level: 1-2*(", index+1, "/", num, ")"))
    return(list(conf.int = CI,
                conf.level.achieved = conf.achieved))
}
```

Next, we define an equivalent function to implement the "naive" approach: 

* Choose a grid of possible values of delta = mu_X-mu_Y to try
* For each delta...
    - subtract delta off of the x's,
    - and carry out a test of H_0: mu_X=mu_Y on the resulting data
* Our confidence interval consists of the range of delta values for which H_0 is NOT rejected

```{r function-forloop}
#### Define a function for "naive" approach with for-loop ####

cint.naive.forloop <- function(x, y, deltas,
                               nmc = 10000,
                               conf.level = 0.95) {
  # Two-tailed CIs only, for now

  sig <- 1 - conf.level
  pvalmeans <- rep(1, length(deltas))

  # Code copied/modified from within CIPerm::dset()
  n <- length(x)
  m <- length(y)
  N <- n + m
  num <- choose(N, n) # number of possible combinations
  # Form a matrix where each column contains indices in new "group1" for that comb or perm
  if(nmc == 0 | num <= nmc) { # take all possible combinations
    dcombn1 <- utils::combn(1:N, n)
  } else { # use Monte Carlo sample of permutations, not all possible combinations
    dcombn1 <- replicate(nmc, sample(N, n))
    dcombn1[,1] <- 1:n # force the 1st "combination" to be original data order
    num <- nmc
  }
  # Form the equivalent matrix for indices in new "group2"
  dcombn2 <- apply(dcombn1, 2, function(x) setdiff(1:N, x))

  for(dd in 1:length(deltas)) {
    xtmp <- x - deltas[dd]

    # Code copied/modified from within CIPerm::dset()
    combined <- c(xtmp, y)
    # Form the corresponding matrices of data values, not data indices
    group1_perm <- matrix(combined[dcombn1], nrow = n)
    group2_perm <- matrix(combined[dcombn2], nrow = m)
    # For each comb or perm, compute difference in group means
    diffmean <- colMeans(group1_perm) - colMeans(group2_perm)

    # Code copied/modified from within CIPerm::pval()
    pvalmeans[dd] <- sum(abs(diffmean - mean(diffmean)) >= abs(diffmean[1] - mean(diffmean)))/length(diffmean)
  }
  print( range(deltas[which(pvalmeans >= sig)]) )
  return(list(cint = range(deltas[which(pvalmeans >= sig)]),
              pvalmeans = pvalmeans,
              deltas = deltas))
}
```

We also wrote another "naive" function that avoids for-loops and takes a "vectorized" approach instead. We created large arrays: each permutation requires a matrix, and the 3rd array dimension is over permutations. Then `apply` and `colSums` work very quickly to carry out all the per-permutations steps on this array at once. However, we found that unless datasets were trivially small, creating and storing the array took a lot of time and memory, and it ended up far slower than the for-loop approach. Consequently we do not report this function or its results here, although curious readers can find it hidden in the vignette's .Rmd file.

```{r function-array, echo=FALSE, eval=FALSE}
## THIS FUNCTION ALSO WORKS,
## BUT REQUIRES SUBSTANTIALLY MORE MEMORY THAN THE PREVIOUS ONE
## AND IS SLOWER FOR LARGE DATASETS,
## SO WE LEFT IT OUT OF THE COMPARISONS

# Possible speedup:
# Could we replace for-loop with something like manipulating a 3D array?
# e.g., where we say
###     group1_perm <- matrix(combined[dcombn1], nrow = n)
# could we first make `combined` into a 2D matrix
#   where each col is like it is now,
#   but each row is for a different value of delta;
# and then group1_perm would be a 3D array,
#   where each 1st+2nd dims matrix is like it is now,
#   but 3rd dim indexes over deltas;
# and then `diffmean` and `pvalmeans`
#   would be colMeans or similar over an array
#   so the output would be right

cint.naive.array <- function(x, y, deltas,
                             nmc = 10000,
                             conf.level = 0.95) {
  # Two-tailed CIs only, for now
  sig <- 1 - conf.level

  # New version where xtmp is a matrix, not a vector;
  # and some stuff below will be arrays, not matrices...
  xtmp <- matrix(x, nrow = length(x), ncol = length(deltas))
  xtmp <- xtmp - deltas[col(xtmp)]
  combined <- rbind(xtmp,
                    matrix(y, nrow = length(y), ncol = length(deltas)))

  # Code copied/modified from within CIPerm::dset()
  n <- length(x)
  m <- length(y)
  N <- n + m
  num <- choose(N, n) # number of possible combinations
  # Form a matrix where each column contains indices in new "group1" for that comb or perm
  if(nmc == 0 | num <= nmc) { # take all possible combinations
    dcombn1 <- utils::combn(1:N, n)
  } else { # use Monte Carlo sample of permutations, not all possible combinations
    dcombn1 <- replicate(nmc, sample(N, n))
    dcombn1[,1] <- 1:n # force the 1st "combination" to be original data order
    num <- nmc
  }
  # Form the equivalent matrix for indices in new "group2"
  dcombn2 <- apply(dcombn1, 2, function(x) setdiff(1:N, x))

  # ARRAYS of data indices, where 3rd dim indexes over deltas
  dcombn1_arr <- outer(dcombn1, N * (0:(length(deltas)-1)), FUN = "+")
  dcombn2_arr <- outer(dcombn2, N * (0:(length(deltas)-1)), FUN = "+")

  # Form the corresponding ARRAYS of data values, not data indices
  group1_perm <- array(combined[dcombn1_arr], dim = dim(dcombn1_arr))
  group2_perm <- array(combined[dcombn2_arr], dim = dim(dcombn2_arr))

  # For each comb or perm, compute difference in group means
  # diffmean <- matrix(colMeans(group1_perm, dim=1), nrow = num) -
  #   matrix(colMeans(group2_perm, dim=1), nrow = num)
  diffmean <- colMeans(group1_perm, dim=1) - colMeans(group2_perm, dim=1)

  # Code copied/modified from within CIPerm::pval()
  pvalmeans <- colSums(abs(diffmean - colMeans(diffmean)[col(diffmean)]) >= abs(diffmean[1,] - colMeans(diffmean))[col(diffmean)])/nrow(diffmean)
  # plot(deltas, pvalmeans)

  print( range(deltas[which(pvalmeans >= sig)]) )
  return(list(cint = range(deltas[which(pvalmeans >= sig)]),
              pvalmeans = pvalmeans,
              deltas = deltas))
}
```

# Speed tests

## Tiny dataset

When we compare the timings, `cint.naive()` will need to have a reasonable search grid for values of `delta`. On this problem, we happen to know the correct CI endpoints are the integers `(-21, 3)`, so we "cheat" by using `(-22):4` as the grid for `cint.naive()`. Of course in practice, if you had only the naive approach, you would probably have to try a wider grid, since you wouldn't already know the answer.

```{r tests-tiny}
#### Speed tests on Nguyen's tiny dataset ####

# Use 1st tiny dataset from Nguyen's paper
library(CIPerm)
x <- c(19, 22, 25, 26)
y <- c(23, 33, 40)

# Actual CIPerm package's approach:
system.time({
  print( cint(dset(x, y), conf.level = 0.95, tail = "Two") )
})

# Streamlined version of CIPerm approach:
system.time({
  print( cint.nguyen(x, y, conf.level = 0.95) )
})

# Naive approach with for-loops:
deltas <- ((-22):4)
system.time({
  pvalmeans <- cint.naive.forloop(x, y, deltas)$pvalmeans
})

# Sanity check to debug `cint.naive`:
# are the p-vals always higher when closer to middle of CI?
# Are they always above 0.05 inside and below it outside CI?
cbind(deltas, pvalmeans)
plot(deltas, pvalmeans)
abline(h = 0.05)
abline(v = c(-21, 3), lty = 2)
# Yes, it's as it should be :)
```

Above we see timings for the original `CIPerm::cint(dset))`, its streamlined version`cint.nguyen()`, and the naive equivalent `cint.naive()`.  
In the case of a tiny dataset like this, all 3 approaches take almost no time. So it can be a tossup as to which one is fastest on this example.


```{r test-tiny-array, echo=FALSE, eval=FALSE}
# Naive approach with arrays:
system.time({
  pvalmeans <- cint.naive.array(x, y, deltas)$pvalmeans
})

# Sanity check again:
cbind(deltas, pvalmeans)
plot(deltas, pvalmeans)
abline(h = 0.05)
abline(v = c(-21, 3), lty = 2)
# Yep, still same results. OK.
# So this "array" version IS faster,
# and nearly as fast as Nguyen
# (though harder to debug...
#  but the same could be said of Nguyen's approach...)

# WAIT WAIT WAIT
# I changed the for-loop approach to stop creating dcombn1 & dcombn2
# inside each step of the loop
# since we can reuse the same ones for every value of deltas...
# AND NOW it's FASTER than the array version,
# even for tiny datasets!
# And later we'll see that array version is unbearably slow
# for large datasets where it takes forever
# to create & store these massive arrays...
# So it's not worth reporting in the final vignette :(
```

## Larger dataset

Try a slightly larger dataset, where the total number of permutations is above the default `nmc=10000` but still manageable on a laptop.

Again, `cint.naive()` will need to have a reasonable search grid for values of `delta`. This time again we will set up a grid that just barely covers the correct endpoints from `cint.nguyen()`, and again we'll try to keep it at around 20 to 25 grid points in all. Then we'll try timing it again with a slightly finer grid.

```{r tests-larger}
choose(18, 9) ## 48620
set.seed(20220528)
x <- rnorm(9, mean = 0, sd = 1)
y <- rnorm(9, mean = 1, sd = 1)

# Actual CIPerm approach
# (with nmc = 0,
#  so that we use ALL of the choose(N,n) combinations
#  instead of a MC sample of permutations)
system.time({
  print( cint(dset(x, y, nmc = 0),
              conf.level = 0.95, tail = "Two") )
})

# Streamlined version of CIPerm approach:
system.time({
  print( cint.nguyen(x, y, nmc = 0, conf.level = 0.95) )
})

# Coarser grid
deltas <- ((-21):(1))/10     # grid steps of 0.1

# Naive with for-loops:
system.time({
  pvalmeans <- cint.naive.forloop(x, y, deltas, nmc = 0)$pvalmeans
})


# Finer grid
deltas <- ((-21*2):(1*2))/20 # grid steps of 0.05

# Naive with for-loops:
system.time({
  pvalmeans <- cint.naive.forloop(x, y, deltas, nmc = 0)$pvalmeans
})
```

Now, `CIPerm::cint(dset())` and the `cint.nguyen` approach take about the same amount of time, and they are noticeably faster than the `cint.naive` approach. (And even more so when we use a finer search grid.) That happens even with the "cheat" of using `CIPerm` first to find the right CI limits so that our naive search grid isn't too wide, just stepping over the minimal required range with a reasonable grid-coarseness.

(For the original data where CI = (-21, 3), we stepped from -22 to 4 in integers.
For latest data where CI = (-2.06, 0.08),
we stepped from -2.1 to 0.1 in units of 0.1 and then 0.05.)

In practice there are probably cleverer ways for the naive method to choose grid discreteness and endpoints...
but still, this seems like a more-than-fair chance for the naive approach.


```{r tests-larger-array, echo=FALSE, eval=FALSE}
# Coarser grid
deltas <- ((-21):(1))/10     # grid steps of 0.1

# Naive with arrays:
system.time({
  pvalmeans <- cint.naive.array(x, y, deltas, nmc = 1e6)$pvalmeans
})

# !!! Indeed it looks like naive.forloop
# is FASTER than naive.array
# now that I've removed the slow dcombn1 & dcombn2 creating from the loop.


# Finer grid
deltas <- ((-21*2):(1*2))/20 # grid steps of 0.05

# Naive with arrays:
system.time({
  pvalmeans <- cint.naive.array(x, y, deltas, nmc = 1e6)$pvalmeans
})
```


## Largest dataset

Try an even larger example,
where we're OK with a smallish total number of permutations (eg `nmc=10000`),
but the dataset itself is "huge",
so that each individual permutation takes a longer time.

We still don't make it all that large, since we want this vignette to knit in a reasonable time. But we've seen similar results when we made even-larger datasets that took longer to run.


Once again, for `cint.naive()` we will set up a search grid that just barely covers the correct endpoints from `cint.nguyen()`, and again we'll try to keep it at around 20 to 25 grid points in all.

```{r tests-largest}
#### Speed tests on much larger dataset ####

set.seed(20220528)
x <- rnorm(5e3, mean = 0, sd = 1)
y <- rnorm(5e3, mean = 1, sd = 1)

# Actual CIPerm approach
# (with nmc = 2000 << choose(N,n),
#  so it only takes a MC sample of permutations
#  instead of running all possible combinations)
system.time({
  print( cint(dset(x, y, nmc = 2e3),
              conf.level = 0.95, tail = "Two") )
})

# Streamlined version of CIPerm approach:
system.time({
  print( cint.nguyen(x, y, nmc = 2e3, conf.level = 0.95) )
})

# Grid of around 20ish steps
deltas <- ((-11*10):(-9*10))/100 # grid steps of 0.01

# Naive with for-loops:
system.time({
  pvalmeans <- cint.naive.forloop(x, y, deltas, nmc = 2e3)$pvalmeans
})
```

```{r tests-largest-array, echo=FALSE, eval=FALSE}
# Naive with arrays:
system.time({
  pvalmeans <- cint.naive.array(x, y, deltas, nmc = 2e3)$pvalmeans
})
```

Here we finally see that the overhead built into `CIPerm::cint(dset())` makes it slightly slower than the streamlined `cint.nguyen()`. However, both are *substantially* faster than `cint.naive()`---they run in less than half the time of the naive approach.


# `bench::mark()` and `profvis`

The results above are from one run per method on each example dataset.

Using the last ("largest") dataset, we also tried using `bench::mark()` to summarize the distribution of runtimes over many runs, as well as memory usage.  
*(Not actually run when the vignette knits, in order to avoid needing `bench` as a dependency.)*

```{r tests-largest-benchmark}
# bench::mark(cint.nguyen(x, y, nmc = 2e3),
#             cint.naive.forloop(x, y, deltas, nmc = 2e3),
#             check = FALSE, min_iterations = 10)
#> # A tibble: 2 × 13
#>   expression                                        min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory                  time            gc               
#>   <bch:expr>                                   <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>                  <list>          <list>           
#> 1 cint.nguyen(x, y, nmc = 2000)                    2.5s    2.71s     0.364    1.58GB     2.92    10    80     27.44s <NULL> <Rprofmem [32,057 × 3]> <bench_tm [10]> <tibble [10 × 3]>
#> 2 cint.naive.forloop(x, y, deltas, nmc = 2000)    7.87s    8.07s     0.120     7.4GB     4.64    10   388      1.39m <NULL> <Rprofmem [32,256 × 3]> <bench_tm [10]> <tibble [10 × 3]>

```

Finally, we also tried using the `profvis` package to check what steps are taking the longest time.  
*(Again, not actually run when the vignette knits, in order to avoid needing `profvis` as a dependency; but screenshots are included below.)*

```{r tests-largest-profvis}
# profvis::profvis(cint.nguyen(x, y, nmc = 2e3))
# profvis::profvis(cint.naive.forloop(x, y, deltas, nmc = 2e3))
```

For Nguyen's method, the initial setup with `combn()` or `sample()` followed by `apply(setdiff())` takes about 80% of the time, and the per-permutation calculations only take about 20% of the time:  
`profvis::profvis(cint.nguyen(x, y, nmc = 2e3))`

![](profvis_nguyen.png)

For the naive method, the initial setup is identical, and each round of per-permutation calculations is slightly faster than in Nguyen's method... but because you have to repeat them many times (*unless you already know the CI endpoints, in which case you wouldn't need to run this at all!*), they can add up to substantially more total time than for Nguyen's method.  
`profvis::profvis(cint.naive.forloop(x, y, deltas, nmc = 2e3))`

![](profvis_naive.png)
