<jupyter_start><jupyter_text>SamplingIn this notebook, picking up where we left off in the "Getting started" tutorial, we are going to walk through a round of sampling. And we are even going to set up everything in a reproducable manner as good as possible.  Imports & SettingsFirst the imports.<jupyter_code>import numpy as np
import scipy as sp 
import pandas as pd

import emcee                      # inference and backends for sample storage
from multiprocessing import Pool  # for parallelization of the inference

import lymph<jupyter_output><empty_output><jupyter_text>Now some settings, e.g. the name of the HDF5 file we would later like to use.<jupyter_code>demo_hdf_file = "./_data/demo.hdf5"<jupyter_output><empty_output><jupyter_text>Set up the ModelFirst, we will set up the model as we would normally. In contrast to the "Getting started" notebook, we will set up a `Bilateral` model here, but that isn't more complicated. Only the data that needs to be provided to this kind of model needs to have information on the contralateral involvement as well, obviously.<jupyter_code>graph = {
    ('tumor', 'primary'): ['I', 'II', 'III', 'IV'],
    ('lnl'  , 'I'):       ['II'], 
    ('lnl'  , 'II'):      ['III'],
    ('lnl'  , 'III'):     ['IV'],
    ('lnl'  , 'IV'):      []
}
original_model = lymph.Bilateral(graph)
diagnostic_spsn = {
    "MRI": [0.63, 0.81],
    "PET": [0.86, 0.79]
}
original_model.modalities = diagnostic_spsn<jupyter_output><empty_output><jupyter_text>Generate synthetic data:::{note} This step can be skipped, as that data is already in the `./_data` directory. But it may also serve as a guide on how to generate synthetic datasets.:::<jupyter_code>max_t = 10
t = np.arange(max_t + 1)

early_p = 0.3
late_p = 0.7

early_time_dist = sp.stats.binom.pmf(t, max_t, early_p)
late_time_dist = sp.stats.binom.pmf(t, max_t, late_p)
time_dists = {"early": early_time_dist, "late": late_time_dist}

original_model.ipsi.base_probs   = [0.05, 0.2 , 0.12, 0.1 ]
original_model.contra.base_probs = [0.01, 0.06, 0.03, 0.01]
original_model.trans_probs = [0.1, 0.3, 0.2]
synthetic_data = original_model.generate_dataset(
    num_patients=200, 
    stage_dist=[0.6, 0.4], 
    time_dists=time_dists
)
synthetic_data.to_csv("./_data/bilateral.csv", index=False)<jupyter_output><empty_output><jupyter_text>Now we load the data into the model instance.<jupyter_code>synthetic_data = pd.read_csv("./_data/bilateral.csv", header=[0,1,2])
original_model.patient_data = synthetic_data<jupyter_output><empty_output><jupyter_text>Store the model in an HDF5 fileAnd before we proceed any further, we store the specifics of this model instance in an HDF5 file. It will basically store the graph, the modalities with their sensitivities & specificities as well as the just loaded data in the HDF5 file and allow us to recreate an instance.<jupyter_code>original_model.to_hdf(
    filename=demo_hdf_file, 
    name="original/model"
)<jupyter_output>/home/rmnldwg/repos/lymph/.venv/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'base_symmetric' in node 'model'. Offending HDF5 class: 8
  value = self._g_getattr(self._v_node, name)
