experiments/convergence_nwfix.R

In edahelsinki/cyclesampler: Randomization Algorithms for Large Sparse Matrices

## -----------------------------------------------------------------------------
## Perform the second convergence experiment
## 
##       Sample networks where both edge and vertex weights are allowed to vary,
##       subject to interval constraints.
## -----------------------------------------------------------------------------
## Example usage:
##
##        Rscript --vanilla convergence_nwfix.R 2000 100 1 datadir outputdir 0 1
##         
## -----------------------------------------------------------------------------

library(cyclesampler)
source("utilities.R")
source("init.R")

## -----------------------------------------------------------------------------
## Get command-line arguments
## -----------------------------------------------------------------------------
args <- commandArgs(trailingOnly = TRUE)

## -----------------------------------------------------------------------------
## Set the random seed
## -----------------------------------------------------------------------------
set.seed(42)

## -----------------------------------------------------------------------------
## Get the dataset
## -----------------------------------------------------------------------------
n_steps      <- as.numeric(args[1])                  ## number of total steps to take
stepsize     <- as.numeric(args[2])                  ## how often to save
ds           <- dataset_list[[as.numeric(args[3])]]  ## index of the dataset

datadir      <- args[4]                              ## data directory
outputdir    <- args[5]                              ## output directory

running_time <- as.numeric(args[6])                  ## estimate running time (Boolean)
convergence  <- as.numeric(args[7])                  ## estimate convergence (Boolean)
## -----------------------------------------------------------------------------

data <- read_data(name = ds$fname, basepath = datadir)

res  <- c(ds, get_data_properties(data))

out <- list()
out$t_normalize <- system.time( data <- normalizedata(data, scale_values = FALSE) )[3]

out$t_addselfloops <- system.time(
{
    w_tmp   <- get_node_weights(data)

    nw_min <- w_tmp * 0.9
    nw_max <- w_tmp * 1.1

    sl_tmp  <- addselfloops(data, A = nw_min, B = nw_max)
    data_sl <- rbind(data, sl_tmp$data)
})[3]

out$t_init <- system.time( X <- cyclesampler(data_sl, a = c(rep(min(data[, 3]), nrow(data)), sl_tmp$a), b = c(rep(max(data[, 3]), nrow(data)), sl_tmp$b)) )[3]

## -----------------------------------------------------------------------------
## -- running time
## -----------------------------------------------------------------------------
if (running_time) {
    out$t_sample         <- replicate(10, system.time( tmp  <- X$samplecycles2(n = X$getncycles()) )[3])
    out$t_sample_average <- mean(out$t_sample)
    saveRDS(out, file = paste0(outputdir, "/", ds$name, "_running_time_nwfix.rds"), compress = "xz")
}

## -----------------------------------------------------------------------------
## -- convergence
## save the dataset after having taken <stepsize> steps
## -----------------------------------------------------------------------------
if (convergence) {
    nullstate <- X$getstate()
    resvec    <- rep(NA, n_steps)

    attr(resvec, "n_steps")       <- n_steps
    attr(resvec, "stepsize")      <- stepsize
    attr(resvec, "ds")            <- ds
    attr(resvec, "time_start")    <- Sys.time()

    for (i in seq.int(n_steps / stepsize)) {
        a <- (i-1) * stepsize + 1
        b <- i * stepsize

        attr(resvec, "t0_current") <- Sys.time()

        resvec[a:b] <- estimate_convergence(X, s_orig = nullstate, n_samples = stepsize, ind = seq.int(nrow(data)))

        attr(resvec, "t1_current") <- Sys.time()
        attr(resvec, "time_stop")  <- Sys.time()

        saveRDS(resvec, file = paste0(outputdir, "/", ds$name, "_convergence_nwfix.rds"), compress = "xz")
    }
}
## -----------------------------------------------------------------------------