QInfer updater demo script

launch/local_launch.sh

@@ -4,7 +4,7 @@
 # QMLA run configuration
 ##### -------------------------------------------------- #####
 num_instances=1 # number of instances in run
-run_qhl=1 # perform QHL on known (true) model
+run_qhl=0 # perform QHL on known (true) model
 run_qhl_multi_model=0 # perform QHL for defined list of models
 experiments=500
 particles=2000
@@ -15,7 +15,7 @@ plot_level=1
 # Choose an exploration strategy 
 # This will determine how QMLA proceeds. 
 ##### -------------------------------------------------- #####
-exploration_strategy="HeisenbergTestDynamicsReproduction"
+exploration_strategy="DemoBayesFactorsByFscore"
 
 
 ##### -------------------------------------------------- #####

qmla/exploration_strategies/demos/bayes_factor_by_fscore.py

@@ -61,9 +61,9 @@ def __init__(self, exploration_rules, true_model=None, **kwargs):
         self.genetic_algorithm.terminate_early_if_top_model_unchanged = True
         self.branch_comparison_strategy = "all"
         self.tree_completed_initially = True
-        self.fraction_particles_for_bf = 1
-        self.fraction_own_experiments_for_bf = 1
-        self.fraction_opponents_experiments_for_bf = 1
+        self.fraction_particles_for_bf = 0.5
+        self.fraction_own_experiments_for_bf = 0.5
+        self.fraction_opponents_experiments_for_bf = 0.5
         self.iqle_mode = True
 
         self.max_time_to_consider = 15
@@ -101,6 +101,11 @@ def __init__(self, exploration_rules, true_model=None, **kwargs):
             "pauliSet_1J2_zJz_d5+pauliSet_1J3_zJz_d5+pauliSet_2J3_zJz_d5+pauliSet_2J5_zJz_d5+pauliSet_3J5_zJz_d5",  # F=1
         ]
 
+        test_fitness_models = [  # TODO reinstate longer list above for full plot
+            "pauliSet_1J2_zJz_d5+pauliSet_2J3_zJz_d5+pauliSet_2J5_zJz_d5+pauliSet_3J5_zJz_d5",
+            "pauliSet_1J2_zJz_d5+pauliSet_1J3_zJz_d5+pauliSet_2J3_zJz_d5+pauliSet_2J5_zJz_d5+pauliSet_3J5_zJz_d5",  # F=1
+        ]
+
         self.initial_models = list(
             np.random.choice(
                 test_fitness_models, len(test_fitness_models), replace=False

scripts/demo_qinfer_updater.py

@@ -0,0 +1,112 @@
+import numpy as np
+from copy import copy
+
+from qinfer import (
+    SMCUpdater,
+    FiniteOutcomeModel,
+    MultivariateNormalDistribution,
+    UniformDistribution,
+)
+
+
+def compute_prob_measuring_zero(particle, time):
+    r"""
+    Function to compute probability of measuring 0 in the gap experiment
+    """
+    return (2 / 3) * (1 + 0.5 * np.cos(particle * time))
+
+
+class DemoModel(FiniteOutcomeModel):
+    def __init__(self, **kwargs):
+        super(DemoModel).__init__(**kwargs)
+
+    def are_models_valid(self):
+        return True
+
+    @property
+    def expparams_dtype(self):
+        return [("t", "float")]
+
+    @property
+    def n_modelparams(self):
+        return 1
+
+    @property
+    def n_outcomes(self):
+        return 2
+
+    def likelihood(self, outcomes, modelparams, expparams):
+
+        # modelparams give particles
+        # outcomes should be a length 1 list with a single outcome for the single experiment performed
+
+        time = expparams["t"]
+        print(
+            f"Likelihood fnc given \nOutcomes:\n {outcomes} \nParticles:\n {modelparams} \nExperiment:\n {expparams}"
+        )
+
+        prob_measuring_zero = compute_prob_measuring_zero(
+            modelparams, time=expparams["t"]  # particles
+        )
+
+        likelihood = FiniteOutcomeModel.pr0_to_likelihood_array(
+            outcomes=outcomes,
+            pr0=prob_measuring_zero,
+        )
+
+        return likelihood
+
+
+# Make instance of our model and generate an SMCUpdater from it
+demo_model = DemoModel()
+uniform_distribution = UniformDistribution(ranges=[0, 1])
+
+updater = SMCUpdater(model=demo_model, n_particles=3, prior=uniform_distribution)
+
+# Track particles to inspect before/after update
+track_particles = [copy(updater.particle_locations)]
+track_weights = [copy(updater.particle_weights)]
+
+# Design experiment
+experiment = np.empty((1,), dtype=demo_model.expparams_dtype)
+experiment["t"] = 7.69
+
+# Update distribution
+updater.update(outcome=0, expparams=experiment)
+track_particles.append(updater.particle_locations)
+track_weights.append(updater.particle_weights)
+
+
+# Within update, it calls hypothetical_update - we can inspect its outputs here
+hypothetical_weights, likelihood_array = updater.hypothetical_update(
+    outcomes=np.array([0]),
+    expparams=experiment,
+    return_likelihood=True,
+)
+likelihood_array = likelihood_array.reshape(
+    3,
+)
+
+
+# Inspect particle calculations
+
+normalisation_factor = np.sum(likelihood_array * track_weights[0])
+
+# check the calculation of norm factor is correct
+if not np.isclose(normalisation_factor, updater.normalization_record[0][0]):
+    print("Normalisation record is not the same as expected")
+
+for particle_idx in range(3):
+
+    likelihood = likelihood_array[particle_idx]
+    weight_before = track_weights[0][particle_idx]
+    weight_after = track_weights[1][particle_idx]
+
+    check_correct = np.isclose(
+        weight_after,
+        (weight_before * likelihood) / normalisation_factor,
+        atol=1e-8,  # due to printing accuracy used in likelihood
+    )
+
+    if not check_correct:
+        print(f"Calculation not the same at particle idx {particle_idx}")