Added MRP for example

.gitignore

@@ -36,7 +36,6 @@
 # !your_data_file.csv
 # !your_data_directory/
 
-
 #####
 # Python
 # https://github.com/github/gitignore/blob/main/Python.gitignore
@@ -326,6 +325,8 @@ replay_pid*
 
 #Emacs
 *~
+*.org
+
 #####
 # Rust
 # https://github.com/github/gitignore/blob/main/Rust.gitignore

packages/example_model/README.md

@@ -11,3 +11,9 @@ uv run python -m example_model.example_model packages/example_model/tests/model_
 
 Which will write `model_output.json` with the results of the run.  You
 can specify a different location with the `-o` command line argument.
+
+To run the calibration example, run
+
+```bash
+uv run python -m example_model.calibrate
+```

packages/example_model/src/example_model/calibrate.py

@@ -0,0 +1,109 @@
+"""Calibrate the example branching process."""
+
+import numpy as np
+from mrp import Environment
+from mrp.api import apply_dict_overrides
+
+from calibrationtools.perturbation_kernel import (
+    IndependentKernels,
+    MultivariateNormalKernel,
+    SeedKernel,
+)
+from calibrationtools.prior_distribution import (
+    IndependentPriors,
+    SeedPrior,
+    UniformPrior,
+)
+from calibrationtools.sampler import ABCSampler
+from calibrationtools.variance_adapter import AdaptMultivariateNormalVariance
+from example_model import Binom_BP_Model
+
+##===================================#
+## Define model
+##===================================#
+env = Environment(
+    {
+        "input": {
+            "seed": 123,
+            "max_gen": 15,
+            "n": 3,
+            "p": 0.5,
+            "max_infect": 500,
+        },
+        "output": {"spec": "filesystem", "dir": "./output"},
+    }
+)
+default_inputs = {
+    "seed": 123,
+    "max_gen": 15,
+    "n": 3,
+    "p": 0.5,
+    "max_infect": 500,
+}
+model = Binom_BP_Model(env=env)
+
+##===================================#
+## Define priors
+##===================================#
+P = IndependentPriors(
+    [
+        UniformPrior("n", 0, 5),
+        UniformPrior("p", 0, 1),
+        SeedPrior("seed"),
+    ]
+)
+
+K = IndependentKernels(
+    [
+        MultivariateNormalKernel(
+            [p.params[0] for p in P.priors if not isinstance(p, SeedPrior)],
+        ),
+        SeedKernel("seed"),
+    ]
+)
+
+V = AdaptMultivariateNormalVariance()
+
+
+##===================================#
+## Run ABC-SMC
+##===================================#
+def particles_to_params(particle, **kwargs):
+    base_inputs = kwargs.get("base_inputs")
+    model_params = apply_dict_overrides(base_inputs, particle)
+    return model_params
+
+
+def outputs_to_distance(model_output, target_data, **kwargs):
+    return abs(np.sum(model_output) - target_data)
+
+
+sampler = ABCSampler(
+    generation_particle_count=500,
+    tolerance_values=[5.0, 1.0],
+    priors=P,
+    perturbation_kernel=K,
+    variance_adapter=V,
+    particles_to_params=particles_to_params,
+    outputs_to_distance=outputs_to_distance,
+    target_data=5,
+    model_runner=model,
+    seed=123,  # Propagation of seed must be SeedSequence not int for proper pseudorandom draws
+)
+
+sampler.run(base_inputs=default_inputs)
+
+##===================================#
+## Get results
+##===================================#
+# Print IQR of param1 in the posterior particles
+posterior_particles = sampler.get_posterior_particles()
+p_values = [p["p"] for p in posterior_particles.particles]
+n_values = [p["n"] for p in posterior_particles.particles]
+
+print(
+    f"param p(25-75):{np.percentile(p_values, 25)} - {np.percentile(p_values, 75)}"
+)
+print(
+    f"param n(25-75):{np.percentile(n_values, 25)} - {np.percentile(n_values, 75)}"
+)

src/calibrationtools/sampler.py

@@ -56,7 +56,13 @@ def particle_population(self) -> ParticlePopulation:
     def particle_population(self, population: ParticlePopulation):
         self._updater.set_particle_population(population)
 
-    def run(self):
+    def run(self, **kwargs: Any):
+        for k in kwargs.keys():
+            if k in self.__class__.__dict__:
+                raise ValueError(
+                    f"Keyword argument '{k}' conflicts with existing attribute. Please choose a different name for the argument. Attributes cannot be set from `.run()`"
+                )
+
         self.particle_population = self.sample_particles_from_priors()
         proposed_population = ParticlePopulation()
 
@@ -77,7 +83,7 @@ def run(self):
                 attempts += 1
                 # Create the parameter inputs for the runner by sampling perturbed value from previous population
                 proposed_particle = self.sample_proposed_particle()
-                params = self.particles_to_params(proposed_particle)
+                params = self.particles_to_params(proposed_particle, **kwargs)
 
                 # Generate the distance metric from model run
                 outputs = self.model_runner.simulate(params)