/home/rmnldwg/repos/lymph/.venv/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'trans_symmetric' in node 'model'. Offending HDF5 class: 8
  value = self._g_getattr(self._v_node, name)<jupyter_text>Prepare the likelihood functionIn the utilities of the `lymph` package we also provide a small wrapper around the awesome [emcee](https://github.com/dfm/emcee) `EnsembleSampler` that allows us to store some inference-specific parameters before sampling and of course the samples themselves after sampling.Let's start with the first part:<jupyter_code># plus one dimension for the late T-stage's time parameter
ndim = len(original_model.spread_probs) + 1

# number of concurrent walkers that sample the space
nwalkers = 10 * ndim

# define the log-likelihood
def log_prob_fn(theta, sys, early_p=0.3, max_t=10):
    spread_probs, late_p = theta[:-1], theta[-1]
    
    if late_p > 1. or late_p < 0.:
        return -np.inf
    
    t = np.arange(max_t + 1)
    time_dists={
        "early": lymph.utils.fast_binomial_pmf(t, max_t, early_p),
        "late" : lymph.utils.fast_binomial_pmf(t, max_t, late_p)
    }
    
    return sys.marginal_log_likelihood(
        spread_probs, t_stages=["early", "late"], time_dists=time_dists
    )<jupyter_output><empty_output><jupyter_text>:::{admonition} Warning:class: warningThe provided log-likelihood function won't be stored anywhere! It is not possible to store arbitrary python code in an HDF5 file and retrieve it automatically in a safe manner.:::  SamplingFor storing the results, we make use of the `HDFBackend` from `emcee`, while the sampling itself can be done any way one pleases. However, we have written a sampling method `run_sampling` that smartly samples until convergence.:::{admonition} See also:class: noteThe creators of the `emcee` package have laid out how this "sampling to convergence" works in a [really nice tutorial](https://emcee.readthedocs.io/en/stable/tutorials/monitor/), which basically served as inspiration to the `run_sampling` method as well as our attempts of storing the model settings in an HDF5 file to begin with.:::<jupyter_code># this chain will surely be too short, but it doesn't matter here
max_steps = 200

# prepare the backend
backend = emcee.backends.HDFBackend(
    filename=demo_hdf_file,
    name="original/samples"
)
backend.reset(nwalkers, ndim)

# use Pool() from multiprocessing for parallelisation
with Pool() as pool:
    original_sampler = lymph.utils.EnsembleSampler(
        nwalkers, ndim, 
        log_prob_fn, args=[original_model], 
        pool=pool, backend=backend
    )
    acor_list = original_sampler.run_sampling(max_steps)<jupyter_output>Starting sampling<jupyter_text>We can make sure the chain of samples is actually stored by trying to retrieve it from the HDF5 file directly:<jupyter_code>test_backend = emcee.backends.HDFBackend(
    filename=demo_hdf_file,
    name="original/samples",
    read_only=True
)
test_backend.shape<jupyter_output><empty_output><jupyter_text>The first round has finished now. Let's see if we can repoduce all of that as intended. Do it all againWhen we load a model instance from the HDF5 storage, all the settings, i.e. the graph, the diagnostic modalities and the loaded data, should still be the same as in the beginning. So let's check that with some `assert`s.<jupyter_code>recovered_model = lymph.utils.system_from_hdf(
    filename=demo_hdf_file,
    name="original/model"
)
assert recovered_model.graph == graph, "Wrong graph!"
recovered_model.graph
assert recovered_model.modalities == diagnostic_spsn, "Wrong diagnostic modalities!"
recovered_model.modalities
assert np.all(recovered_model.patient_data == synthetic_data), "Wrong data!"
recovered_model.patient_data<jupyter_output><empty_output><jupyter_text>The recovery worked! Since we want to do another sampling round which itself should be reproducable as well, we can immediately store the recovered model in a new group of the HDF5 file.<jupyter_code>recovered_model.to_hdf(
    filename=demo_hdf_file,
    name="recovered/model"
)<jupyter_output>/home/rmnldwg/repos/lymph/.venv/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'base_symmetric' in node 'model'. Offending HDF5 class: 8
  value = self._g_getattr(self._v_node, name)
/home/rmnldwg/repos/lymph/.venv/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'trans_symmetric' in node 'model'. Offending HDF5 class: 8
  value = self._g_getattr(self._v_node, name)<jupyter_text>Now for the `EnsembleSampler`. We can recover the number of walkers and dimension from the previously stored HDF5 file. Note that I use two backends: One for retrieving the shape of the stored chain, which accesses the HDF5 group of the original sampling round. The however, I set up a new group for the second sampling round. The reason for this is that I don't want to call the `reset` method of the backend on my stored samples, thereby deleting them.<jupyter_code>tmp_backend = emcee.backends.HDFBackend(
    filename=demo_hdf_file,
    name="original/samples"
)
nwalkers, ndim = tmp_backend.shape

recovered_backend = emcee.backends.HDFBackend(
    filename=demo_hdf_file,
    name="recovered/samples"
)
recovered_backend.reset(nwalkers, ndim)

with Pool() as pool:
    recovered_sampler = lymph.utils.EnsembleSampler(
        nwalkers, ndim,
        log_prob_fn, args=[recovered_model],
        pool=pool, backend=recovered_backend
    )
    acor_list = recovered_sampler.run_sampling(max_steps)<jupyter_output>Starting sampling
{'fraction_non_alphanumeric': 0.0933046821985976, 'fraction_numerical': 0.012666817462112644, 'mean_word_length': 4.860172299536116, 'pattern_counts': {'":': 5, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '635915', 'license': 'MIT', 'n_tokens_mistral': 2975, 'n_tokens_neox': 2813, 'n_words': 893, 'path': 'docs/source/sampling.ipynb', 'repo_name': 'lfranceschetti/lymph', 'size': '8842'}