particle updater tests

Author: KOVALW

src/calibrationtools/__init__.py

@@ -1,5 +1,6 @@
 from .particle import Particle
 from .particle_population import ParticlePopulation
+from .particle_updater import _ParticleUpdater
 from .perturbation_kernel import (
     IndependentKernels,
     MultivariateNormalKernel,
@@ -16,6 +17,7 @@
 )
 from .sampler import ABCSampler
 from .variance_adapter import (
+    AdaptIdentityVariance,
     AdaptMultivariateNormalVariance,
     AdaptNormalVariance,
     VarianceAdapter,
@@ -25,6 +27,7 @@
     "ABCSampler",
     "Particle",
     "ParticlePopulation",
+    "_ParticleUpdater",
     "PriorDistribution",
     "UniformPrior",
     "SeedPrior",
@@ -38,4 +41,5 @@
     "VarianceAdapter",
     "AdaptNormalVariance",
     "AdaptMultivariateNormalVariance",
+    "AdaptIdentityVariance",
 ]

src/calibrationtools/particle_updater.py

@@ -21,9 +21,6 @@ def __init__(
         self.priors = priors
         self.variance_adapter = variance_adapter
         self.seed_sequence = seed_sequence
-        self.calculate_transition_probability = np.vectorize(
-            self.perturbation_kernel.transition_probability
-        )
 
     def set_particle_population(self, particle_population: ParticlePopulation):
         self.particle_population = particle_population
@@ -58,8 +55,13 @@ def sample_perturbed_particle(
 
     def calculate_weight(self, particle: Particle) -> float:
         numerator = self.priors.probability_density(particle)
-        transition_probs = self.calculate_transition_probability(
-            particle, self.particle_population.particles
+        transition_probs = np.array(
+            [
+                self.perturbation_kernel.transition_probability(
+                    to_particle=particle, from_particle=p
+                )
+                for p in self.particle_population.particles
+            ]
         )
 
         denominator = np.sum(
@@ -76,6 +78,3 @@ def adapt_variance(self):
         self.variance_adapter.adapt(
             self.particle_population, self.perturbation_kernel
         )
-        self.calculate_transition_probability = np.vectorize(
-            self.perturbation_kernel.transition_probability
-        )

src/calibrationtools/sampler.py

@@ -63,7 +63,7 @@ def run(self):
         for generation in range(len(self.tolerance_values)):
             if self.verbose:
                 print(
-                    f"Running generation {generation + 1} with tolerance {self.tolerance_values[generation]}..."
+                    f"Running generation {generation + 1} with tolerance {self.tolerance_values[generation]}... Previous population is {self.particle_population}"
                 )
 
             # Rejection sampling algorithm

tests/conftest.py

@@ -2,9 +2,11 @@
 from numpy.random import SeedSequence
 
 from calibrationtools import (
+    AdaptIdentityVariance,
     IndependentKernels,
     NormalKernel,
     Particle,
+    ParticlePopulation,
     SeedKernel,
 )
 from calibrationtools.prior_distribution import (
@@ -66,6 +68,21 @@ def Kc() -> IndependentKernels:
     return Kc
 
 
+@pytest.fixture
+def particle_population() -> ParticlePopulation:
+    particles = [
+        {"p": 0.1, "seed": 0},
+        {"p": 0.5, "seed": 1},
+        {"p": 0.9, "seed": 2},
+    ]
+    return ParticlePopulation(states=particles, weights=[0.2, 0.3, 0.5])
+
+
+@pytest.fixture
+def proposed_particle() -> Particle:
+    return Particle({"p": 0.6, "seed": 3})
+
+
 @pytest.fixture
 def N() -> int:
     return 10
@@ -89,3 +106,8 @@ def P() -> IndependentPriors:
 @pytest.fixture
 def Pc() -> IndependentPriors:
     return IndependentPriors([UniformPrior("p", 0.0, 1.0), SeedPrior("seed")])
+
+
+@pytest.fixture
+def V() -> AdaptIdentityVariance:
+    return AdaptIdentityVariance()

tests/test_particle_updater.py

@@ -0,0 +1,74 @@
+import pytest
+from calibrationtools import ParticlePopulation, _ParticleUpdater
+
+@pytest.fixture
+def particle_updater(seed_sequence, K, P, V):
+    updater = _ParticleUpdater(K, P, V, seed_sequence)
+    return updater
+
+def test_set_particle_population_normalizes_weights(particle_updater, particle_population):
+    # Set the particle population in the updater
+    particle_updater.set_particle_population(particle_population)
+
+    # Check that the weights are normalized
+    assert particle_updater.particle_population.total_weight == pytest.approx(1.0)
+    
+    particle_population_unnormalized = ParticlePopulation(states=particle_population.particles, weights=[0.1, 0.1, 0.1])  # Not normalized
+    particle_updater.set_particle_population(particle_population_unnormalized)
+    
+    # Check that the weights are normalized
+    assert particle_updater.particle_population.total_weight == pytest.approx(1.0)
+    assert particle_updater.particle_population.weights == pytest.approx([1/3, 1/3, 1/3])
+
+def test_sample_particle(particle_updater, particle_population):
+    particle_updater.set_particle_population(particle_population)
+    sampled_particle = particle_updater.sample_particle()
+    assert sampled_particle in particle_population.particles
+
+def test_sample_perturbed_particle(particle_updater, particle_population):
+    particle_updater.set_particle_population(particle_population)
+    perturbed_particle = particle_updater.sample_perturbed_particle()
+    assert perturbed_particle not in particle_population.particles  # Perturbed particle should not be the same as any in the population
+    assert particle_updater.priors.probability_density(perturbed_particle) > 0  # Perturbed particle should have non-zero prior density
+
+def test_sample_perturbed_particle_max_attempts(particle_updater, particle_population):
+    # Create a perturbation kernel that always produces invalid particles
+    class InvalidPerturbationKernel:
+        def perturb(self, current_particle, seed_sequence):
+            return {"p": -1.0, "seed": 0}  # Invalid particle outside the prior support
+
+    particle_updater.perturbation_kernel = InvalidPerturbationKernel()
+    particle_updater.set_particle_population(particle_population)
+
+    with pytest.raises(RuntimeError):
+        particle_updater.sample_perturbed_particle(max_attempts=5)
+
+def test_calculate_weight(particle_updater, particle_population, proposed_particle):
+    particle_updater.set_particle_population(particle_population)
+    weight = particle_updater.calculate_weight(proposed_particle)
+    assert weight >= 0  # Weights should be non-negative
+
+    # Check that weight was calculated correctly for nomrla perturbation
+    states = particle_population.particles
+    weights = particle_population.weights
+    transition_probs = [
+        particle_updater.perturbation_kernel.transition_probability(to_particle=proposed_particle, from_particle=p)
+        for p in states
+    ]
+    weighted_probs = [w*p for w, p in zip(weights, transition_probs)]
+
+    expected_weight = particle_updater.priors.probability_density(proposed_particle) / sum(weighted_probs)
+    assert weight == pytest.approx(expected_weight)
+
+def test_calculate_weight_zero_prob_perturbation(particle_updater, particle_population):
+    # Create a perturbation kernel that always produces zero transition probability
+    class ZeroTransitionPerturbationKernel:
+        def transition_probability(self, to_particle, from_particle):
+            return 0.0  # Zero transition probability
+
+    particle_updater.perturbation_kernel = ZeroTransitionPerturbationKernel()
+    particle_updater.set_particle_population(particle_population)
+
+    proposed_particle = {"p": 0.5, "seed": 1}  # A valid particle with non-zero prior density
+    weight = particle_updater.calculate_weight(proposed_particle)
+    assert weight == 0.0  # Weight should be zero due to zero transition probabilities