Multi-RNG is not random number safe

[cliffckerr]

Due to the use of ss.uniform.rvs() in _make_new_cases_multi_RNG():

import starsim as ss
import sciris as sc
import numpy as np

class use_np_rands(ss.Intervention):
    def apply(self, sim):
        print(np.random.rand())

class use_ss_rands(ss.Intervention):
    def apply(self, sim):
        print(ss.uniform.rvs(size=1))

pars = dict(diseases='sir', networks='random')

sims = [
    ss.Sim(pars=pars),
    ss.Sim(pars=pars, interventions=use_np_rands()),
    ss.Sim(pars=pars, interventions=use_ss_rands()),
]

for sim in sims:
    sim.run()

for i,sim in enumerate(sims):
    sc.heading(f'Sim {i}')
    print(sim.summarize())

——————————
Sim 0
——————————

#0. 'n_alive':            7463.694444444444
#1. 'new_deaths':         74.5
#2. 'sir_n_susceptible':  344.1388888888889
#3. 'sir_n_infected':     246.33333333333334
#4. 'sir_n_recovered':    6873.222222222223
#5. 'sir_prevalence':     0.02730965008537685
#6. 'sir_new_infections': 275.05555555555554
#7. 'sir_cum_infections': 9380.805555555555

——————————
Sim 1
——————————

#0. 'n_alive':            7432.416666666667 # eek, should match sim 0
#1. 'new_deaths':         75.44444444444444
#2. 'sir_n_susceptible':  360.6111111111111
#3. 'sir_n_infected':     245.88888888888889
#4. 'sir_n_recovered':    6825.916666666667
#5. 'sir_prevalence':     0.027285846815105312
#6. 'sir_new_infections': 274.6111111111111
#7. 'sir_cum_infections': 9364.777777777777

——————————
Sim 2
——————————

#0. 'n_alive':            7432.416666666667 # np.random.rand() and ss.uniform.rvs() give the same results
#1. 'new_deaths':         75.44444444444444
#2. 'sir_n_susceptible':  360.6111111111111
#3. 'sir_n_infected':     245.88888888888889
#4. 'sir_n_recovered':    6825.916666666667
#5. 'sir_prevalence':     0.027285846815105312
#6. 'sir_new_infections': 274.6111111111111
#7. 'sir_cum_infections': 9364.777777777777

[daniel-klein]

@cliffckerr , I believe this is because multirng is disabled by default. Adding the following line to the script:

ss.options(multirng=True)

resulted in the following output

——————————
Sim 0
——————————

#0. 'n_alive':            8169.583333333333
#1. 'new_deaths':         53.5
#2. 'sir_n_susceptible':  587.5555555555555
#3. 'sir_n_infected':     239.86111111111111
#4. 'sir_n_recovered':    7342.166666666667
#5. 'sir_prevalence':     0.02612392034310466
#6. 'sir_new_infections': 268.1388888888889
#7. 'sir_cum_infections': 9144.305555555555

——————————
Sim 1
——————————

#0. 'n_alive':            8169.583333333333
#1. 'new_deaths':         53.5
#2. 'sir_n_susceptible':  587.5555555555555
#3. 'sir_n_infected':     239.86111111111111
#4. 'sir_n_recovered':    7342.166666666667
#5. 'sir_prevalence':     0.02612392034310466
#6. 'sir_new_infections': 268.1388888888889
#7. 'sir_cum_infections': 9144.305555555555

——————————
Sim 2
——————————

#0. 'n_alive':            8169.583333333333
#1. 'new_deaths':         53.5
#2. 'sir_n_susceptible':  587.5555555555555
#3. 'sir_n_infected':     239.86111111111111
#4. 'sir_n_recovered':    7342.166666666667
#5. 'sir_prevalence':     0.02612392034310466
#6. 'sir_new_infections': 268.1388888888889
#7. 'sir_cum_infections': 9144.305555555555

[cliffckerr]

Ah, sorry, it was a poorly constructed example! It is not random number safe, but it only appeared that it was because in the default example beta=1, so the random numbers on infection actually didn't matter at all :joy:

Here's a better example:

import starsim as ss
import sciris as sc
import numpy as np

ss.options(multirng=True)

class use_np_rands(ss.Intervention):
    def apply(self, sim):
        print(np.random.rand())

class use_ss_rands(ss.Intervention):
    def apply(self, sim):
        print(ss.uniform.rvs(size=1))

pars = dict(diseases=dict(type='sir', beta=dict(random=[0.05, 0.05])), networks='random', n_agents=10e3)

sims = [
    ss.Sim(pars=pars, label='a'),
    ss.Sim(pars=pars, label='b', interventions=use_np_rands()),
    ss.Sim(pars=pars, label='c', interventions=use_ss_rands()),
]

for sim in sims:
    sim.run()

for i,sim in enumerate(sims):
    sc.heading(f'Sim {i}')
    print(sim.summarize())

which gives me

——————————
Sim 0
——————————

#0. 'n_alive':            9739.638888888889
#1. 'new_deaths':         7.277777777777778
#2. 'sir_n_susceptible':  8793.222222222223
#3. 'sir_n_infected':     5.277777777777778
#4. 'sir_n_recovered':    941.1388888888889
#5. 'sir_prevalence':     0.0005398673333522812
#6. 'sir_new_infections': 33.55555555555556
#7. 'sir_cum_infections': 1173.2222222222222

——————————
Sim 1
——————————

#0. 'n_alive':            9739.416666666666
#1. 'new_deaths':         7.277777777777778
#2. 'sir_n_susceptible':  8775.944444444445
#3. 'sir_n_infected':     5.75
#4. 'sir_n_recovered':    957.7222222222222
#5. 'sir_prevalence':     0.0005881118385117705
#6. 'sir_new_infections': 34.02777777777778
#7. 'sir_cum_infections': 1190.0277777777778

——————————
Sim 2
——————————

#0. 'n_alive':            9739.416666666666
#1. 'new_deaths':         7.277777777777778
#2. 'sir_n_susceptible':  8775.944444444445
#3. 'sir_n_infected':     5.75
#4. 'sir_n_recovered':    957.7222222222222
#5. 'sir_prevalence':     0.0005881118385117705
#6. 'sir_new_infections': 34.02777777777778
#7. 'sir_cum_infections': 1190.0277777777778

They're very similar but they don't quite match (1173 vs 1190 infections).

[daniel-klein]

Thanks, will take a second look.

[cliffckerr]

Fixed!