Script to run single coverage experiment; based on Delta vs Alpha but with s=+25%

e5.array.R

#~ '
#~ using sim_replicate3 -- includes breakthrough & updated models
#~ Based on situation 17 sep 2021
#~ Including hier Bayes 
#~ smaller rho0
#~ s = +25%
#~ '

selcoef <- .25 
VAR_IMPORTS <- 650
rho1 = .66 

st0 = Sys.time() 

library( variantAnalysis ) # https://github.com/emvolz-phylodynamics/variantAnalysis

jobindex = as.numeric( commandArgs(trailingOnly=TRUE)[1] ) 

print( Sys.time() )

nreps = 100

tfins <- seq( lubridate::decimal_date( as.Date( "2021-09-15"))
			, lubridate::decimal_date( as.Date( "2021-11-15"))
			, by = 7/365 
			)
rho0s <-  c( .01, .05, .10, .20) / rho1 # note the stuff in c() is coverage of infections not coverage of cases  
repparms = expand.grid( tfins, rho0s, 1:nreps )


fn = system.file( 'stochastic_seir2.R', package = 'variantAnalysis' )
suppressMessages( {seir_gen = odin::odin( fn ) } )



parms <- unlist( repparms[ jobindex, ]  ) 

print( 'parms, job index: ' )
print( parms )
print( jobindex ) 

#~ source('/home/erikvolz/git/variantAnalysis/R/power1.R' )
st0 = Sys.time() 
res <- sim_replicate3( 
	  TFIN = parms[1] , RHO0 =  parms[2]
	  , rho1 = rho1
	  , MU = lubridate::decimal_date( as.Date( "2021-09-15"))
	  , seir_gen = seir_gen 
	  , E_imports0 = VAR_IMPORTS 
	  , E_imports1 = VAR_IMPORTS*2
	  , Rvariant = 1.05*(1+selcoef)
	  , Rancestral = 1.05 #
	  , k = 25
	  , ncpu = 6
	  , bs_replicates = 100
	)
st1 <- Sys.time() 

#~ # dbg
#~ variantAnalysis::variable_frequency_day( res$data, variable = 'lineage' , value = 'variant') -> vfd 
#~ print( vfd$plot )
#~ print( dim( res$data ) )
#~ print( res$fit_logistic )
#~ print( res$fit_importation_adjusted )
#~ print( res$fit_breakthrough )

print( Sys.time() )
print( st1 - st0 )

saveRDS( res, file = paste0('e5-res-', jobindex, '.rds' ) )