# timing_script.R -- regenerate the TIMING table of rsobol_fast.R on the
# current machine.  Usage, from a directory containing rsobol_fast.R:
#     Rscript timing_script.R
# If the original rsobol.R and fiftysobol.col are also present in the
# working directory, the original is timed directly at m = 10 and m = 15
# and extrapolated (x 32 per unit of m, marked est) at m = 20, giving
# speedup factors; otherwise only the new code's times are printed.
# The output is a ready-to-paste block of R comments.
#
# Warning: the m = 20 cells generate 2^20 points and take up to about a
# minute each; the whole grid can take several minutes, longer if the
# original is present (its m = 15, s = 50 cells run near a minute each).

source("rsobol_fast.R")
invisible(rsobol(m = 3, s = 1))                  # warm the table parse
invisible(rsobol(m = 3, s = 1, rand = FALSE))    # warm the lookup table

haveorig <- file.exists("rsobol.R") && file.exists("fiftysobol.col")
if (haveorig) {
  orig <- new.env()
  sys.source("rsobol.R", envir = orig)
}

timeit <- function(f, reps) {
  ts <- replicate(reps, system.time(f())[3])
  median(ts)
}

mset <- c(10, 15, 20)
sset <- c(5, 10, 20, 50)
types <- c("nestu", "truncmato", "fullmato", "plain")

newt <- array(NA_real_, c(length(mset), length(sset), length(types)),
              dimnames = list(mset, sset, types))
oldt <- newt

for (mi in seq_along(mset)) for (si in seq_along(sset)) {
  m <- mset[mi]; s <- sset[si]
  reps <- if (m >= 20) 1 else 3
  newt[mi, si, "nestu"] <- timeit(function() rsobol(m = m, s = s, seed = 3), reps)
  newt[mi, si, "truncmato"] <- timeit(function() rsobol(m = m, s = s, seed = 3,
                                                        type = "truncmato"), reps)
  newt[mi, si, "fullmato"]  <- timeit(function() rsobol(m = m, s = s, seed = 3,
                                                        type = "fullmato"), reps)
  newt[mi, si, "plain"] <- timeit(function() rsobol(m = m, s = s,
                                                    rand = FALSE), reps)
  cat(sprintf("new  m=%2d s=%2d done\n", m, s)); flush(stdout())
  if (haveorig && m <= 15) {
    oldt[mi, si, "nestu"] <- timeit(function()
      orig$rsobol(fn = "fiftysobol.col", m = m, s = s, seed = 3), 1)
    oldt[mi, si, "truncmato"] <- timeit(function()
      orig$rsobol(fn = "fiftysobol.col", m = m, s = s, seed = 3,
                  type = "mato"), 1)   # the original's name for truncmato
    oldt[mi, si, "plain"] <- timeit(function()
      orig$rsobol(fn = "fiftysobol.col", m = m, s = s, rand = FALSE), 1)
    cat(sprintf("orig m=%2d s=%2d done\n", m, s)); flush(stdout())
  }
}
if (haveorig && all(c("15", "20") %in% dimnames(oldt)[[1]]))
  oldt["20", , ] <- oldt["15", , ] * 32   # extrapolate original: x 32 per 5 in m, marked est

fs <- function(x) if (x < 10) sprintf("%7.3f", x) else sprintf("%7.1f", x)
cat("\nPaste-ready table (machine: EDIT THIS to describe your hardware):\n\n")
cat("#            type=nestu       type=truncmato     type=fullmato    rand=FALSE\n")
cat("#  m   s    sec    speedup    sec    speedup    sec (no prior)   sec    speedup\n")
for (mi in seq_along(mset)) {
  for (si in seq_along(sset)) {
    row <- sprintf("# %2d %3d ", mset[mi], sset[si])
    for (ty in types) {
      nt <- newt[mi, si, ty]; ot <- oldt[mi, si, ty]
      tag <- if (is.na(ot)) "" else
        sprintf("(%.0fx%s)", ot / nt, if (mset[mi] == 20) " est" else "")
      row <- paste0(row, sprintf(" %s %-10s", fs(nt), tag))
    }
    cat(sub(" +$", "", row), "\n")
  }
  if (mset[mi] != 20) cat("#\n")
}