---
title: "Performance"
vignette: >
  %\VignetteIndexEntry{Performance}
  %\VignetteEngine{litedown::vignette}
  %\VignetteEncoding{UTF-8}
---



The sequential Poisson method, and order sampling methods more generally, are
simple and consequently not computationally expensive. This makes them suitable
for a range of different applications, especially when drawing a sample from a
large population. Despite this, there are
two optimizations in this package to keep drawing samples fast.

## Calculating inclusion probabilities

The first optimization concerns the inclusion probabilities when there are
take-all units. As seen in `vignette("take-all")`, the algorithm finds
take-all units one at a time without recomputing the inclusion probabilities
many times. This is much faster than the naive approach when drawing a large
sample, and is on par with the usual algorithm that finds take-all units in 
batches.

``` {.r}
library(sps)

# Make a population with 200 take-all units.
x <- c(rep(1, 1e6 - 200), rep(1e6, 200))
n <- 1e3


# Naive implementation.
ip <- function(x, n, alpha = 0.001) {
  p <- \(x, n) x * (n / sum(x))
  ta_units <- integer(0)
  pi <- p(x, n)
  max_ts <- which.max(pi)
  while (pi[max_ts] > 1 - alpha) {
    ta_units <- c(ta_units, max_ts)
    pi <- p(replace(pi, max_ts, 0), n - length(ta_units))
    max_ts <- which.max(pi)
  }
  replace(pi, ta_units, 1)
}

cols <- c("expression", "median", "mem_alloc", "n_itr")
 
bench::mark(
  inclusion_prob(x, n),
  ip(x, n),
  sampling::inclusionprobabilities(x, n)
)[cols]
```

```
#> # A tibble: 3 × 4
#>   expression                               median mem_alloc n_itr
#>   <bch:expr>                             <bch:tm> <bch:byt> <int>
#> 1 inclusion_prob(x, n)                    24.18ms    61.1MB    21
#> 2 ip(x, n)                                  1.41s       3GB     1
#> 3 sampling::inclusionprobabilities(x, n)   70.1ms     134MB     8
```

## Partial sorting

The second optimization is recognizing that both the computation of inclusion
probabilities and the sequential Poisson method can benefit from partial sorting
algorithms. In both cases, only the $n$ largest/smallest elements of a vector
are needed and, when installed, the `topn()` function from the **kit** package
is used to avoid a complete sort. This is a drop-in
replacement for `order(x)[1:n]` that can be faster
when `n` is much smaller than the length of `x`, and can have a modest impact on
performance when drawing a sample from a large population.

``` {.r}
options(sps.usekit = TRUE)
bench::mark(sps(x, n))[cols]
```

```
#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 sps(x, n)    75.5ms     111MB     6
```

``` {.r}
options(sps.usekit = FALSE)
bench::mark(sps(x, n))[cols]
```

```
#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 sps(x, n)     153ms     126MB     4
```

Partial sorting generally speeds up drawing smaller samples as well, but the
effect is not as large because the sequential Poisson method is already quick.

## Iteratively drawing a sample

Both optimizations above can be seen when iteratively drawing a sequential
Poisson sample as creating an iterator involves finding the point at which each
unit enters the take-all stratum. Although this is more expensive than drawing
a single sample, it quickly pays off because of the small cost to iterating over
units in the sample.

``` {.r}
bench::mark(s <- sps_iterator(x, n), min_iterations = 10)[cols]
```

```
#> # A tibble: 1 × 4
#>   expression                median mem_alloc n_itr
#>   <bch:expr>              <bch:tm> <bch:byt> <int>
#> 1 s <- sps_iterator(x, n)    931ms     336MB    10
```

``` {.r}
bench::mark(s())[cols]
```

```
#> # A tibble: 1 × 4
#>   expression   median mem_alloc n_itr
#>   <bch:expr> <bch:tm> <bch:byt> <int>
#> 1 s()          2.82ms    15.2MB    80
```

``` {.r}
bench::mark(sps(x, n))
```

```
#> # A tibble: 1 × 13
#>   expression      min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time
#>   <bch:expr> <bch:tm> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm>
#> 1 sps(x, n)     158ms  159ms      6.13     126MB     7.66     4     5      653ms
#> # ℹ 4 more variables: result <list>, memory <list>, time <list>, gc <list>
```