Rework to make use of frozen SciPy distributions

[daniel-klein]

This PR includes major updates.

SciPy distributions as parameters

We often used parameters that were later fed into distributions. For example, a parameter named mean_partners could be later fed into a normal distribution to determine the number of partners that each individual would seek.

self.pars['mean_partners'] = 5
...
self.partners = rng.normal(mean=self.pars['mean_partners'], std=3)

We had a nascent effort to create our own "frozen" distribution class to enable

import stisim as ss

self.pars['partner_distrib'] = ss.Normal(mean=5, std=3)
...
self.partners = self.pars.partner_distrib.sample(uids)

But this required us to wrap each numpy distribution individually.

New in this PR we can directly use frozen distributions from SciPy:

import scipy.stats as sps
self.pars['mean_partners'] = sps.norm(loc=5, scale=3)
...
self.partners = self.pars.partner_distrib.rvs(uids)

Callable parameters!

This is so cool... the parameters of SciPy distributions can be callable! This enables users to easily make any distribution parameter vary as function of anything else. Callables can be lambda or functions. For example (from test_base.py):

def test_ppl_construction():

    def init_debut(self, sim, uids):
        loc = np.full(len(uids), 16) # Mean of 16 for male
        loc[sim.people.female[uids]] = 21 # Mean of 21 for female
        return loc

    mf_pars = {
        'debut_dist': sps.norm(loc=init_debut, scale=2),  # Age of debut is normally distributed with mean that varies by sex
    }
    sim_pars = {'networks': [ss.mf(mf_pars)], 'n_agents': 100}
    gon_pars = {'beta': {'mf': [0.08, 0.04]}, 'p_death': 0.2}
    gon = ss.Gonorrhea(pars=gon_pars)

    sim = ss.Sim(pars=sim_pars, diseases=[gon])
    sim.initialize()
    sim.run()
    plt.figure()
    plt.plot(sim.tivec, sim.results.gonorrhea.n_infected)
    plt.title('Number of gonorrhea infections')

    return sim

Here you can see that the call signature is self, sim, uids. While a static method, self is the context of the calling instance, in this case an instance of mf - this is super handy! Users also gain access to the sim and the uids of the agents in question. The return value should be a scalar or array of length len(uids).

This is very powerful! Users can make parameters vary by time, perhaps using interpolation of data or other functional forms, or by any other attribute in the simulation. The callable has access to data in self, in case data or other information is resident on the class object. (This is really just a shortcut to allow self.cd4 in place of e.g. sim.people.hiv.cd4. I am still finding clever ways to use this functionality, one is to enable users to directly override specific scientific methods, like this example from hiv.py:

class HIV(STI):
    def __init__(self, pars=None):
        ...
        self.death_prob_per_dt = sps.bernoulli(p=self.death_prob)
        ...
        return

    @staticmethod
    def death_prob(self, sim, uids):
        return 0.05 / (self.pars.cd4_min - self.pars.cd4_max)**2 *  (self.cd4[uids] - self.pars.cd4_max)**2

Also, this is amazing for initialization. While the default in many modules is sps.bernoulli on the whole population, users could quickly "subtarget" a desired subgroup without changing model code, subclassing, or writing a custom "intervention" acting at time 0 only.

Oh, and users don't have to completely overwrite the distribution. From text_examples.py:

sir.pars['dur_inf'].kwds['loc'] = 5

Finally, I expect this new functionality to be super powerful when used in combination with extra_states.

Common random numbers

SciPy distributions in __dict__ and Module pars are automatically connected to a MultiRNG or SingleRNG automatically. The MultiRNG is now thin as it is now only primarily responsible for setting the seed, jumping, and resets. The SingleRNG is improved as well, and now uses the numpy random singleton directly. The user no longer needs to create random number generators and doesn't even have to specific the (unique) RNG names as these too are determined automatically using the .name and .label. Subclasses should change the label. Please note #167.

Bernoulli filtering

For Bernoulli draws, we previously had the bernoulli_filter. This functionality is now available on top of SciPy distributions:

import scipy.stats as sps
self.will_recover = sps.bernoulli(p=0.8)
...
recover_uids = self.will_recover.filter(uids)

Automatic initialization

I mention this briefly above, but distributions of type rv_frozen are automatically initialized and connected to appropriate RNGs. However, this only works if they are in __dict__ or pars. Other distributions would need to be manually initialized.

Status

Still a work in progress. As of the current state, the code will print lots of "DJK TODO" messages where more attention is needed. I'm working on it! FEEDBACK WELCOME, this is a big one, so let's discuss.