Repository: starsimhub/starsim (language: python)

Issue:
Modulus-based random numbers induces unwanted correlations
We use r = (r1 + r2) % 1 for transmission, where r1 is a random number from the first agent and r2 is a random number from the second agent. r is indeed uniformly distributed on the interval from [0,1). However, when reusing the same per-agent random draws to test transmission on multiple network edges, we encounter correlations. This makes intuitive sense because a network can have order N^2 edges, but we're only drawing order N random numbers.

The correlation can be seen with as few as N=4 nodes:
```python
import numpy as np

N = 4
p = 0.5
R = 1000000

r = np.random.rand(4, R)
for method in ['indep', 'modulo']:
    if method=='modulo':
        p01 = (r[0,:] + r[1,:]) % 1
        p02 = (r[0,:] + r[2,:]) % 1
        p23 = (r[2,:] + r[3,:]) % 1
        p13 = (r[1,:] + r[3,:]) % 1
    else:
        p01 = r[0,:]
        p02 = r[1,:]
        p23 = r[2,:]
        p13 = r[3,:]

    alone = np.count_nonzero(p13<p) / R
    combined = np.count_nonzero((p13<p) & (p01<p) & (p02<p) & (p23<p)) / np.count_nonzero((p01<p) & (p02<p) & (p23<p))
    print(method, alone, combined)
```

Output shows the bias:
```bash
indep 0.50018 0.5008594843094144
modulo 0.49979 0.668084936217859
```
Revert transmission to "acquisition-based" from "modulo"
### Description
Reverting transmission. Performance looks VERY BAD in the benchmark, 10x more time, but a practical application in TBsim changed from 107s in main to 110s with this branch. Not sure what's going on with the benchmark (or TBsim), but otherwise happy with the changes. Because this part of the model is heavily used, we might consider number or even cython for further performance optimization.

### Checklist
- [ ] Code commented & docstrings added
- [ ] New tests were needed and have been added
- [ ] A new version number was needed & changelog has been updated
- [ ] A new PyPI version needs to be released
Revert transmission example
Illustrative PR -- ready for review, but do not merge!! -- showing a different way of implementing the "roulette" approach (that actually does away with the roulette part). The thinking here is that most of the time all we really need is to know the probability of infection per target agent, which isn't too expensive to calculate (but still needs to be numbafied, since I couldn't think of a way to do it without a for loop or `np.apply_along_axis()`). The source determination could be done if and only if the user requests the infection log as well, since I don't think it's used/stored otherwise (afaik?).

Not super pretty. Also, for some reason this implementation randomly segfaults sometimes 🤣 Still hoping we can get to the bottom of this properly!


Discussion / hints:
I propose that we restore the previous node-based acquisition method.
I think this is a tricky issue. There are three things I'm still unclear on:
1. Whether this would be an issue in practice for large networks;
2. Whether the node-based acquisition method actually solves the issue;
3. How this issue relates to constructing CRN-safe networks in the first place.

