Using surmise’s global RNG (set_RNG)

1. Using surmise’s global RNG (set_RNG)#

surmise stores a single, user-provided random number generator for all of its internal sampling of random variables via scipy.stats.

  1. Users must call surmise.set_RNG(...) at least once before using any surmise functionality. Otherwise, an exception is raised.

  2. All emulation/calibration workflows are reproducible if the user repeats a workflow on the same system and provides surmise with the same RNG.

This notebook demonstrates these usages using the classic borehole example.

import numpy as np
import itertools as it
import scipy as sp
import scipy.stats as sps

import surmise
from surmise.emulation import emulator

print(f"surmise v{surmise.__version__}")
print(f"scipy   v{sp.__version__}")

surmise v0.1.dev1+g024e2a44d
scipy   v1.18.0

1.1. Borehole function setup#

The borehole function provides a fast-to-evaluate function as a classic emulation/calibration example. We make use of the function coded in surmise’s test suite. The only important detail here is that the function takes any parameters within the unit-cube, where \(x\in [0, 1]^3, \theta \in [0, 1]^4\).

Since this is a fresh session, surmise does not yet have an RNG set and any surmise call that needs randomness will raise an exception with instructions.

from surmise.tests.shared_scenario import (
    borehole_model,
    x_bh as x,
    thetatot_bh as theta
)

try:
    emulator(x=x, theta=theta, f=borehole_model(x, theta), method='PCGP')
except Exception as e:
    print(e)

Please use set_RNG before using surmise

1.2. set_RNG provides one generator#

Pass any type of RNG compatible with your scipy.stats installation (e.g., np.random.default_rng(seed)). Users are free to perform all other RNG draws in their application as they see fit, including using the same RNG set into surmise. The resulting workflow will be dependent on the collection of generator(s).

# import secrets
# SEED = secrets.randbits(128)

SEED = 111848137687551523431846058163015350939
SURMISE_SEED, DATA_SEED = np.random.SeedSequence(SEED).spawn(2)

# set RNG in surmise
surmise_rng = np.random.default_rng(SURMISE_SEED)
surmise.set_RNG(surmise_rng)

data_rng = np.random.default_rng(DATA_SEED)

x = sps.uniform.rvs(0, 1, size=(15, 3), random_state=data_rng)
theta = sps.uniform.rvs(0, 1, size=(50, 4), random_state=data_rng)

# surmise uses the surmise RNG under the hood
emu = emulator(x=x, theta=theta, f=borehole_model(x, theta), method='PCGP')
pred = emu.predict(x=x, theta=theta)
print("prediction mean shape:", pred.mean().shape)

prediction mean shape: (15, 50)

set_RNG does not accept old RNG types or seeds, but rather only the type of RNG compatible with the version of scipy.stats used to establish the surmise release.

print()
print(surmise_rng)
print(type(surmise_rng))
print()

for bad_rng in [np.random.RandomState(0), 12345]:
    try:
        surmise.set_RNG(bad_rng)
    except TypeError as e:
        print(f"{type(bad_rng).__name__!s:>12}: TypeError: {e}")

surmise.set_RNG(surmise_rng)  # restore the valid RNG

Generator(PCG64)
<class 'numpy.random._generator.Generator'>

 RandomState: TypeError: Given RNG cannot be used with scipy.stats
         int: TypeError: Given RNG cannot be used with scipy.stats

1.3. Using the same seed for reproducible results#

All randomness, namely both the user data generation and surmise’s, are controlled via these user-managed RNGs. Repeating a workflow with the same collection of RNGs should reproduce the results.

def run_workflow(seeds):
    seed, data_seed = seeds
    surmise.set_RNG(np.random.default_rng(seed))
    data_rng = np.random.default_rng(data_seed)
    x = sps.uniform.rvs(0, 1, size=(50, 3), random_state=data_rng)
    x[:, 2] = x[:, 2] > 0.5
    thetas = sps.uniform.rvs(0, 1, size=(15, 4), random_state=data_rng)
    emu = emulator(x=x, theta=thetas, f=borehole_model(x, thetas), method='PCGP')
    return emu.predict(x=x, theta=thetas).mean()

surmise_1, surmise_2, data_1, data_2 = np.random.SeedSequence(SEED).spawn(4)

runs = {
    "m1 (surmise 1, data 1)": run_workflow((surmise_1, data_1)),
    "m2 (surmise 1, data 1)": run_workflow((surmise_1, data_1)),
    "m3 (surmise 2, data 1)": run_workflow((surmise_2, data_1)),
    "m4 (surmise 1, data 2)": run_workflow((surmise_1, data_2)),
    "m5 (surmise 2, data 2)": run_workflow((surmise_2, data_2)),
}

for (run_name_a, run_result_a), (run_name_b, run_result_b) in it.combinations(runs.items(), 2):
    print(f"{run_name_a}, {run_name_b}: {np.array_equal(run_result_a, run_result_b)}")

m1 (surmise 1, data 1), m2 (surmise 1, data 1): True
m1 (surmise 1, data 1), m3 (surmise 2, data 1): False
m1 (surmise 1, data 1), m4 (surmise 1, data 2): False
m1 (surmise 1, data 1), m5 (surmise 2, data 2): False
m2 (surmise 1, data 1), m3 (surmise 2, data 1): False
m2 (surmise 1, data 1), m4 (surmise 1, data 2): False
m2 (surmise 1, data 1), m5 (surmise 2, data 2): False
m3 (surmise 2, data 1), m4 (surmise 1, data 2): False
m3 (surmise 2, data 1), m5 (surmise 2, data 2): False
m4 (surmise 1, data 2), m5 (surmise 2, data 2): False

1.4. Changing RNG at any time#

While surmise stores exactly one RNG at a time, users may change that RNG as needed for their study (e.g., to restart a numerical study from a known seed).

default_rng = np.random.default_rng(SEED)
surmise.set_RNG(default_rng)
# do something

pc64dxsm_rng = np.random.Generator(np.random.PCG64DXSM(SEED + 1))
surmise.set_RNG(pc64dxsm_rng)
# do something else