### 1. Does it matter in practice?
The induced correlations between pairs of agents are proportional to 1/√n. With 1000 agents, they're on the order of 1-3%:
![image](https://github.com/starsimhub/starsim/assets/3239256/0b51aa3e-1ab1-4c5b-9102-1d2f4d1d5d9d)
```py
import numpy as np
import sciris as sc
import pylab as pl

n = 1000

r1 = np.random.rand(n,n)
v2a = np.random.rand(n)
v2b = np.random.rand(n).reshape((n,1))
r2a = np.tile(v2a, (n,1))
r2b = np.tile(v2b, (1,n))
r2 = (r2a + r2b) % 1.0

sc.options(dpi=200)
pl.figure(figsize=(12,8))
for i,r in enumerate([r1, r2]):

    mean0 = r.mean(axis=0)
    mean1 = r.mean(axis=1)
    std0 = mean0.std()
    std1 = mean1.std()
    
    if i == 0:
        x0 = x1 = np.arange(n)
        title = 'Truly random'
        xlabel = 'Index'
    else:
        x0 = v2a
        x1 = v2b
        title = 'Modulus approach'
        xlabel = 'Value'
    
    pl.subplot(2,3,1+i*3)
    pl.imshow(r)
    pl.title(title)
    
    pl.subplot(2,3,2+i*3)
    pl.scatter(x0, mean0)
    pl.xlabel(xlabel)
    pl.title(f'Array mean, axis=0; std={std0:n}')
    
    pl.subplot(2,3,3+i*3)
    pl.scatter(x1, mean1)
    pl.xlabel(xlabel)
    pl.title(f'Array mean, axis=1; std={std1:n}')

sc.figlayout()
pl.show()
```
So with a more realistic sim size -- 100,000 agents, say -- the correlations would only be about 0.2%. That might still matter if it were affecting the mean, but I'm not sure that I can think of a case where a small correlation between pairs of agents would matter. (Also: I confess I still don't understand the mathematical reason for where these correlations are coming from; I don't see a reason why, say, `(0.1+r)%1`, `(0.5+r)%1`, or `(0.9+r)%1` should be distributed diffrently from each other.)

### 2. Would the node-based approach solve this?
The modulus method is drawing the same number of bits' worth of random numbers as the node-based approach -- 2N rather than N². So I believe that the node-based approach would have the same induced correlations; it's just that they're more hidden since the per-edge probability is never explicitly calculated.

To be more specific, going back to the node-based transmission code (commit 46addb1), the probability per node is calculated as `p_acq_node = 1-np.prod(1-probs, axis=1)`, and then a single random number is drawn against this (`cases = r[people.slot[uids]] < p_acq_node`). So this will induce correlations between pairs of agents, just as with the modulus-based approach -- in fact the correlations should be even greater, since we're only using N random numbers to calculate transmission events, instead of 2N. (The other N random numbers are used to allocate a source agent if a transmission occurs.)

### 3. Can we solve networks?
CRN is based on the assumption that each agent carries around their own personal RNGs, which they use to collapse probabilities into events. But this assumption fails for interactions, since there are O(N²) pairwise interactions, whereas we have only O(N) RNGs. This is why the only CRN-safe network currently is `EmbeddingNet`, which is deterministic.

Is this a problem? I'm not sure. Where does the source of randomness in an interaction between two agents come from, other than the two agents themselves? We intentionally induce correlations all the time (e.g. using `np.random.choice()` to guarantee a node degree). So while assuming that the stochasticity between pairs of agents is due to 2N rather than N² independent sources is a _different_ assumption, it's not clear to me that it maps less well onto the real world. In fact: could use the modulus approach to make the networks CRN safe??

### Appendix: quantifying the difference
It's a bit hard to find an at-scale statistic that shows the issue; the clearest illustration I've come up with is that if you compute a histogram of the differences between adjacent points, it's very clearly different:
![image](https://github.com/starsimhub/starsim/assets/3239256/2ef1d512-2094-4d01-a574-ca8b7151a32d)
```py
import numpy as np
import sciris as sc
import pylab as pl

nreps = 3
n = 100
density = 1.0
sc.options(dpi=200)
pl.figure(figsize=(3*nreps, 10))
modulo = 'Modulus'
truly = 'Truly random'
kw = dict(alpha=0.7)

def diff_hist(arr, thresh=1e-3):
    a1diff = np.abs(np.diff(arr, axis=0))
    a2diff = np.abs(np.diff(arr, axis=1))
    vals = sc.cat(a1diff.flatten(), a2diff.flatten())
    vals = sc.sanitize(vals, replacenans=False)
    vals = vals[vals>thresh]
    hist, edges = np.histogram(vals, bins=n)
    return edges[:-1], hist

for rep in range(nreps):
    r1 = np.random.rand(n,n)
    v2a = np.random.rand(n)
    v2b = np.random.rand(n).reshape((n,1))
    r2a = np.tile(v2a, (n,1))
    r2b = np.tile(v2b, (1,n))
    r2 = (r2a + r2b) % 1.0
    
    mask = np.random.rand(n,n) > density
    r1[mask] = np.nan
    r2[mask] = np.nan
    
    pl.subplot(4,nreps,rep+nreps*0+1)
    pl.plot(v2a, 'o-', **kw)
    pl.plot(v2b, 'o-', **kw)
    pl.title('Input arrays for modulus')
    
    pl.subplot(4,nreps,rep+nreps*1+1)
    pl.imshow(r1)
    pl.title(truly)
    
    pl.subplot(4,nreps,rep+nreps*2+1)
    pl.imshow(r2)
    pl.title(modulo)
    
    pl.subplot(4,nreps,rep+nreps*3+1)
    pl.plot(*diff_hist(r2), 'o-', label=modulo, **kw)
    pl.plot(*diff_hist(r1), 'o-', label=truly, **kw)
    pl.legend()
    pl.title('Histograms of differences')
    
sc.figlayout()
pl.show()
```
The truly-random case is a straight line, while the modulus case has large deviations above and below the line. But how best to quantify that?
# Yet another investigation

I tried to capture the three methods a different way: in terms of repeated draws, and the probability after all edges are summed (which more accurately captures the transmission process).

## Code

Run with `python -i rnd_test.py` for best results; you might want to change/comment out the `sc.options(dpi=200)` line.
```py
import numpy as np
import sciris as sc
import pylab as pl

n = 50
seed = 1
reps = 1001
pause = 0.5
nframes = 20
nbins = n

np.random.seed(seed)
all_diffs = sc.objdict(Random=[], Modulo=[], Roulette=[])
row_diffs = sc.dcp(all_diffs)

# Choose logarithmically spaced frames to plot
frames = np.linspace(0, np.log(reps), nframes)
frames = np.exp(frames).astype(int) - 1

sc.options(dpi=200)
fig = pl.figure(figsize=(12,10))
sc.figlayout()


for rep in range(reps):
    
    # Truly random
    A = np.random.rand(n,n)
    
    # Modulo
    b1 = np.random.rand(n,1)
    b2 = np.random.rand(1,n)
    B1 = np.tile(b1, (1,n))
    B2 = np.tile(b2, (n,1))
    B = (B1 + B2) % 1.0
    
    # Roulette
    c1 = np.random.rand(n,1)
    C1 = np.tile(c1, (1,n))
    C2 = np.random.rand(n,n)
    C2 = C2*n/C2.sum(axis=0)
    C = C1*C2
    
    arrdict = sc.objdict(Random=A, Modulo=B, Roulette=C)
    
    for i,key,arr in arrdict.enumitems():
    
        v1 = np.diff(arr, axis=0)
        v2 = np.diff(arr, axis=1)
        v = sc.cat(v1.flatten(), v2.flatten()).tolist()
        all_diffs[key] += v
        
        rd = np.diff(arr.sum(axis=1)).tolist()
        row_diffs[key] += rd
    
    if rep in frames:
        fig.clear()
        for i,key,arr in arrdict.enumitems():
            
            pl.subplot(3,3,1+i*3)
            pl.imshow(arr, cmap='turbo')
            t = f'{key}: rep={rep}\nmean={arr.mean():0.3f}, std={arr.std():0.3f}\nmin={arr.min():0.3f}, max={arr.max():0.3f}'
            pl.title(t)
            
            pl.subplot(3,3,2+i*3)
            hist, bins = np.histogram(all_diffs[key], bins=nbins)
            w = np.diff(bins).mean()
            pl.bar(bins[:-1], hist, width=w, facecolor='k', lw=0)
            pl.title(f'{key}: all diffs rep={rep}')
            sc.commaticks()
            
            pl.subplot(3,3,3+i*3)
            rhist, rbins = np.histogram(row_diffs[key], bins=nbins)
            w = np.diff(rbins).mean()
            pl.bar(rbins[:-1], rhist, width=w, facecolor='k', lw=0)
            pl.title(f'{key}: row diffs rep={rep}')
            sc.commaticks()
        
        pl.pause(pause)
```

## Results
(On the final iteration)
![image](https://github.com/starsimhub/starsim/assets/3239256/7604f930-d03f-4fb3-85bb-07402709129c)

## Summary
I believe this shows that in the limit, all three approaches have different statistical distributions. However, all three are similar in terms of first-order statistics. The modulo method has higher-order statistics closer to the truly random method compared to the roulette method.

Alight, I can confirm that the modulo-based approach to transmission does indeed induce unwanted correlations. Code in the [crn_paper](https://github.com/starsimhub/crn_paper/blob/main/explore_tx_correlation.py) and copied below:
```python
import os
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from pandas.util import hash_pandas_object
import hashlib
import sciris as sc
from pathlib import Path

# PARAMETERS
n = 8 # Number of nodes
n_sources = 2 # Number of sources (seed infections)
reps = 10_000_000
edge_prob = 0.7 # Edge probability
trans_prob = 0.4
seed = 1

lbl = f'n{n}_s{n_sources}_r{reps}_ep{edge_prob}_tp{trans_prob}_seed{seed}'
figdir = os.path.join('figs', 'TransExploration', lbl)
Path(figdir).mkdir(parents=True, exist_ok=True)

sc.options(dpi=200)

def hash(df):
    #return int(hashlib.sha256(hash_pandas_object(df, index=True).values).hexdigest(), 16)
    return hashlib.sha256(hash_pandas_object(df, index=True).values).hexdigest()[:6]

def random(G):
    src = nx.get_node_attributes(G, 'infected')
    infected = []
    tx = []

    # Randomize in case a node is infected multiple times
    el = list(G.edges())
    np.random.shuffle(el)

    for (n1,n2) in el:
        # n1 --> n2
        if src[n1] and (not src[n2]) and (n2 not in infected) and (np.random.rand() < trans_prob):
            tx.append((n1, n2))
            infected.append(n2)
        # n2 --> n1
        if src[n2] and (not src[n1]) and (n1 not in infected) and (np.random.rand() < trans_prob):
            tx.append((n2, n1))
            infected.append(n1)
    return tx


def modulo(G):
    src = nx.get_node_attributes(G, 'infected')
    infected = []
    tx = []
    r1 = np.random.rand(n)
    r2 = np.random.rand(n)
    el = list(G.edges())
    np.random.shuffle(el)
    for (n1,n2) in el:
        # n1 --> n2
        if (src[n1]) and (not src[n2]) and (n2 not in infected) and (((r1[n1]+r2[n2])%1) < trans_prob):
            tx.append((n1, n2))
            infected.append(n2)
        # n2 --> n1
        if (src[n2]) and (not src[n1]) and (n1 not in infected) and (((r1[n2]+r2[n1])%1) < trans_prob):
            tx.append((n2, n1))
            infected.append(n1)
    return tx


def roulette(G):
    src = nx.get_node_attributes(G, 'infected')
    tx = []
    for n2 in G.nodes():
        # All nodes --> n2
        if not src[n2]:
            srcs = [n for n in G.neighbors(n2) if src[n]] #[src[n] for n in G.neighbors(n2)] # Count infected neighbors
            cnt = len(srcs)
            if np.random.rand() < 1-(1-trans_prob)**cnt:
                n1 = np.random.choice(srcs)
                tx.append((n1, n2))
    return tx
 

def transmit(G, trans_fn):
    # Align transmissions from tx_in if provided
    txs = {}
    counts = {}
 
    for _ in np.arange(reps):
        txl = trans_fn(G)
        tx = pd.DataFrame(txl, columns=['src', 'dst']).sort_values(['src', 'dst']).reset_index(drop=True)
        
        h = hash(tx)

        if h not in txs:
            txs[h] = tx

        counts[h] = counts.get(h, 0) + 1

    df = pd.DataFrame(counts.values(), index=pd.Index(counts.keys(), name='Hash'), columns=['Counts'])
    return txs, df


# Build the graph
G = nx.random_graphs.erdos_renyi_graph(n=n, p=edge_prob, seed=seed)

# Seed infections
infected = {i:False for i in range(n)}
sources = np.random.choice(a=range(n), size=n_sources, replace=False)
for source in sources:
    infected[source] = True
nx.set_node_attributes(G, infected, 'infected')
 
# Do transmissions via each method in parallel
results = sc.parallelize(transmit, iterkwargs=[{'trans_fn':random}, {'trans_fn':modulo}, {'trans_fn':roulette}], kwargs={'G':G}, die=True, serial=False)
tx, cnt = zip(*results)

df = pd.concat(cnt, axis=1) \
    .fillna(0) \
    .astype(int)
df.columns = ['Random', 'Modulo', 'Roulette']

# Manipulate results
df.reset_index(inplace=True)
dfm = df.melt(id_vars='Hash', var_name='Method', value_name='Count')

# Plot
g = sns.barplot(data=dfm, x='Hash', y='Count', hue='Method')
plt.xticks(rotation=90)
g.figure.tight_layout()
sc.savefig(os.path.join(figdir, 'Bar.png'), g.figure)

txc = tx[0].copy()
for i in range(1, len(tx)):
    txc.update(tx[i])
for h, v in txc.items():
    print(f'\nUnique transmission tree #{h}')
    print(v)
sc.savejson(os.path.join(figdir, 'Transmissions.json'), txc)

f = plt.figure(figsize=(12,8))
pos = nx.spring_layout(G, seed=3113794652)  # positions for all nodes
ec = nx.draw_networkx_edges(G, pos, alpha=0.2)
colors = ['red' if infected else 'blue' for infected in nx.get_node_attributes(G, 'infected').values()]
nc = nx.draw_networkx_nodes(G, pos, nodelist=G.nodes(), node_color=colors, node_size=100)
nx.draw_networkx_edges(G, pos, width=1.0, alpha=0.5)
nx.draw_networkx_labels(G, pos)
sc.savefig(os.path.join(figdir, 'Network.png'), f)

plt.show()
```

This code produces the following figures:
![Network](https://github.com/starsimhub/starsim/assets/15837372/b3782f2c-051e-47b2-ba42-bbe199d17d15)
Fig 1 shows the transmission network consisting of 8 nodes, two of which are initially infected

![BarHighRes](https://github.com/starsimhub/starsim/assets/15837372/94dbba91-7c1b-4d8c-976c-a76b9becd8f3)
Fig 2 shows the frequency of transmission trees for each of three transmission algorithms: Random, Modulo, and Roulette. The x-axis is a hash representing each unique transmission pattern, e.g. (1-->4, 5-->0, 5-->7). The transmissions associated with each hash are available in [Transmissions.json](https://github.com/user-attachments/files/15519378/Transmissions.json).

What we observe in the bar frequencies is that Random (blue) and Roulette (green) have very similar frequencies across all possible transmission patterns, but Modulo (orange) has a different frequency pattern. This difference indicates a statistical bias in transmission.

These results come from repeating transmission on the above network 10M times, so it is unlikely that these differences are caused by stochastic fluctuations. (We could perform a statistical test, but no need.)

Note that in my testing, there are some simple configurations in which the modulo approach matches the other two. For example, if testing with just 5 nodes, the output matches the other transmission algorithms. But even with `n=8` nodes as shown here, the bias is already clear.

Unless I have made a mistake in the implementation, I fear this study suggests that we need to revert to the "roulette-based" approach to transmission.

cc @cliffckerr 
TLDR: you're right, something is going on here!

I adapted your code to plot in a different way; definitely some performance optimizations to be had here (e.g. parallelize by reps instead of by function), but haven't implemented those for now (but did do some Sciris simplifications that I couldn't resist :D):

### Results
![image](https://github.com/starsimhub/starsim/assets/3239256/09532e74-93be-4c8a-9acd-c74daee4f21c)

### Code
```py
import numpy as np
import pandas as pd
import networkx as nx
import matplotlib.pyplot as plt
from pandas.util import hash_pandas_object
import sciris as sc

# PARAMETERS
n = 8 # Number of nodes
n_sources = 2 # Number of sources (seed infections)
reps = 100_000
edge_prob = 0.7 # Edge probability
trans_prob = 0.4
seed = 3

lbl = f'n{n}_s{n_sources}_r{reps}_ep{edge_prob}_tp{trans_prob}_seed{seed}'
figdir = sc.path('figs', 'TransExploration', lbl)
figdir.mkdir(parents=True, exist_ok=True)

sc.options(dpi=200)

def hash(df):
    #return int(hashlib.sha256(hash_pandas_object(df, index=True).values).hexdigest(), 16)
    return sc.sha(hash_pandas_object(df, index=True).values).hexdigest()[:6]

def random(G):
    src = nx.get_node_attributes(G, 'infected')
    infected = []
    tx = []

    # Randomize in case a node is infected multiple times
    el = list(G.edges())
    np.random.shuffle(el)

    for (n1,n2) in el:
        # n1 --> n2
        if src[n1] and (not src[n2]) and (n2 not in infected) and (np.random.rand() < trans_prob):
            tx.append((n1, n2))
            infected.append(n2)
        # n2 --> n1
        if src[n2] and (not src[n1]) and (n1 not in infected) and (np.random.rand() < trans_prob):
            tx.append((n2, n1))
            infected.append(n1)
    return tx


def modulo(G):
    src = nx.get_node_attributes(G, 'infected')
    infected = []
    tx = []
    r1 = np.random.rand(n)
    r2 = np.random.rand(n)
    el = list(G.edges())
    np.random.shuffle(el)
    for (n1,n2) in el:
        # n1 --> n2
        if (src[n1]) and (not src[n2]) and (n2 not in infected) and (((r1[n1]+r2[n2])%1) < trans_prob):
            tx.append((n1, n2))
            infected.append(n2)
        # n2 --> n1
        if (src[n2]) and (not src[n1]) and (n1 not in infected) and (((r1[n2]+r2[n1])%1) < trans_prob):
            tx.append((n2, n1))
            infected.append(n1)
    return tx


def roulette(G):
    src = nx.get_node_attributes(G, 'infected')
    tx = []
    for n2 in G.nodes():
        # All nodes --> n2
        if not src[n2]:
            srcs = [n for n in G.neighbors(n2) if src[n]] #[src[n] for n in G.neighbors(n2)] # Count infected neighbors
            cnt = len(srcs)
            if np.random.rand() < 1-(1-trans_prob)**cnt:
                n1 = np.random.choice(srcs)
                tx.append((n1, n2))
    return tx
 

def transmit(G, trans_fn):
    # Align transmissions from tx_in if provided
    txs = {}
    counts = {}
 
    for _ in np.arange(reps):
        txl = trans_fn(G)
        tx = pd.DataFrame(txl, columns=['src', 'dst']).sort_values(['src', 'dst']).reset_index(drop=True)
        
        h = hash(tx)

        if h not in txs:
            txs[h] = tx

        counts[h] = counts.get(h, 0) + 1

    df = pd.DataFrame(counts.values(), index=pd.Index(counts.keys(), name='Hash'), columns=['Counts'])
    return txs, df

T = sc.timer()

# Build the graph
G = nx.random_graphs.erdos_renyi_graph(n=n, p=edge_prob, seed=seed)

# Seed infections
infected = {i:False for i in range(n)}
sources = np.random.choice(a=range(n), size=n_sources, replace=False)
for source in sources:
    infected[source] = True
nx.set_node_attributes(G, infected, 'infected')
 
# Do transmissions via each method in parallel
results = sc.parallelize(transmit, iterkwargs=[{'trans_fn':random}, {'trans_fn':modulo}, {'trans_fn':roulette}], kwargs={'G':G}, die=True, serial=False)
tx, cnt = zip(*results)

df = pd.concat(cnt, axis=1) \
    .fillna(0) \
    .astype(int)
df.columns = ['Random', 'Modulo', 'Roulette']

# Manipulate results
df.reset_index(inplace=True)
dfm = df.melt(id_vars='Hash', var_name='Method', value_name='Count')

# Alternate plotting
keys = ['Random', 'Modulo', 'Roulette']
colors = {'Random':'k', 'Modulo':'slategray', 'Roulette':'red'}
z = sc.objdict()
for k in keys:
    z[k] = dfm[dfm.Method==k].Count.values

order = np.argsort(z.Random)

fig = plt.figure()
for k in keys:
    plt.plot(z[k][order], 'o-', alpha=0.7, label=k, c=colors[k])
plt.legend()

f = plt.figure(figsize=(12,8))
pos = nx.spring_layout(G, seed=3113794652)  # positions for all nodes
ec = nx.draw_networkx_edges(G, pos, alpha=0.2)
colors = ['red' if infected else 'blue' for infected in nx.get_node_attributes(G, 'infected').values()]
nc = nx.draw_networkx_nodes(G, pos, nodelist=G.nodes(), node_color=colors, node_size=100)
nx.draw_networkx_edges(G, pos, width=1.0, alpha=0.5)
nx.draw_networkx_labels(G, pos)
sc.savefig(figdir/'Network.png', f)

T.toc()
plt.show()
```
I still wasn't sure whether these differences would be enough to affect the overall results, but it seems like they are, unfortunately. With a (super)simple ABM approach, the modulo method results in fewer infections, which would be as expected if there are induced correlations between different possible transmission events:

### Results
![image](https://github.com/starsimhub/starsim/assets/3239256/a9530748-ecbb-4841-ab0c-ffcbd3ca330a)

### Code
```py
import numpy as np
import sciris as sc
import pylab as pl
# import numba as nb

# @nb.njit # No faster -- not sure why not
def get_probs(n, inds, trg, trans_prob):
    probs = np.ones(n)
    trg_inds = trg[inds]
    for trg_ind in trg_inds:
        probs[trg_ind] *= (1-trans_prob)
    probs = 1 - probs
    return probs

# PARAMETERS
p = sc.objdict(
    n = 100, # Number of nodes
    n_sources = 2, # Number of sources (seed infections)
    edge_prob = 0.1, # Edge probability
    trans_prob = 0.05,
    seed = 6,
    npts = 10,
)


class Sim:
    def __init__(self, p, method='random', repseed=0):
        self.p = p
        self.repseed = repseed
        self.method = method
        self.initialize()
        self.make_network()
    
    def initialize(self):
        np.random.seed(self.p.seed+self.repseed)
        n = self.p.n
        self.uids = np.arange(n)
        self.inf = np.zeros(n, dtype=int)
        self.inf[0:self.p.n_sources] = 1
        self.xpos = np.random.rand(n)
        self.ypos = np.random.rand(n)
        self.t = 0
        self.tvec = np.arange(p.npts+1)
        self.n_infected = p.n_sources*np.ones(p.npts+1)

    def make_network(self):
        p = self.p
        n_possible = p.n**2
        n_edges = np.random.binomial(n_possible, p.edge_prob)
        p1 = np.random.choice(self.uids, size=n_edges)
        p2 = np.random.choice(self.uids, size=n_edges)
        self.src = sc.cat(p1, p2) # Bidirectional transmission
        self.trg = sc.cat(p2, p1)
        
    def get_edges(self):
        possible = (self.inf[self.src]==1) * (self.inf[self.trg]==0)
        inds = sc.findinds(possible)
        return inds
    
    def trans_rand(self, inds):
        r = np.random.rand(len(inds))
        trans = self.p.trans_prob > r
        trans_inds = inds[trans]
        inf_inds = self.trg[trans_inds]
        return inf_inds
    
    def trans_modulo(self, inds):
        r1 = np.random.rand(self.p.n)
        r2 = np.random.rand(self.p.n)
        # np.random.shuffle(r2)
        r1 = r1[self.src[inds]]
        r2 = r2[self.trg[inds]]
        r = (r1 + r2) % 1.0
        trans = self.p.trans_prob > r
        trans_inds = inds[trans]
        inf_inds = self.trg[trans_inds]
        return inf_inds
    
    def trans_roulette(self, inds):
        trans_prob = self.p.trans_prob
        n = self.p.n
        trg = self.trg
        probs = get_probs(n, inds, trg, trans_prob)
        r = np.random.rand(self.p.n)
        inf_inds = probs > r
        return inf_inds
        
    def step(self):
        self.t += 1
        inds = self.get_edges()
        if self.method == 'random':
            inf_inds = self.trans_rand(inds)
        elif self.method == 'roulette':
            inf_inds = self.trans_roulette(inds)
        elif self.method == 'modulo':
            inf_inds = self.trans_modulo(inds)
        else:
            raise Exception('Method not found')
        self.inf[inf_inds] = True
        self.n_infected[self.t] = self.inf.sum()
            
    def run(self):
        for t in range(self.p.npts):
            self.step()
            if self.n_infected[self.t] == self.p.n:
                self.n_infected[self.t+1:] = self.p.n
                break
        return
            
    def plot_network(self):
        fig = pl.figure()
        x = self.xpos
        y = self.ypos
        inf = self.inf == 1
        ninf = self.inf == 0
        pl.scatter(x[ninf], y[ninf], c='green')
        pl.scatter(x[inf], y[inf], c='red')
        for p1,p2 in zip(self.src, self.trg):
            pl.plot([x[p1], x[p2]], [y[p1], y[p2]], alpha=0.2, c='k', lw=0.1)
        return fig
    
    def plot(self):
        fig = pl.figure()
        pl.plot(self.tvec, self.n_infected)
        return fig


def run_sim(sim):
    sim.run()
    return sim


class MultiSim:
    def __init__(self, p, method='random', reps=10):
        self.p = p
        self.method = method
        self.reps = reps
        self.initialize()
    
    def initialize(self):
        self.sims = [Sim(p=self.p, method=self.method, repseed=rep) for rep in range(self.reps)]
    
    def run(self):
        self.sims = sc.parallelize(run_sim, self.sims)
        self.results = np.zeros((self.p.npts+1, self.reps))
        for i,sim in enumerate(self.sims):
            self.results[:,i] = sim.n_infected
        
    def plot_all(self):
        pl.figure()
        pl.plot(self.results, alpha=0.2)
    
    def plot_stats(self, ax=None):
        if ax is None:
            fig,ax = pl.subplots(1)
        pl.sca(ax)
        qs = [0, 0.0025, 0.01, 0.05, 0.1, 0.25]
        res = self.results
        mean = res.mean(axis=1)
        std = res.std(axis=1)
        x = np.arange(len(mean))
        for q in qs:
            low = np.quantile(res, q, axis=1)
            high = np.quantile(res, 1-q, axis=1)
            pl.fill_between(x, low, high, alpha=0.1, facecolor='red')
        pl.plot(mean, 'o-', c='k', lw=2)
        pl.xlabel('Timestep')
        pl.ylabel('Number infected')
        pl.title(f'Method: {self.method}; stats: μ={mean.mean():0.4f}, σ={std.mean():0.4f}')
        pl.xlim(left=0, right=p.npts)
        pl.ylim(bottom=0, top=p.n)
        
        
if __name__ == '__main__':
    sc.options(dpi=200)
    T = sc.timer()
    fig, axes = pl.subplots(1,3, figsize=(18,8))
    for m,method in enumerate(['random', 'roulette', 'modulo']):
        msim = MultiSim(p=p, method=method, reps=1_000)
        T.tic()
        msim.run()
        T.toc(method)
        msim.plot_stats(ax=axes[m])
    sc.figlayout()
    pl.show()
```

Next step: figure out where exactly these correlations are coming from, and if there's an easy way to prevent them. If not, next step is trying to find a version of the roulette method that has negligible performance penalty.
Thanks for the deep investigation, @cliffckerr!

Agree something is off with the modulo approach for transmission on complex networks (it seems fine on small/simple networks). My explorations to determine the source of the correlation led to pyramid and larger structures of connected nodes wherein non-random patterns result from reusing random draws along closed loops of 4+ nodes. But it wasn't a full investigation.

Regarding performance, my final attempt at the "roulette" (that's not the right name for this approach) implementation seemed to have performance on par with random transmission. I don't recall seeing much of a performance penalty, as with the earlier implementations. Again, wasn't a full analysis, but it didn't seem too bad. My understanding is that the desire to replace it with the modulo approach was more about code simplicity, but it has been long enough now that I don't recall if there was also a performance difference.

I have started reimplementing the "acquisition-based" approach (better name?), should have time this quite week to finish and then we can look at performance in PR.
Now I'm confused again. I would've expected the matrix to have the right mean, but the row-sum counts to be off due to induced correlations. But they match!

These results compute the NxN transmission matrix for the truly-random and modulo case. Rows are targets, columns are sources. A target agent is "infected" if at least one edge (entry in that row) is above the threshold (i.e. if the sum of the row is greater than 0).

Based on the investigations above, I had expected that the statistics would match for the full matrix in terms of the number of above-threshold edges, but that the number of infected agents _wouldn't_ match due to induced correlations (N^2 vs. 2N). But in fact they both match:
![image](https://github.com/starsimhub/starsim/assets/3239256/58bed0df-3b3b-4de3-96a9-4c862ceaf3f4)

This makes me very confused, because now I'm not sure what was causing the difference in the "ABM" results above. Conceptually I thought this case and that case were equivalent.

```py
"""
What -- this shows no difference in the number of infections!!
"""

import numpy as np
import pylab as pl
import sciris as sc

n = 100
reps = 10_000
nc = 4
thresh = 1-(0.5/n)
sc.options(dpi=200)
pl.figure(figsize=(16,12))

d = sc.objdict(defaultdict=sc.autolist)

def matsort(arr):
    s0 = arr.sum(axis=0)
    s1 = arr.sum(axis=1)
    o0 = np.argsort(s0)
    o1 = np.argsort(s1)
    arr = arr[o1,:][:,o0]
    return arr

for rep in range(reps):

    # Truly random
    A = np.random.rand(n,n)
    A = matsort(A)
    
    # Modulo
    b1 = np.random.rand(n,1)
    b2 = np.random.rand(1,n)
    B1 = np.tile(b1, (1,n))
    B2 = np.tile(b2, (n,1))
    B = (B1 + B2) % 1.0
    B = matsort(B)
    
    # Thresholded
    At = A>thresh
    Bt = B>thresh
    v = sc.objdict()
    v.Asum = A.sum()
    v.Bsum = B.sum()
    v.Atsum = (At.sum(axis=1) > 0).sum()
    v.Btsum = (Bt.sum(axis=1) > 0).sum()
    
    for k,val in v.items():
        d[k] += val
    
    
    
#%% Plotting
for i,k,val in d.enumitems():
    pl.subplot(2,2,i+1)
    pl.hist(val, bins=int(np.sqrt(reps)))
    pl.title(f'k {np.array(val).mean():n}')

sc.figlayout()
pl.show()
```
And indeed, there is no measurable impact on Starsim results between v0.5.2 (modulo) and v0.5.3 (roulette).

This script stores 10,000 simulations from each branch (modulo vs. roulette) and compares the outputs (specifically `sim.summary`). All the results match to within ±2 SEM (even without doing a correction for multiple comparisons):
![image](https://github.com/starsimhub/starsim/assets/3239256/563a924a-a231-4447-99b5-ed5f058ad5ed)

Code:
```py
"""
To use:
1. Run this script on `main` to generate `rng_stat_modulo.json`
2. Change to `revert-transmission` and run the script again to generate `rng_stat_roulette.json`
3. Run a third time to plot the results

This script should take about 2 minutes to run on `main` and 10 minutes to run on `revert-transmission`.
"""
import os
import numpy as np
import sciris as sc
import starsim as ss
import pylab as pl

reps = 10_000
keys = {'0.5.2':'modulo', '0.5.3':'roulette'}
filenames = {k:f'rng_stat_{k}.json' for k in keys.values()}


def make_sim(seed):
    pars = dict(
        rand_seed = seed,
        n_agents = 1000,
        end = 2020,
        diseases = ['sir', 'sis'],
        networks = 'random',
        verbose = 0,
    )
    sim = ss.Sim(pars)
    return sim

def run_sim_summary(seed):
    sim = make_sim(seed)
    sim.run()
    return sim.summary

def plot_sim():
    sim = make_sim()
    fig = sim.run().plot()
    return fig

def run_sims():
    ver = ss.__version__
    sc.heading(f'Running {ver}...')
    T = sc.timer()
    reslist = sc.parallelize(run_sim_summary, np.arange(reps))
    res = dict()
    for k in reslist[0].keys():
        res[k] = []
    for entry in reslist:
        for k in res.keys():
            res[k].append(entry[k])
    key = keys[ver]
    filename = filenames[key]
    sc.savejson(filename, res)
    T.toc()
    return

def get_raw():
    raw = sc.objdict()
    for k,filename in filenames.items():
        if not os.path.exists(filename):
            run_sims()
            print('Now change branch (if needed) and re-run')
            break
        raw[k] = sc.loadjson(filename)
    return raw

def get_res(raw):
    res = sc.objdict()
    for k,v in raw.items():
        res[k] = sc.objdict()
        for k2,v2 in v.items():
            mean = np.mean(v2)
            sem = np.std(v2)/np.sqrt(len(v2))
            res[k][k2] = sc.objdict(mean=mean, sem=sem)
    return res

def plot(raw, res):
    keys = res[0].keys()
    nkeys = len(keys)
    
    sc.options(dpi=150)
    fig = pl.figure(figsize=(16,16))
    count = 0
    for i,k in enumerate(keys):
        mod = res.modulo[k]
        rou = res.roulette[k]
        z = 2
        modmin = mod.mean - z*mod.sem
        modmax = mod.mean + z*mod.sem
        roumin = rou.mean - z*rou.sem
        roumax = rou.mean + z*rou.sem
        statsig = True if (modmax < roumin) or (modmin > roumax) else False
        for i2,k2 in enumerate(res.keys()):
            count += 1
            this = res[k2][k]
            title = f'{k2}\n{k}\n{this.mean:n}±{2*this.sem:n}'
            if statsig:
                title += '\n(mismatch!)'
            pl.subplot(nkeys//2,4,count)
            pl.hist(raw[k2][k], bins=int(np.sqrt(reps)))
            pl.title(title)
        
    sc.figlayout()
    pl.show()
    return fig
    
    
if __name__ == '__main__':
    raw = get_raw()
    res = get_res(raw)
    fig = plot(raw, res)
```

I feel we're still not quite at the bottom of this, but I'm feeling reassured that there doesn't seem to be any difference in the final results, even though I agree that some of the intermediate calculation results differ.
I suspect the correlations of the modulo approach are lost in the random mixing network in the above examples.
I had a quick look at this, my impression is that the issue is that the output of `(a+b)%1` is continuous in the sense that a small change to `a` means there is a small change in `(a+b)%1`. So if you take Cliff's matrix, and sort it by the values of the first random, you get something like

![image](https://github.com/starsimhub/starsim/assets/755796/07033977-4657-4494-b246-808ad2c149d7)
 
If that's the case, then we might have better results with a function that behaves more like a hash, where a small difference in `a` or `b` results in large change to the output. Just by trial and error, something like this might be suitable:

```python
def f(a,b):
    a1 = (a*1e9).astype(int)
    b1 = (b*1e9).astype(int)
    c = np.bitwise_xor(a1*b1, a1-b1)
    c = c.astype(float)
    c = c-c.min()
    c = c/c.max()
    return c
```

For Cliff's test this gives

![image](https://github.com/starsimhub/starsim/assets/755796/05899d7d-f228-48ac-888a-d4fc48775d8b)

The sorted matrix looks like

![image](https://github.com/starsimhub/starsim/assets/755796/b1f6bc80-57d4-4b59-959b-d1298f7f506d)


For Dan's original test

```
indep 0.49965 0.49917657968789975
modulo 0.500639 0.6689087176373789
xor 0.500523 0.5012711661150285
```

For Cliff's network

![image](https://github.com/starsimhub/starsim/assets/755796/ed90c852-688f-416a-b365-0ef1a00982b7)

The distribution does seem pretty uniform

```
    a = np.random.rand(1000000)
    b = np.random.rand(1000000)
    a1 = (a*1e9).astype(int)
    b1 = (b*1e9).astype(int)
    c = np.bitwise_xor(a1*b1, a1-b1)
    c = c.astype(float)
    c = c-c.min()
    c = c/c.max()
    x,y = np.histogram(a, 1000)
    pl.plot(y[1:],x, label='a', alpha=0.5)
    x,y = np.histogram(b, 1000)
    pl.plot(y[1:],x, label='b', alpha=0.5)
    x,y = np.histogram(c, 1000)
    pl.plot(y[1:],x, label='c', alpha=0.5)
    pl.legend()
```

![image](https://github.com/starsimhub/starsim/assets/755796/a442ccc9-3942-452d-9553-f2043a88a27c)

And performance wise it's actually a little faster than the original modulus code. I haven't had time to try all of the other tests above and admittedly the selection of that function is a bit arbitrary - but this might be a promising avenue to explore



@RomeshA , this is amazing.
I'm curious how much performance difference there really is in practice. On TBsim, the difference between modulo and acquisition-based (roulette) is negligible. But benchmark.json shows an order of magnitude difference. Did you get to the bottom of this difference with any of your follow-up testing?


Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 58a3458f..77c86f8f 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -7,6 +7,12 @@ What's new
 All notable changes to the codebase are documented in this file. Changes that may result in differences in model output, or are required in order to run an old parameter set with the current version, are flagged with the term "Regression information".
 
 
+Version 0.5.7 (2024-06-27)
+--------------------------
+- Implemented a new ``ss.combine_rands()`` function based on a bitwise-XOR, since the previous modulo-based approach could introduce correlations between pairs of agents.
+- *GitHub info*: PR `546 <https://github.com/starsimhub/starsim/pull/546>`
+
+
 Version 0.5.6 (2024-06-22)
 --------------------------
 - ``ss.Infection.make_new_cases()`` now returns the index of the network associated with each transmission event
@@ -15,6 +21,8 @@ Version 0.5.6 (2024-06-22)
 - ``ss.uids.to_numpy()`` now returns a view rather than a copy
 - ``ss.bernoulli.filter()`` now supports ``ss.BoolArr`` as an input, where the filtering will operate on the ``uids`` returned by ``ss.BoolArr.uids``
 - ``ss.uids()`` supports construction from ``set`` objects (via ``np.fromiter()``)
+- *GitHub info*: PR `565 <https://github.com/starsimhub/starsim/pull/555>`_
+
 
 Version 0.5.5 (2024-06-19)
 --------------------------
diff --git a/starsim/disease.py b/starsim/disease.py
index a3611522..027a5b37 100644
--- a/starsim/disease.py
+++ b/starsim/disease.py
@@ -152,8 +152,8 @@ def __init__(self, *args, **kwargs):
         )
 
         # Define random number generators for make_new_cases
-        self.rng_target = ss.random(name='target')
-        self.rng_source = ss.random(name='source')
+        self.rng_target = ss.randint(name='target')
+        self.rng_source = ss.randint(name='source')
         return
     
     def init_pre(self, sim):
@@ -277,12 +277,10 @@ def make_new_cases(self):
                 beta_per_dt = net.beta_per_dt(disease_beta=beta, dt=self.sim.dt)
                 p_transmit = rel_trans[src] * rel_sus[trg] * beta_per_dt
 
-                # Generate a new random number based on the two other random numbers -- 3x faster than `rvs = np.remainder(rvs_s + rvs_t, 1)`
+                # Generate a new random number based on the two other random numbers
                 rvs_s = self.rng_source.rvs(src)
                 rvs_t = self.rng_target.rvs(trg)
-                rvs = rvs_s + rvs_t
-                inds = np.where(rvs>1.0)[0]
-                rvs[inds] -= 1
+                rvs = ss.combine_rands(rvs_s, rvs_t)
                 
                 new_cases_bool = rvs < p_transmit
                 new_cases.append(trg[new_cases_bool])
diff --git a/starsim/diseases/hiv.py b/starsim/diseases/hiv.py
index db8fdc06..f97c3748 100644
--- a/starsim/diseases/hiv.py
+++ b/starsim/diseases/hiv.py
@@ -105,7 +105,7 @@ def apply(self, sim):
             return
 
         ti_delay = 1 # 1 time step delay TODO
-        recently_infected = ss.true((sim.people.hiv.ti_infected == sim.ti-ti_delay) & sim.people.alive)
+        recently_infected = (sim.people.hiv.ti_infected == sim.ti-ti_delay).uids
 
         n_added = 0
         if len(recently_infected) > 0:
@@ -115,7 +115,7 @@ def apply(self, sim):
             n_added = len(inds)
 
         # Add result
-        self.results['n_art'][sim.ti] = np.count_nonzero(sim.people.alive & sim.people.hiv.on_art)
+        self.results['n_art'][sim.ti] = np.count_nonzero(sim.people.hiv.on_art)
 
         return n_added
 
diff --git a/starsim/distributions.py b/starsim/distributions.py
index 3568fe58..07134865 100644
--- a/starsim/distributions.py
+++ b/starsim/distributions.py
@@ -473,7 +473,7 @@ def call_pars(self):
         for key,val in self._pars.items():
             
             # If the parameter is callable, then call it
-            if callable(val): 
+            if callable(val) and not isinstance(val, type): # Types can appear as callable
                 size_par = uids if uids is not None else size
                 out = val(self.module, self.sim, size_par) # TODO: swap order to sim, module, size?
                 val = np.asarray(out) # Necessary since UIDArrays don't allow slicing # TODO: check if this is correct
@@ -718,11 +718,21 @@ def sync_pars(self):
 
 
 class randint(Dist):
-    """ Random integers, values on the interval [low, high-1] (i.e. "high" is excluded) """
-    def __init__(self, low=0, high=2,  **kwargs):
-        super().__init__(distname='integers', dist=sps.randint, low=low, high=high, **kwargs)
+    """ Randint distribution, integer values on interval [low, high) using numpy.random.Generator.integers """
+    def __init__(self, low=0, high=np.iinfo(ss.dtypes.int).max, dtype=int, **kwargs):
+        if ss.options._centralized:
+            # randint because we're accessing via numpy.random
+            super().__init__(distname='randint', low=low, high=high, dtype=dtype, **kwargs)
+        else:
+            # integers instead of randint because interfacing a numpy.random.Generator
+            super().__init__(distname='integers', low=low, high=high, dtype=dtype, **kwargs)
         return
     
+    def ppf(self, rands):
+        p = self._pars
+        rvs = rands * (p.high + 1 - p.low) + p.low
+        rvs = rvs.astype(self.pars['dtype'])
+        return rvs
 
 class weibull(Dist):
     """ Weibull distribution -- NB, uses SciPy rather than NumPy """
diff --git a/starsim/utils.py b/starsim/utils.py
index 932cb5f8..71b5e75f 100644
--- a/starsim/utils.py
+++ b/starsim/utils.py
@@ -9,12 +9,10 @@
 import numba as nb
 import pandas as pd
 
-ss_float_ = ss.dtypes.float
-
 # %% Helper functions
 
 # What functions are externally visible -- note, this gets populated in each section below
-__all__ = ['ndict', 'warn', 'unique', 'find_contacts']
+__all__ = ['ndict', 'warn', 'unique', 'find_contacts', 'combine_rands']
 
 
 class ndict(sc.objdict):
@@ -289,9 +287,9 @@ def standardize_data(data=None, metadata=None, min_year=1800, out_of_range=0, de
 
     # Add defaults for year and age
     if 'year' not in index:
-        index['year'] = np.full(values.shape, fill_value=default_year, dtype=ss_float_)
+        index['year'] = np.full(values.shape, fill_value=default_year, dtype=ss.dtypes.float)
     if 'age' not in index:
-        index['age'] = np.full(values.shape, fill_value=default_age, dtype=ss_float_)
+        index['age'] = np.full(values.shape, fill_value=default_age, dtype=ss.dtypes.float)
 
     # Reorder the index so that it starts with age first
     index.insert(0, 'age', index.pop('age'))
@@ -317,4 +315,20 @@ def standardize_data(data=None, metadata=None, min_year=1800, out_of_range=0, de
     # Perform a final sort to prevent "indexing past lexsort depth" warnings
     output = output.sort_index()
 
-    return output
\ No newline at end of file
+    return output
+
+
+def combine_rands(a, b):
+    """
+    Efficient algorithm for combining two arrays of random numbers into one
+    
+    Args:
+        a (array): array of random integers between np.iinfo(np.int64).min and np.iinfo(np.int64).max
+        b (array): ditto, same size as a
+        
+    Returns:
+        A new array of random numbers the same size as a and b
+    """
+    c = np.bitwise_xor(a*b, a-b).astype(np.uint64)
+    u = c / np.iinfo(np.uint64).max
+    return u
diff --git a/starsim/version.py b/starsim/version.py
index 72ba654b..3132f8f6 100644
--- a/starsim/version.py
+++ b/starsim/version.py
@@ -4,6 +4,6 @@
 
 __all__ = ['__version__', '__versiondate__', '__license__']
 
-__version__ = '0.5.6'
-__versiondate__ = '2024-06-22'
+__version__ = '0.5.7'
+__versiondate__ = '2024-06-27'
 __license__ = f'Starsim {__version__} ({__versiondate__}) — © 2023-2024 by IDM'


Test patch:
diff --git a/tests/baseline.json b/tests/baseline.json
index a5760dcc..40513b4e 100644
--- a/tests/baseline.json
+++ b/tests/baseline.json
@@ -1,26 +1,26 @@
 {
   "summary": {
     "yearvec": 2010.1000000000026,
-    "births_new": 65.23529411764706,
-    "births_cumulative": 3261.401960784314,
-    "births_cbr": 29.71230613735741,
-    "deaths_new": 44.64705882352941,
-    "deaths_cumulative": 2225.294117647059,
-    "deaths_cmr": 20.330395920260354,
-    "sir_n_susceptible": 2437.1960784313724,
-    "sir_n_infected": 3560.6372549019607,
-    "sir_n_recovered": 4979.85294117647,
-    "sir_prevalence": 0.3333337529382368,
-    "sir_new_infections": 127.94117647058823,
-    "sir_cum_infections": 13050.0,
-    "sis_n_susceptible": 3938.6470588235293,
-    "sis_n_infected": 7039.0392156862745,
-    "sis_prevalence": 0.6287209006949334,
-    "sis_new_infections": 209.7450980392157,
-    "sis_cum_infections": 21394.0,
-    "sis_rel_sus": 0.4606516409026911,
-    "n_alive": 10977.686274509804,
-    "new_deaths": 45.705882352941174,
-    "cum_deaths": 4625.0
+    "births_new": 65.22549019607843,
+    "births_cumulative": 3261.1862745098038,
+    "births_cbr": 29.718240091091186,
+    "deaths_new": 44.61764705882353,
+    "deaths_cumulative": 2225.8529411764707,
+    "deaths_cmr": 20.325548477328535,
+    "sir_n_susceptible": 2412.774509803922,
+    "sir_n_infected": 3579.901960784314,
+    "sir_n_recovered": 4981.098039215686,
+    "sir_prevalence": 0.33508936820289537,
+    "sir_new_infections": 128.97058823529412,
+    "sir_cum_infections": 13155.0,
+    "sis_n_susceptible": 4001.872549019608,
+    "sis_n_infected": 6971.901960784314,
+    "sis_prevalence": 0.6229119856234305,
+    "sis_new_infections": 209.09803921568627,
+    "sis_cum_infections": 21328.0,
+    "sis_rel_sus": 0.4649922085611042,
+    "n_alive": 10973.774509803921,
+    "new_deaths": 45.76470588235294,
+    "cum_deaths": 4629.0
   }
 }
\ No newline at end of file
diff --git a/tests/benchmark.json b/tests/benchmark.json
index 9d7ab625..25d1fe6c 100644
--- a/tests/benchmark.json
+++ b/tests/benchmark.json
@@ -1,12 +1,12 @@
 {
   "time": {
-    "initialize": 0.015,
-    "run": 0.927
+    "initialize": 0.018,
+    "run": 1.264
   },
   "parameters": {
     "n_agents": 10000,
     "n_years": 20,
     "dt": 0.2
   },
-  "cpu_performance": 0.7750438800086998
+  "cpu_performance": 0.9766497596430002
 }
\ No newline at end of file
diff --git a/tests/devtests/devtest_rngtrans.py b/tests/devtests/devtest_rngtrans.py
new file mode 100644
index 00000000..2d3c6ac9
--- /dev/null
+++ b/tests/devtests/devtest_rngtrans.py
@@ -0,0 +1,166 @@
+import networkx as nx
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sns
+from pandas.util import hash_pandas_object
+import hashlib
+import sciris as sc
+
+from starsim.utils import combine_rands
+
+import warnings
+warnings.filterwarnings("ignore", "overflow encountered in scalar subtract")
+warnings.filterwarnings("ignore", "overflow encountered in scalar multiply")
+ 
+n = 8 # Number of nodes
+n_sources = 2 # Number of sources (seed infections)
+reps = 1_000_000
+edge_prob = 0.7 # Edge probability
+trans_prob = 0.4
+seed = 1
+
+def hash(df):
+    #return int(hashlib.sha256(hash_pandas_object(df, index=True).values).hexdigest(), 16)
+    return hashlib.sha256(hash_pandas_object(df, index=True).values).hexdigest()[:6]
+
+def random(G):
+    src = nx.get_node_attributes(G, 'infected')
+    infected = []
+    tx = []
+
+    # Randomize in case a node is infected multiple times
+    el = list(G.edges())
+    np.random.shuffle(el)
+
+    for (n1,n2) in el:
+        # n1 --> n2
+        if src[n1] and (not src[n2]) and (n2 not in infected) and (np.random.rand() < trans_prob):
+            tx.append((n1, n2))
+            infected.append(n2)
+        # n2 --> n1
+        if src[n2] and (not src[n1]) and (n1 not in infected) and (np.random.rand() < trans_prob):
+            tx.append((n2, n1))
+            infected.append(n1)
+    return tx
+
+
+def modulo(G):
+    src = nx.get_node_attributes(G, 'infected')
+    infected = []
+    tx = []
+    r1 = np.random.rand(n)
+    r2 = np.random.rand(n)
+    el = list(G.edges())
+    np.random.shuffle(el)
+    for (n1,n2) in el:
+        # n1 --> n2
+        if (src[n1]) and (not src[n2]) and (n2 not in infected) and (((r1[n1]+r2[n2])%1) < trans_prob):
+            tx.append((n1, n2))
+            infected.append(n2)
+        # n2 --> n1
+        if (src[n2]) and (not src[n1]) and (n1 not in infected) and (((r1[n2]+r2[n1])%1) < trans_prob):
+            tx.append((n2, n1))
+            infected.append(n1)
+    return tx
+
+
+def combine_hash(G):
+    src = nx.get_node_attributes(G, 'infected')
+    infected = []
+    tx = []
+    el = list(G.edges())
+    r1 = np.random.randint(low=0, high=np.iinfo(np.uint64).max, dtype=np.uint64, size=n)
+    r2 = np.random.randint(low=0, high=np.iinfo(np.uint64).max, dtype=np.uint64, size=n)
+    np.random.shuffle(el)
+    for (n1,n2) in el:
+        # n1 --> n2
+        r = combine_rands(r1[n1], r2[n2])
+        if (src[n1]) and (not src[n2]) and (n2 not in infected) and (r < trans_prob):
+            tx.append((n1, n2))
+            infected.append(n2)
+        # n2 --> n1
+        r = combine_rands(r1[n2], r2[n1])
+        if (src[n2]) and (not src[n1]) and (n1 not in infected) and (r < trans_prob):
+            tx.append((n2, n1))
+            infected.append(n1)
+    return tx
+
+def roulette(G):
+    src = nx.get_node_attributes(G, 'infected')
+    tx = []
+    for n2 in G.nodes():
+        # All nodes --> n2
+        if not src[n2]:
+            srcs = [n for n in G.neighbors(n2) if src[n]] #[src[n] for n in G.neighbors(n2)] # Count infected neighbors
+            cnt = len(srcs)
+            if np.random.rand() < 1-(1-trans_prob)**cnt:
+                n1 = np.random.choice(srcs)
+                tx.append((n1, n2))
+    return tx
+ 
+
+def transmit(G, trans_fn):
+    # Align transmissions from tx_in if provided
+    txs = {}
+    counts = {}
+ 
+    for _ in np.arange(reps):
+        txl = trans_fn(G)
+        tx = pd.DataFrame(txl, columns=['src', 'dst']).sort_values(['src', 'dst']).reset_index(drop=True)
+        
+        h = hash(tx)
+
+        if h not in txs:
+            txs[h] = tx
+
+        counts[h] = counts.get(h, 0) + 1
+
+    df = pd.DataFrame(counts.values(), index=pd.Index(counts.keys(), name='Hash'), columns=['Counts'])
+    return txs, df
+
+
+# Build the graph
+G = nx.random_graphs.erdos_renyi_graph(n=n, p=edge_prob, seed=seed)
+
+# Seed infections
+infected = {i:False for i in range(n)}
+sources = np.random.choice(a=range(n), size=n_sources, replace=False)
+for source in sources:
+    infected[source] = True
+nx.set_node_attributes(G, infected, 'infected')
+ 
+# Do transmissions via each method in parallel
+results = sc.parallelize(transmit, iterkwargs=[{'trans_fn':random}, {'trans_fn':modulo}, {'trans_fn':combine_hash}, {'trans_fn':roulette}], kwargs={'G':G}, die=True, serial=False)
+tx, cnt = zip(*results)
+
+df = pd.concat(cnt, axis=1) \
+    .fillna(0) \
+    .astype(int)
+df.columns = ['Random', 'Modulo', 'CombineHash', 'Roulette']
+
+# Manipulate results
+df.reset_index(inplace=True)
+dfm = df.melt(id_vars='Hash', var_name='Method', value_name='Count')
+
+# Plot
+g = sns.barplot(data=dfm, x='Hash', y='Count', hue='Method')
+plt.xticks(rotation=90)
+g.figure.tight_layout()
+
+txc = tx[0].copy()
+for i in range(1, len(tx)):
+    txc.update(tx[i])
+for h, v in txc.items():
+    print(f'\nUnique transmission tree #{h}')
+    print(v)
+
+plt.figure()
+pos = nx.spring_layout(G, seed=3113794652)  # positions for all nodes
+ec = nx.draw_networkx_edges(G, pos, alpha=0.2)
+colors = ['red' if infected else 'blue' for infected in nx.get_node_attributes(G, 'infected').values()]
+nc = nx.draw_networkx_nodes(G, pos, nodelist=G.nodes(), node_color=colors, node_size=100)
+nx.draw_networkx_edges(G, pos, width=1.0, alpha=0.5)
+nx.draw_networkx_labels(G, pos)
+
+plt.show()
\ No newline at end of file
diff --git a/tests/test_plugins.py b/tests/test_plugins.py
index fc94e5d4..0bc9030f 100644
--- a/tests/test_plugins.py
+++ b/tests/test_plugins.py
@@ -50,7 +50,7 @@ def update(self):
 
 class Penicillin(ss.Intervention):
     """ Create a penicillin (BPG) intervention for treating syphilis """
-    def __init__(self, year=2020, prob=0.5):
+    def __init__(self, year=2020, prob=0.8):
         super().__init__() # Initialize the intervention
         self.prob = prob # Store the probability of treatment
         self.year = year
@@ -76,7 +76,7 @@ def apply(self, sim):
 
 def make_args():
     """ Make people, HIV, syphilis, and network """
-    pars = dict(n_agents=1000, verbose=0)
+    pars = dict(n_agents=2000, verbose=0)
     mf = ss.MFNet(duration=ss.lognorm_ex(mean=5, stdev=0.5))
     hiv = ss.HIV(beta={'mf': [0.0008, 0.0004]}, init_prev=0.2)
     syph = ss.Syphilis(beta={'mf': [0.1, 0.05]}, init_prev=0.05)
diff --git a/tests/test_randomness.py b/tests/test_randomness.py
index 2b38f672..0f543112 100644
--- a/tests/test_randomness.py
+++ b/tests/test_randomness.py
@@ -6,6 +6,7 @@
 import numpy as np
 import sciris as sc
 import starsim as ss
+import scipy.stats as sps
 import matplotlib.pyplot as pl
 
 n = 5 # Default number of samples
@@ -285,6 +286,29 @@ def test_independence(do_plot=False, thresh=0.1):
     return sim
     
 
+def test_combine_rands(do_plot=False):
+    n = int(1e6)
+    atol = 1e-3
+    target = 0.5
+    np.random.seed(2)
+    a = np.random.randint(low=0, high=np.iinfo(np.uint64).max, dtype=np.uint64, size=n)
+    b = np.random.randint(low=0, high=np.iinfo(np.uint64).max, dtype=np.uint64, size=n)
+    c = ss.combine_rands(a, b)
+    if do_plot:
+        pl.figure()
+        for i,k,v in sc.objdict(a=a,b=b,combined=c).enumitems():
+            pl.subplot(3,1,i+1)
+            pl.hist(v)
+            pl.title(k)
+        sc.figlayout()
+        pl.show()
+    
+    mean = c.mean()
+    assert np.isclose(mean, target, atol=atol), f'Expected value to be 0.5±{atol}, not {mean}'
+    ks = sps.kstest(c, sps.uniform(0,1).cdf)
+    assert ks.pvalue > 0.05, f'Distribution does not seem to be uniform, p={ks.pvalue}<0.05'
+    return c
+
 
 # %% Run as a script
 if __name__ == '__main__':
@@ -297,6 +321,7 @@ def test_independence(do_plot=False, thresh=0.1):
     o4 = test_order(n)
     o5 = test_worlds(do_plot=do_plot)
     o6 = test_independence(do_plot=do_plot)
+    o7 = test_combine_rands(do_plot=do_plot)
 
     T.toc()
 
diff --git a/tests/test_syphilis.py b/tests/test_syphilis.py
index 731eaf91..815f50e6 100644
--- a/tests/test_syphilis.py
+++ b/tests/test_syphilis.py
@@ -4,8 +4,9 @@
 
 # %% Imports and settings
 import numpy as np
-import starsim as ss
 import pandas as pd
+import sciris as sc
+import starsim as ss
 import matplotlib.pyplot as plt
 
 quick_run = True
@@ -83,7 +84,7 @@ def apply(self, sim):
         return
 
 
-def test_syph(dt=1, n_agents=500):
+def test_syph(dt=1, n_agents=500, do_plot=False):
 
     sim_kwargs = make_syph_sim(dt=dt, n_agents=n_agents)
     sim = ss.Sim(analyzers=check_states(), **sim_kwargs)
@@ -92,32 +93,33 @@ def test_syph(dt=1, n_agents=500):
     # Check plots
     burnin = 0
     pi = int(burnin/sim.dt)
-
-    fig, ax = plt.subplots(2, 2)
-    ax = ax.ravel()
-    ax[0].stackplot(
-        sim.yearvec[pi:],
-        # sim.results.syphilis.n_susceptible[pi:],
-        sim.results.syphilis.n_congenital[pi:],
-        sim.results.syphilis.n_exposed[pi:],
-        sim.results.syphilis.n_primary[pi:],
-        sim.results.syphilis.n_secondary[pi:],
-        (sim.results.syphilis.n_latent_temp[pi:]+sim.results.syphilis.n_latent_long[pi:]),
-        sim.results.syphilis.n_tertiary[pi:],
-    )
-    ax[0].legend(['Congenital', 'Exposed', 'Primary', 'Secondary', 'Latent', 'Tertiary'], loc='lower right')
-
-    ax[1].plot(sim.yearvec[pi:], sim.results.syphilis.prevalence[pi:])
-    ax[1].set_title('Syphilis prevalence')
-
-    ax[2].plot(sim.yearvec[pi:], sim.results.n_alive[pi:])
-    ax[2].set_title('Population')
-
-    ax[3].plot(sim.yearvec[pi:], sim.results.syphilis.new_infections[pi:])
-    ax[3].set_title('New infections')
-
-    fig.tight_layout()
-    plt.show()
+    
+    if do_plot:
+        fig, ax = plt.subplots(2, 2)
+        ax = ax.ravel()
+        ax[0].stackplot(
+            sim.yearvec[pi:],
+            # sim.results.syphilis.n_susceptible[pi:],
+            sim.results.syphilis.n_congenital[pi:],
+            sim.results.syphilis.n_exposed[pi:],
+            sim.results.syphilis.n_primary[pi:],
+            sim.results.syphilis.n_secondary[pi:],
+            (sim.results.syphilis.n_latent_temp[pi:]+sim.results.syphilis.n_latent_long[pi:]),
+            sim.results.syphilis.n_tertiary[pi:],
+        )
+        ax[0].legend(['Congenital', 'Exposed', 'Primary', 'Secondary', 'Latent', 'Tertiary'], loc='lower right')
+    
+        ax[1].plot(sim.yearvec[pi:], sim.results.syphilis.prevalence[pi:])
+        ax[1].set_title('Syphilis prevalence')
+    
+        ax[2].plot(sim.yearvec[pi:], sim.results.n_alive[pi:])
+        ax[2].set_title('Population')
+    
+        ax[3].plot(sim.yearvec[pi:], sim.results.syphilis.new_infections[pi:])
+        ax[3].set_title('New infections')
+    
+        fig.tight_layout()
+        plt.show()
 
     return sim
 
@@ -171,10 +173,14 @@ def test_syph_intvs(dt=1, n_agents=500, do_plot=False):
 
 
 if __name__ == '__main__':
+    T = sc.timer()
+    do_plot = True
 
     dt = [1/12, 1][quick_run]
     n_agents = [20e3, 500][quick_run]
 
-    sim = test_syph(dt=dt, n_agents=n_agents)
+    sim = test_syph(dt=dt, n_agents=n_agents, do_plot=do_plot)
     sim_base, sim_intv = test_syph_intvs(dt=dt, n_agents=n_agents, do_plot=True)
+    
+    T.toc()