Advanced History Matching Configuration¶
This notebook demonstrates advanced configuration options and customisation capabilities of the history matching library.
Overview¶
Building on the basic and manual workflows, we'll explore:
- Single-constructor configuration and DataFrame-based inputs
- Custom sampling strategies and emulator types
- Feature selection strategies and switching between iterations
- Progress callbacks and monitoring
- Performance comparison across configurations
- Checkpoint and resume functionality
We reuse the stochastic SIR model from model.py. For the basic automated workflow see
01_basic_workflow.ipynb; for the manual step-by-step workflow see 02_manual_workflow.ipynb.
In [1]:
Copied!
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import historymatching as hm
from model import SIR, generate_observed_data
%matplotlib inline
np.random.seed(42)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import historymatching as hm
from model import SIR, generate_observed_data
%matplotlib inline
np.random.seed(42)
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR I0000 00:00:1780533786.485116 1037745 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`. I0000 00:00:1780533786.547206 1037745 cpu_feature_guard.cc:227] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: AVX2 AVX_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR I0000 00:00:1780533787.792721 1037745 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
/software/conda/envs/hm/lib/python3.13/site-packages/gpflow/versions.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81. import pkg_resources
Generate Synthetic Observed Data¶
In [2]:
Copied!
beta_true = 1.3
gamma_true = 0.5
population_size = 10_000
initial_infected = 100
incidence_obs, true_model = generate_observed_data(
beta_true=beta_true,
gamma_true=gamma_true,
population_size=population_size,
n_seed_infections=initial_infected,
)
incidence_obs = incidence_obs.values # numpy array for convenience below
print(f"True parameters: β={beta_true}, γ={gamma_true}, R₀={beta_true/gamma_true:.2f}")
print(f"Peak incidence: {incidence_obs.max():.0f}, Total cases: {incidence_obs.sum():.0f}")
beta_true = 1.3
gamma_true = 0.5
population_size = 10_000
initial_infected = 100
incidence_obs, true_model = generate_observed_data(
beta_true=beta_true,
gamma_true=gamma_true,
population_size=population_size,
n_seed_infections=initial_infected,
)
incidence_obs = incidence_obs.values # numpy array for convenience below
print(f"True parameters: β={beta_true}, γ={gamma_true}, R₀={beta_true/gamma_true:.2f}")
print(f"Peak incidence: {incidence_obs.max():.0f}, Total cases: {incidence_obs.sum():.0f}")
True parameters: β=1.3, γ=0.5, R₀=2.60 Peak incidence: 1576, Total cases: 9091
In [3]:
Copied!
def enhanced_sir_simulation(samples: pd.DataFrame) -> pd.DataFrame:
"""
Enhanced simulation function that returns more features.
"""
results = []
for _, row in samples.iterrows():
model = SIR(
beta=row['beta'],
gamma=row['gamma'],
s0=population_size - initial_infected,
i0=initial_infected
)
S, I, R = model.run()
incidence = model.get_incidence()
# Create comprehensive feature set
result = {}
# Time series features (every other day to reduce dimensionality)
for i in range(0, len(incidence), 2):
result[f'incidence_day_{i}'] = incidence[i]
result[f'infected_day_{i}'] = I[i]
# Summary statistics
result['peak_incidence'] = max(incidence)
result['peak_incidence_day'] = np.argmax(incidence)
result['total_cases'] = sum(incidence)
result['attack_rate'] = sum(incidence) / population_size
result['final_size'] = R[-1]
result['max_infected'] = max(I)
# Epidemic characteristics
result['time_to_peak'] = np.argmax(incidence)
result['epidemic_duration'] = len([x for x in incidence[1:] if x > 0])
result['growth_rate'] = np.log(incidence[3] / incidence[1]) if incidence[1] > 0 and incidence[3] > 0 else 0
# Cumulative features
cumulative_cases = np.cumsum(incidence)
result['cases_week_1'] = cumulative_cases[6] if len(cumulative_cases) > 6 else cumulative_cases[-1]
result['cases_week_2'] = cumulative_cases[13] if len(cumulative_cases) > 13 else cumulative_cases[-1]
results.append(result)
return pd.DataFrame(results)
# Test the enhanced simulation
test_samples = pd.DataFrame({
'beta': [1.0, 1.5],
'gamma': [0.3, 0.5]
})
test_results = enhanced_sir_simulation(test_samples)
print(f"Enhanced simulation produces {len(test_results.columns)} features:")
print(f"Features: {list(test_results.columns[:10])}...")
def enhanced_sir_simulation(samples: pd.DataFrame) -> pd.DataFrame:
"""
Enhanced simulation function that returns more features.
"""
results = []
for _, row in samples.iterrows():
model = SIR(
beta=row['beta'],
gamma=row['gamma'],
s0=population_size - initial_infected,
i0=initial_infected
)
S, I, R = model.run()
incidence = model.get_incidence()
# Create comprehensive feature set
result = {}
# Time series features (every other day to reduce dimensionality)
for i in range(0, len(incidence), 2):
result[f'incidence_day_{i}'] = incidence[i]
result[f'infected_day_{i}'] = I[i]
# Summary statistics
result['peak_incidence'] = max(incidence)
result['peak_incidence_day'] = np.argmax(incidence)
result['total_cases'] = sum(incidence)
result['attack_rate'] = sum(incidence) / population_size
result['final_size'] = R[-1]
result['max_infected'] = max(I)
# Epidemic characteristics
result['time_to_peak'] = np.argmax(incidence)
result['epidemic_duration'] = len([x for x in incidence[1:] if x > 0])
result['growth_rate'] = np.log(incidence[3] / incidence[1]) if incidence[1] > 0 and incidence[3] > 0 else 0
# Cumulative features
cumulative_cases = np.cumsum(incidence)
result['cases_week_1'] = cumulative_cases[6] if len(cumulative_cases) > 6 else cumulative_cases[-1]
result['cases_week_2'] = cumulative_cases[13] if len(cumulative_cases) > 13 else cumulative_cases[-1]
results.append(result)
return pd.DataFrame(results)
# Test the enhanced simulation
test_samples = pd.DataFrame({
'beta': [1.0, 1.5],
'gamma': [0.3, 0.5]
})
test_results = enhanced_sir_simulation(test_samples)
print(f"Enhanced simulation produces {len(test_results.columns)} features:")
print(f"Features: {list(test_results.columns[:10])}...")
Enhanced simulation produces 31 features: Features: ['incidence_day_0', 'infected_day_0', 'incidence_day_2', 'infected_day_2', 'incidence_day_4', 'infected_day_4', 'incidence_day_6', 'infected_day_6', 'incidence_day_8', 'infected_day_8']...
Advanced configuration¶
The single hm.HistoryMatching(...) constructor provides fine-grained control over all aspects of the history matching workflow.
In [4]:
Copied!
# Create parameter space using DataFrame for more control
parameter_df = pd.DataFrame({
'parameter': ['beta', 'gamma'],
'minimum': [0.5, 0.1],
'maximum': [3.0, 1.0],
'description': ['Transmission rate per day', 'Recovery rate per day']
})
# Create comprehensive observations - mix of summary stats and time series
observations_df = pd.DataFrame({
'feature': [
'peak_incidence',
'total_cases',
'attack_rate',
'time_to_peak',
'incidence_day_4',
'incidence_day_8',
'incidence_day_12',
'cases_week_1',
'cases_week_2'
],
'mean': [
max(incidence_obs),
sum(incidence_obs),
sum(incidence_obs) / population_size,
np.argmax(incidence_obs),
incidence_obs[4],
incidence_obs[8],
incidence_obs[12],
sum(incidence_obs[:7]),
sum(incidence_obs[:14])
],
'std': [
50, # Peak incidence uncertainty
200, # Total cases uncertainty
0.02, # Attack rate uncertainty
1, # Time to peak uncertainty (days)
30, # Daily incidence uncertainty
40,
25,
100, # Weekly totals uncertainty
150
]
})
print(f"Advanced parameter space ({len(parameter_df)} parameters):")
print(parameter_df)
print(f"\nComprehensive observations ({len(observations_df)} features):")
print(observations_df)
# Create parameter space using DataFrame for more control
parameter_df = pd.DataFrame({
'parameter': ['beta', 'gamma'],
'minimum': [0.5, 0.1],
'maximum': [3.0, 1.0],
'description': ['Transmission rate per day', 'Recovery rate per day']
})
# Create comprehensive observations - mix of summary stats and time series
observations_df = pd.DataFrame({
'feature': [
'peak_incidence',
'total_cases',
'attack_rate',
'time_to_peak',
'incidence_day_4',
'incidence_day_8',
'incidence_day_12',
'cases_week_1',
'cases_week_2'
],
'mean': [
max(incidence_obs),
sum(incidence_obs),
sum(incidence_obs) / population_size,
np.argmax(incidence_obs),
incidence_obs[4],
incidence_obs[8],
incidence_obs[12],
sum(incidence_obs[:7]),
sum(incidence_obs[:14])
],
'std': [
50, # Peak incidence uncertainty
200, # Total cases uncertainty
0.02, # Attack rate uncertainty
1, # Time to peak uncertainty (days)
30, # Daily incidence uncertainty
40,
25,
100, # Weekly totals uncertainty
150
]
})
print(f"Advanced parameter space ({len(parameter_df)} parameters):")
print(parameter_df)
print(f"\nComprehensive observations ({len(observations_df)} features):")
print(observations_df)
Advanced parameter space (2 parameters):
parameter minimum maximum description
0 beta 0.5 3.0 Transmission rate per day
1 gamma 0.1 1.0 Recovery rate per day
Comprehensive observations (9 features):
feature mean std
0 peak_incidence 1576.0000 50.00
1 total_cases 9091.0000 200.00
2 attack_rate 0.9091 0.02
3 time_to_peak 5.0000 1.00
4 incidence_day_4 1246.0000 30.00
5 incidence_day_8 810.0000 40.00
6 incidence_day_12 130.0000 25.00
7 cases_week_1 5696.0000 100.00
8 cases_week_2 8939.0000 150.00
In [5]:
Copied!
# Create advanced history matching with full configuration control
print("Building advanced history matching engine...")
advanced_engine = hm.HistoryMatching(
bounds=parameter_df,
observations=observations_df,
function=enhanced_sir_simulation,
sampling_strategy={
'type': 'lhs',
'criterion': 'maximin', # Maximize minimum distance between points
'iterations': 10 # LHS optimization iterations
},
feature_selection={
'method': 'fano', # Fano factor-based selection
'max_features': 3, # Limit features per iteration
'correlation_threshold': 0.8 # Avoid highly correlated features
},
emulator_type='gpr',
n_samples=800,
max_iterations=6,
implausibility_threshold=2.8,
auto_reduce_space=False,
oversample_factor=5.0,
random_seed=789,
)
# Summarize the configuration
print("Configuration summary:")
print(f" parameters: {advanced_engine.parameters}")
print(f" observations: {advanced_engine.outputs}")
print(f" emulator type: gpr")
print(f" n_samples: 800")
print(f" max_iterations: 6")
print(f" implausibility_threshold: 2.8")
print(f" oversample_factor: 5.0")
# Create advanced history matching with full configuration control
print("Building advanced history matching engine...")
advanced_engine = hm.HistoryMatching(
bounds=parameter_df,
observations=observations_df,
function=enhanced_sir_simulation,
sampling_strategy={
'type': 'lhs',
'criterion': 'maximin', # Maximize minimum distance between points
'iterations': 10 # LHS optimization iterations
},
feature_selection={
'method': 'fano', # Fano factor-based selection
'max_features': 3, # Limit features per iteration
'correlation_threshold': 0.8 # Avoid highly correlated features
},
emulator_type='gpr',
n_samples=800,
max_iterations=6,
implausibility_threshold=2.8,
auto_reduce_space=False,
oversample_factor=5.0,
random_seed=789,
)
# Summarize the configuration
print("Configuration summary:")
print(f" parameters: {advanced_engine.parameters}")
print(f" observations: {advanced_engine.outputs}")
print(f" emulator type: gpr")
print(f" n_samples: 800")
print(f" max_iterations: 6")
print(f" implausibility_threshold: 2.8")
print(f" oversample_factor: 5.0")
Building advanced history matching engine... Configuration summary: parameters: ['beta', 'gamma'] observations: ['peak_incidence', 'total_cases', 'attack_rate', 'time_to_peak', 'incidence_day_4', 'incidence_day_8', 'incidence_day_12', 'cases_week_1', 'cases_week_2'] emulator type: gpr n_samples: 800 max_iterations: 6 implausibility_threshold: 2.8 oversample_factor: 5.0
In [6]:
Copied!
# The engine was configured and built in the previous cell
print(f" Advanced engine built successfully!")
print(f" Parameters: {len(advanced_engine.parameter_space.get_parameter_names())}")
print(f" Observations: {len(observations_df)}")
print(f" Sampling strategy: {advanced_engine.sampling_strategy.get_strategy_name()}")
print(f" Feature selection: {advanced_engine.feature_selection_strategy.get_strategy_name()}")
print(f" Emulator type: {type(advanced_engine.emulator_factory).__name__}")
print(f" Implausibility threshold: {advanced_engine.implausibility_threshold}")
print(f" Oversample factor: {advanced_engine.oversample_factor}")
# The engine was configured and built in the previous cell
print(f" Advanced engine built successfully!")
print(f" Parameters: {len(advanced_engine.parameter_space.get_parameter_names())}")
print(f" Observations: {len(observations_df)}")
print(f" Sampling strategy: {advanced_engine.sampling_strategy.get_strategy_name()}")
print(f" Feature selection: {advanced_engine.feature_selection_strategy.get_strategy_name()}")
print(f" Emulator type: {type(advanced_engine.emulator_factory).__name__}")
print(f" Implausibility threshold: {advanced_engine.implausibility_threshold}")
print(f" Oversample factor: {advanced_engine.oversample_factor}")
Advanced engine built successfully! Parameters: 2 Observations: 9 Sampling strategy: Latin Hypercube Sampling (criterion=maximin)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[6], line 6 2 print(f" Advanced engine built successfully!") 3 print(f" Parameters: {len(advanced_engine.parameter_space.get_parameter_names())}") 4 print(f" Observations: {len(observations_df)}") 5 print(f" Sampling strategy: {advanced_engine.sampling_strategy.get_strategy_name()}") ----> 6 print(f" Feature selection: {advanced_engine.feature_selection_strategy.get_strategy_name()}") 7 print(f" Emulator type: {type(advanced_engine.emulator_factory).__name__}") 8 print(f" Implausibility threshold: {advanced_engine.implausibility_threshold}") 9 print(f" Oversample factor: {advanced_engine.oversample_factor}") AttributeError: 'HistoryMatching' object has no attribute 'feature_selection_strategy'
Interactive Workflow with Advanced Controls¶
Let's use the interactive workflow to demonstrate advanced features like strategy switching and diagnostics.
In [7]:
Copied!
# Add progress callback to monitor the workflow
def progress_callback(progress):
print(f" Progress update: Iteration {progress.current_iteration}, "
f"Acceptance rate: {progress.acceptance_rate:.3f}, "
f"Total emulators: {progress.total_emulators_trained}")
advanced_engine.add_progress_callback(progress_callback)
print("Starting interactive advanced workflow...")
print(f"Initial parameter ranges:")
for param_name in advanced_engine.parameter_space.get_parameter_names():
bounds = advanced_engine.parameter_space.get_bounds(param_name)
print(f" {param_name}: [{bounds[0]}, {bounds[1]}]")
# Add progress callback to monitor the workflow
def progress_callback(progress):
print(f" Progress update: Iteration {progress.current_iteration}, "
f"Acceptance rate: {progress.acceptance_rate:.3f}, "
f"Total emulators: {progress.total_emulators_trained}")
advanced_engine.add_progress_callback(progress_callback)
print("Starting interactive advanced workflow...")
print(f"Initial parameter ranges:")
for param_name in advanced_engine.parameter_space.get_parameter_names():
bounds = advanced_engine.parameter_space.get_bounds(param_name)
print(f" {param_name}: [{bounds[0]}, {bounds[1]}]")
Starting interactive advanced workflow... Initial parameter ranges: beta: [0.5, 3.0] gamma: [0.1, 1.0]
In [8]:
Copied!
# Iteration 1: Automatic feature selection
print("\n Iteration 1: Using automatic feature selection...")
result1 = advanced_engine.step()
print(f" Iteration 1 results:")
print(f" Samples generated: {len(result1.samples)}")
print(f" Features selected: {result1.emulated_outputs}")
print(f" Parameter ranges after filtering:")
for param in ['beta', 'gamma']:
values = result1.samples[param]
print(f" {param}: [{values.min():.3f}, {values.max():.3f}]")
# Check emulator quality (in real workflow, you'd inspect diagnostics)
print(f" Emulators trained: {len(result1.emulators)}")
# Accept iteration 1
print(" Accepting iteration 1...")
advanced_engine.commit_step()
# Iteration 1: Automatic feature selection
print("\n Iteration 1: Using automatic feature selection...")
result1 = advanced_engine.step()
print(f" Iteration 1 results:")
print(f" Samples generated: {len(result1.samples)}")
print(f" Features selected: {result1.emulated_outputs}")
print(f" Parameter ranges after filtering:")
for param in ['beta', 'gamma']:
values = result1.samples[param]
print(f" {param}: [{values.min():.3f}, {values.max():.3f}]")
# Check emulator quality (in real workflow, you'd inspect diagnostics)
print(f" Emulators trained: {len(result1.emulators)}")
# Accept iteration 1
print(" Accepting iteration 1...")
advanced_engine.commit_step()
Iteration 1: Using automatic feature selection...
I0000 00:00:1780533791.183240 1037745 gpu_device.cc:2043] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 2175 MB memory: -> device: 0, name: NVIDIA GeForce RTX 3050 Laptop GPU, pci bus id: 0000:01:00.0, compute capability: 8.6
I0000 00:00:1780533795.814762 1037944 cuda_solvers.cc:175] Creating GpuSolver handles for stream 0x6215b44b69a0
Iteration 1 results:
Samples generated: 800
Features selected: ['incidence_day_4', 'cases_week_1', 'time_to_peak']
Parameter ranges after filtering:
beta: [0.502, 2.999]
gamma: [0.101, 0.999]
Emulators trained: 3
Accepting iteration 1...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204310/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204310/wave1
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_emulators_trained'
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204310
INFO:historymatching.engine:============================================================
In [9]:
Copied!
# Iteration 2: Switch to manual feature selection
print("\n Iteration 2: Switching to manual feature selection...")
# Update strategy to focus on key summary statistics
advanced_engine.feature_selection = ['peak_incidence', 'attack_rate']
print(f" Updated to manual feature selection: ['peak_incidence', 'attack_rate']")
result2 = advanced_engine.step()
print(f" Iteration 2 results:")
print(f" Samples generated: {len(result2.samples)}")
print(f" Acceptance rate: {advanced_engine.acceptance_rate:.3f}")
print(f" Features selected: {result2.emulated_outputs}")
print(f" Parameter ranges:")
for param in ['beta', 'gamma']:
values = result2.samples[param]
print(f" {param}: [{values.min():.3f}, {values.max():.3f}]")
print(" Accepting iteration 2...")
advanced_engine.commit_step()
# Iteration 2: Switch to manual feature selection
print("\n Iteration 2: Switching to manual feature selection...")
# Update strategy to focus on key summary statistics
advanced_engine.feature_selection = ['peak_incidence', 'attack_rate']
print(f" Updated to manual feature selection: ['peak_incidence', 'attack_rate']")
result2 = advanced_engine.step()
print(f" Iteration 2 results:")
print(f" Samples generated: {len(result2.samples)}")
print(f" Acceptance rate: {advanced_engine.acceptance_rate:.3f}")
print(f" Features selected: {result2.emulated_outputs}")
print(f" Parameter ranges:")
for param in ['beta', 'gamma']:
values = result2.samples[param]
print(f" {param}: [{values.min():.3f}, {values.max():.3f}]")
print(" Accepting iteration 2...")
advanced_engine.commit_step()
INFO:historymatching.engine:Feature selection: Manual Selection (2 features)
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 800 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 800 simulations...
Iteration 2: Switching to manual feature selection... Updated to manual feature selection: ['peak_incidence', 'attack_rate']
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.6s]
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['peak_incidence', 'attack_rate'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 2 emulators...
INFO:historymatching.emulators.factory: Emulator 1/2 [peak_incidence]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/2 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.1891, lengthscale range=[2.595, 3.234] [2.2s]
INFO:historymatching.emulators.factory: Emulator 1/2 [peak_incidence]: trained [2.2s]
INFO:historymatching.emulators.factory: Emulator 2/2 [attack_rate]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: attack_rate
INFO:historymatching.emulators.factory: Emulator 2/2 [attack_rate]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0077, lengthscale range=[0.743, 1.435] [1.5s]
INFO:historymatching.emulators.factory: Emulator 2/2 [attack_rate]: trained [1.5s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [3.7s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 800 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/800 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 18/800 (2%) | 880 tested | rate=2.0455% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 723/800 (90%) | 42,934 tested | rate=1.6840% [2s]
INFO:historymatching.nroy_sampling: LHS rejection: 809/800 (101%) | 47,963 tested | rate=1.6867% [3s]
INFO:historymatching.nroy_sampling: LHS rejection: 800/800 [2.8s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 800/800 [OK] — sources: LHS=800 [2.8s]
DEBUG:historymatching.engine:NROY fraction: 0.016680
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 800 candidates [2.8s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [7.1s total]. Committing...
Iteration 2 results:
Samples generated: 800
Acceptance rate: 0.034
Features selected: ['peak_incidence', 'attack_rate']
Parameter ranges:
beta: [0.904, 1.915]
gamma: [0.100, 1.000]
Accepting iteration 2...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204310/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204310/wave2
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_emulators_trained'
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204310
INFO:historymatching.engine:============================================================
In [10]:
Copied!
# Iteration 3: Switch emulator type and continue
print("\n Iteration 3: Switching to linear emulator for comparison...")
# Switch to a simpler emulator to see the difference
advanced_engine.emulator_type = 'linear'
print(f" Updated emulator type to: linear")
# Use time series features this time
advanced_engine.feature_selection = ['incidence_day_4', 'incidence_day_8']
print(f" Updated to time series features: ['incidence_day_4', 'incidence_day_8']")
result3 = advanced_engine.step()
print(f" Iteration 3 results:")
print(f" Samples generated: {len(result3.samples)}")
print(f" Acceptance rate: {advanced_engine.acceptance_rate:.3f}")
print(f" Features selected: {result3.emulated_outputs}")
# Check if linear emulator performed poorly
if advanced_engine.acceptance_rate < 0.1:
print(f" Low acceptance rate with linear emulator - reverting to GPR")
advanced_engine.revert_step() # Revert this iteration
# Switch back to GPR and retry
advanced_engine.emulator_type = 'gpr'
result3 = advanced_engine.step()
print(f" Retried with GPR: {len(result3.samples)} samples, rate: {advanced_engine.acceptance_rate:.3f}")
print(" Accepting iteration 3...")
advanced_engine.commit_step()
print(f"\n Interactive workflow completed!")
print(f" Total iterations: {advanced_engine.current_iteration}")
print(f" Final acceptance rate: {advanced_engine.acceptance_rate:.3f}")
# Iteration 3: Switch emulator type and continue
print("\n Iteration 3: Switching to linear emulator for comparison...")
# Switch to a simpler emulator to see the difference
advanced_engine.emulator_type = 'linear'
print(f" Updated emulator type to: linear")
# Use time series features this time
advanced_engine.feature_selection = ['incidence_day_4', 'incidence_day_8']
print(f" Updated to time series features: ['incidence_day_4', 'incidence_day_8']")
result3 = advanced_engine.step()
print(f" Iteration 3 results:")
print(f" Samples generated: {len(result3.samples)}")
print(f" Acceptance rate: {advanced_engine.acceptance_rate:.3f}")
print(f" Features selected: {result3.emulated_outputs}")
# Check if linear emulator performed poorly
if advanced_engine.acceptance_rate < 0.1:
print(f" Low acceptance rate with linear emulator - reverting to GPR")
advanced_engine.revert_step() # Revert this iteration
# Switch back to GPR and retry
advanced_engine.emulator_type = 'gpr'
result3 = advanced_engine.step()
print(f" Retried with GPR: {len(result3.samples)} samples, rate: {advanced_engine.acceptance_rate:.3f}")
print(" Accepting iteration 3...")
advanced_engine.commit_step()
print(f"\n Interactive workflow completed!")
print(f" Total iterations: {advanced_engine.current_iteration}")
print(f" Final acceptance rate: {advanced_engine.acceptance_rate:.3f}")
INFO:historymatching.engine:Emulator type: linear
INFO:historymatching.engine:Feature selection: Manual Selection (2 features)
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 3 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 3] Phase 1/5 SAMPLING: using 800 pre-computed samples from wave 2 [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: running 800 simulations...
Iteration 3: Switching to linear emulator for comparison... Updated emulator type to: linear Updated to time series features: ['incidence_day_4', 'incidence_day_8']
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: complete — 31 outputs [0.7s]
INFO:historymatching.engine:[Wave 3] Phase 3/5 FEATURE SELECTION: ['incidence_day_4', 'incidence_day_8'] [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: training 2 emulators...
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: creating (800 samples, 2 params, type=linear)...
DEBUG:historymatching.emulators.factory:Created linear emulator for feature: incidence_day_4
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: training...
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: trained [0.0s]
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: creating (800 samples, 2 params, type=linear)...
DEBUG:historymatching.emulators.factory:Created linear emulator for feature: incidence_day_8
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: training...
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: trained [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: complete [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: finding 800 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/800 (0%) — starting
WARNING:historymatching.nroy_sampling:Filter failed for 'incidence_day_8': operands could not be broadcast together with shapes (16,) (30,) (16,)
INFO:historymatching.nroy_sampling: LHS rejection: 8/800 (1%) | 880 tested | rate=0.9091% [1s]
WARNING:historymatching.nroy_sampling:Filter failed for 'incidence_day_8': operands could not be broadcast together with shapes (2271,) (4491,) (2271,)
INFO:historymatching.nroy_sampling: LHS rejection: 734/800 (91%) | 96,712 tested | rate=0.7590% [2s]
WARNING:historymatching.nroy_sampling:Filter failed for 'incidence_day_8': operands could not be broadcast together with shapes (241,) (478,) (241,)
INFO:historymatching.nroy_sampling: LHS rejection: 815/800 (101%) | 106,277 tested | rate=0.7669% [2s]
INFO:historymatching.nroy_sampling: LHS rejection: 800/800 [2.2s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 800/800 [OK] — sources: LHS=800 [2.2s]
DEBUG:historymatching.engine:NROY fraction: 0.007527
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: found 800 candidates [2.2s]
INFO:historymatching.engine:[Wave 3] ALL PHASES COMPLETE [2.9s total]. Committing...
INFO:historymatching.engine:Iteration 3 reverted
INFO:historymatching.engine:Emulator type: gpr
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 3 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 3] Phase 1/5 SAMPLING: using 800 pre-computed samples from wave 2 [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: running 800 simulations...
Iteration 3 results: Samples generated: 800 Acceptance rate: 0.023 Features selected: ['incidence_day_4', 'incidence_day_8'] Low acceptance rate with linear emulator - reverting to GPR
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: complete — 31 outputs [0.7s]
INFO:historymatching.engine:[Wave 3] Phase 3/5 FEATURE SELECTION: ['incidence_day_4', 'incidence_day_8'] [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: training 2 emulators...
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: incidence_day_4
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.4889, lengthscale range=[6.081, 7.771] [2.5s]
INFO:historymatching.emulators.factory: Emulator 1/2 [incidence_day_4]: trained [2.5s]
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: incidence_day_8
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.2593, lengthscale range=[3.038, 5.895] [2.0s]
INFO:historymatching.emulators.factory: Emulator 2/2 [incidence_day_8]: trained [2.0s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: complete [4.5s]
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: finding 800 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/800 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 17/800 (2%) | 880 tested | rate=1.9318% [1s]
INFO:historymatching.nroy_sampling: LHS rejection: 771/800 (96%) | 45,464 tested | rate=1.6958% [3s]
INFO:historymatching.nroy_sampling: LHS rejection: 804/800 (100%) | 47,345 tested | rate=1.6982% [4s]
INFO:historymatching.nroy_sampling: LHS rejection: 800/800 [3.7s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 800/800 [OK] — sources: LHS=800 [3.7s]
DEBUG:historymatching.engine:NROY fraction: 0.016897
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: found 800 candidates [3.7s]
INFO:historymatching.engine:[Wave 3] ALL PHASES COMPLETE [8.9s total]. Committing...
Retried with GPR: 800 samples, rate: 0.027 Accepting iteration 3...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204310/checkpoint.pkl
INFO:historymatching.engine:Wave 3 output saved to hm_output/run_20260603_204310/wave3
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_emulators_trained'
INFO:historymatching.engine:[Wave 3] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204310
INFO:historymatching.engine:============================================================
Interactive workflow completed! Total iterations: 3 Final acceptance rate: 0.027
Advanced Analysis and Diagnostics¶
In [11]:
Copied!
# Get all results for analysis
all_results = advanced_engine.get_all_results()
final_samples = advanced_engine.get_nroy_samples()
print(f" Advanced Analysis Results:")
print(f" Total iterations completed: {len(all_results)}")
print(f" Final plausible samples: {len(final_samples)}")
print(f" Total emulators trained: {advanced_engine.emulators_trained}")
print(f" Total samples evaluated: {advanced_engine.samples_generated}")
# Parameter estimation summary
print(f"\n Parameter Recovery Assessment:")
for param, true_value in [('beta', beta_true), ('gamma', gamma_true)]:
values = final_samples[param]
median_est = values.median()
ci_lower = values.quantile(0.025)
ci_upper = values.quantile(0.975)
in_ci = ci_lower <= true_value <= ci_upper
relative_error = abs(median_est - true_value) / true_value
print(f" {param}: Estimated {median_est:.3f} (95% CI: [{ci_lower:.3f}, {ci_upper:.3f}])")
print(f" True value: {true_value} {'' if in_ci else '❌'} ({'in' if in_ci else 'out of'} CI)")
print(f" Relative error: {relative_error:.1%}")
# Check R0 recovery
R0_samples = final_samples['beta'] / final_samples['gamma']
R0_true = beta_true / gamma_true
R0_median = R0_samples.median()
R0_ci_lower = R0_samples.quantile(0.025)
R0_ci_upper = R0_samples.quantile(0.975)
R0_in_ci = R0_ci_lower <= R0_true <= R0_ci_upper
print(f"\n Derived Parameter (R₀):")
print(f" Estimated: {R0_median:.3f} (95% CI: [{R0_ci_lower:.3f}, {R0_ci_upper:.3f}])")
print(f" True value: {R0_true:.3f} {'' if R0_in_ci else '❌'}")
# Get all results for analysis
all_results = advanced_engine.get_all_results()
final_samples = advanced_engine.get_nroy_samples()
print(f" Advanced Analysis Results:")
print(f" Total iterations completed: {len(all_results)}")
print(f" Final plausible samples: {len(final_samples)}")
print(f" Total emulators trained: {advanced_engine.emulators_trained}")
print(f" Total samples evaluated: {advanced_engine.samples_generated}")
# Parameter estimation summary
print(f"\n Parameter Recovery Assessment:")
for param, true_value in [('beta', beta_true), ('gamma', gamma_true)]:
values = final_samples[param]
median_est = values.median()
ci_lower = values.quantile(0.025)
ci_upper = values.quantile(0.975)
in_ci = ci_lower <= true_value <= ci_upper
relative_error = abs(median_est - true_value) / true_value
print(f" {param}: Estimated {median_est:.3f} (95% CI: [{ci_lower:.3f}, {ci_upper:.3f}])")
print(f" True value: {true_value} {'' if in_ci else '❌'} ({'in' if in_ci else 'out of'} CI)")
print(f" Relative error: {relative_error:.1%}")
# Check R0 recovery
R0_samples = final_samples['beta'] / final_samples['gamma']
R0_true = beta_true / gamma_true
R0_median = R0_samples.median()
R0_ci_lower = R0_samples.quantile(0.025)
R0_ci_upper = R0_samples.quantile(0.975)
R0_in_ci = R0_ci_lower <= R0_true <= R0_ci_upper
print(f"\n Derived Parameter (R₀):")
print(f" Estimated: {R0_median:.3f} (95% CI: [{R0_ci_lower:.3f}, {R0_ci_upper:.3f}])")
print(f" True value: {R0_true:.3f} {'' if R0_in_ci else '❌'}")
Advanced Analysis Results:
Total iterations completed: 3
Final plausible samples: 800
Total emulators trained: 7
Total samples evaluated: 116949
Parameter Recovery Assessment:
beta: Estimated 1.299 (95% CI: [1.109, 1.528])
True value: 1.3 (in CI)
Relative error: 0.0%
gamma: Estimated 0.482 (95% CI: [0.295, 0.665])
True value: 0.5 (in CI)
Relative error: 3.6%
Derived Parameter (R₀):
Estimated: 2.699 (95% CI: [2.277, 3.785])
True value: 2.600
In [12]:
Copied!
# Advanced visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# Plot 1: Parameter evolution across iterations
ax = axes[0, 0]
colors = ['lightblue', 'orange', 'lightgreen', 'red']
for i, result in enumerate(all_results):
if i < len(colors):
ax.scatter(result.samples['beta'], result.samples['gamma'],
alpha=0.6, s=20, color=colors[i], label=f'Iteration {i+1}')
ax.axvline(beta_true, color='red', linestyle='--', linewidth=2, alpha=0.8)
ax.axhline(gamma_true, color='red', linestyle='--', linewidth=2, alpha=0.8)
ax.set_xlabel('β (transmission rate)')
ax.set_ylabel('γ (recovery rate)')
ax.set_title('Parameter Space Evolution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 2: Final parameter distributions
ax = axes[0, 1]
ax.hist(final_samples['beta'], bins=25, alpha=0.7, color='blue', density=True, label='β')
ax.axvline(beta_true, color='red', linestyle='--', linewidth=2, label=f'True β = {beta_true}')
ax.axvline(final_samples['beta'].median(), color='blue', linestyle='-', linewidth=2,
label=f'Est. β = {final_samples["beta"].median():.2f}')
ax.set_xlabel('β')
ax.set_ylabel('Density')
ax.set_title('β Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
ax = axes[0, 2]
ax.hist(final_samples['gamma'], bins=25, alpha=0.7, color='orange', density=True, label='γ')
ax.axvline(gamma_true, color='red', linestyle='--', linewidth=2, label=f'True γ = {gamma_true}')
ax.axvline(final_samples['gamma'].median(), color='orange', linestyle='-', linewidth=2,
label=f'Est. γ = {final_samples["gamma"].median():.2f}')
ax.set_xlabel('γ')
ax.set_ylabel('Density')
ax.set_title('γ Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 3: R0 distribution
ax = axes[1, 0]
ax.hist(R0_samples, bins=25, alpha=0.7, color='purple', density=True)
ax.axvline(R0_true, color='red', linestyle='--', linewidth=2, label=f'True R₀ = {R0_true:.2f}')
ax.axvline(R0_median, color='purple', linestyle='-', linewidth=2, label=f'Est. R₀ = {R0_median:.2f}')
ax.set_xlabel('R₀')
ax.set_ylabel('Density')
ax.set_title('R₀ Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 4: Acceptance rate evolution
ax = axes[1, 1]
iterations = range(1, len(all_results) + 1)
sample_counts = [len(result.samples) for result in all_results]
ax.bar(iterations, sample_counts, alpha=0.7, color='green')
ax.set_xlabel('Iteration')
ax.set_ylabel('Plausible Samples')
ax.set_title('Sample Count per Iteration')
ax.grid(True, alpha=0.3)
# Plot 5: Feature selection summary
ax = axes[1, 2]
all_features = []
for result in all_results:
all_features.extend(result.emulated_outputs)
feature_counts = pd.Series(all_features).value_counts()
feature_counts.plot(kind='bar', ax=ax, color='teal', alpha=0.7)
ax.set_xlabel('Features')
ax.set_ylabel('Times Selected')
ax.set_title('Feature Selection Frequency')
ax.tick_params(axis='x', rotation=45)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n Feature Selection Summary:")
for feature, count in feature_counts.head().items():
print(f" {feature}: used {count} time(s)")
# Advanced visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# Plot 1: Parameter evolution across iterations
ax = axes[0, 0]
colors = ['lightblue', 'orange', 'lightgreen', 'red']
for i, result in enumerate(all_results):
if i < len(colors):
ax.scatter(result.samples['beta'], result.samples['gamma'],
alpha=0.6, s=20, color=colors[i], label=f'Iteration {i+1}')
ax.axvline(beta_true, color='red', linestyle='--', linewidth=2, alpha=0.8)
ax.axhline(gamma_true, color='red', linestyle='--', linewidth=2, alpha=0.8)
ax.set_xlabel('β (transmission rate)')
ax.set_ylabel('γ (recovery rate)')
ax.set_title('Parameter Space Evolution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 2: Final parameter distributions
ax = axes[0, 1]
ax.hist(final_samples['beta'], bins=25, alpha=0.7, color='blue', density=True, label='β')
ax.axvline(beta_true, color='red', linestyle='--', linewidth=2, label=f'True β = {beta_true}')
ax.axvline(final_samples['beta'].median(), color='blue', linestyle='-', linewidth=2,
label=f'Est. β = {final_samples["beta"].median():.2f}')
ax.set_xlabel('β')
ax.set_ylabel('Density')
ax.set_title('β Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
ax = axes[0, 2]
ax.hist(final_samples['gamma'], bins=25, alpha=0.7, color='orange', density=True, label='γ')
ax.axvline(gamma_true, color='red', linestyle='--', linewidth=2, label=f'True γ = {gamma_true}')
ax.axvline(final_samples['gamma'].median(), color='orange', linestyle='-', linewidth=2,
label=f'Est. γ = {final_samples["gamma"].median():.2f}')
ax.set_xlabel('γ')
ax.set_ylabel('Density')
ax.set_title('γ Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 3: R0 distribution
ax = axes[1, 0]
ax.hist(R0_samples, bins=25, alpha=0.7, color='purple', density=True)
ax.axvline(R0_true, color='red', linestyle='--', linewidth=2, label=f'True R₀ = {R0_true:.2f}')
ax.axvline(R0_median, color='purple', linestyle='-', linewidth=2, label=f'Est. R₀ = {R0_median:.2f}')
ax.set_xlabel('R₀')
ax.set_ylabel('Density')
ax.set_title('R₀ Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot 4: Acceptance rate evolution
ax = axes[1, 1]
iterations = range(1, len(all_results) + 1)
sample_counts = [len(result.samples) for result in all_results]
ax.bar(iterations, sample_counts, alpha=0.7, color='green')
ax.set_xlabel('Iteration')
ax.set_ylabel('Plausible Samples')
ax.set_title('Sample Count per Iteration')
ax.grid(True, alpha=0.3)
# Plot 5: Feature selection summary
ax = axes[1, 2]
all_features = []
for result in all_results:
all_features.extend(result.emulated_outputs)
feature_counts = pd.Series(all_features).value_counts()
feature_counts.plot(kind='bar', ax=ax, color='teal', alpha=0.7)
ax.set_xlabel('Features')
ax.set_ylabel('Times Selected')
ax.set_title('Feature Selection Frequency')
ax.tick_params(axis='x', rotation=45)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n Feature Selection Summary:")
for feature, count in feature_counts.head().items():
print(f" {feature}: used {count} time(s)")
Feature Selection Summary: incidence_day_4: used 2 time(s) cases_week_1: used 1 time(s) time_to_peak: used 1 time(s) peak_incidence: used 1 time(s) attack_rate: used 1 time(s)
Performance Comparison: Different Configurations¶
Let's compare different configuration strategies to understand their impact.
In [13]:
Copied!
def run_quick_comparison(config_name, sampling_strategy='lhs', emulator_type='gpr',
implausibility_threshold=3.0):
"""Run a quick 2-iteration comparison with different settings."""
import time
print(f"\nTesting {config_name}...")
simple_obs = {
'peak_incidence': (max(incidence_obs), 50),
'total_cases': (sum(incidence_obs), 200),
}
engine = hm.HistoryMatching(
bounds={'beta': (0.5, 3.0), 'gamma': (0.1, 1.0)},
observations=simple_obs,
function=enhanced_sir_simulation,
sampling_strategy=sampling_strategy,
emulator_type=emulator_type,
n_samples=200,
max_iterations=2,
implausibility_threshold=implausibility_threshold,
random_seed=999,
)
start_time = time.time()
results = engine.run()
end_time = time.time()
final_samples = engine.get_nroy_samples()
beta_error = abs(final_samples['beta'].median() - beta_true) / beta_true
gamma_error = abs(final_samples['gamma'].median() - gamma_true) / gamma_true
print(f" Runtime: {end_time - start_time:.1f}s")
print(f" Final samples: {len(final_samples)}")
print(f" beta error: {beta_error:.1%}, gamma error: {gamma_error:.1%}")
print(f" Final acceptance rate: {engine.acceptance_rate:.3f}")
return {
'config': config_name,
'runtime': end_time - start_time,
'samples': len(final_samples),
'beta_error': beta_error,
'gamma_error': gamma_error,
'acceptance_rate': engine.acceptance_rate,
}
# Compare configurations
comparisons = []
comparisons.append(run_quick_comparison("LHS + GPR (baseline)", sampling_strategy='lhs', emulator_type='gpr'))
comparisons.append(run_quick_comparison("LHS + Linear", sampling_strategy='lhs', emulator_type='linear'))
comparisons.append(run_quick_comparison("Random + GPR", sampling_strategy='random', emulator_type='gpr'))
comparisons.append(run_quick_comparison("LHS + GPR (conservative)", sampling_strategy='lhs', emulator_type='gpr',
implausibility_threshold=2.5))
def run_quick_comparison(config_name, sampling_strategy='lhs', emulator_type='gpr',
implausibility_threshold=3.0):
"""Run a quick 2-iteration comparison with different settings."""
import time
print(f"\nTesting {config_name}...")
simple_obs = {
'peak_incidence': (max(incidence_obs), 50),
'total_cases': (sum(incidence_obs), 200),
}
engine = hm.HistoryMatching(
bounds={'beta': (0.5, 3.0), 'gamma': (0.1, 1.0)},
observations=simple_obs,
function=enhanced_sir_simulation,
sampling_strategy=sampling_strategy,
emulator_type=emulator_type,
n_samples=200,
max_iterations=2,
implausibility_threshold=implausibility_threshold,
random_seed=999,
)
start_time = time.time()
results = engine.run()
end_time = time.time()
final_samples = engine.get_nroy_samples()
beta_error = abs(final_samples['beta'].median() - beta_true) / beta_true
gamma_error = abs(final_samples['gamma'].median() - gamma_true) / gamma_true
print(f" Runtime: {end_time - start_time:.1f}s")
print(f" Final samples: {len(final_samples)}")
print(f" beta error: {beta_error:.1%}, gamma error: {gamma_error:.1%}")
print(f" Final acceptance rate: {engine.acceptance_rate:.3f}")
return {
'config': config_name,
'runtime': end_time - start_time,
'samples': len(final_samples),
'beta_error': beta_error,
'gamma_error': gamma_error,
'acceptance_rate': engine.acceptance_rate,
}
# Compare configurations
comparisons = []
comparisons.append(run_quick_comparison("LHS + GPR (baseline)", sampling_strategy='lhs', emulator_type='gpr'))
comparisons.append(run_quick_comparison("LHS + Linear", sampling_strategy='lhs', emulator_type='linear'))
comparisons.append(run_quick_comparison("Random + GPR", sampling_strategy='random', emulator_type='gpr'))
comparisons.append(run_quick_comparison("LHS + GPR (conservative)", sampling_strategy='lhs', emulator_type='gpr',
implausibility_threshold=2.5))
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:HISTORY MATCHING — CONFIGURATION
INFO:historymatching.engine:============================================================
INFO:historymatching.engine: Emulator type: gpr
INFO:historymatching.engine: Parameters: 2
INFO:historymatching.engine: Observation targets: 2
INFO:historymatching.engine: Samples per wave: 200
INFO:historymatching.engine: Max iterations: 2
INFO:historymatching.engine: Implausibility threshold: 3.0
INFO:historymatching.engine: Auto space reduction: disabled
INFO:historymatching.engine: Random seed: 999
INFO:historymatching.engine: Oversample factor: 1.1
INFO:historymatching.engine: Max batch size: 10000
INFO:historymatching.engine: Output directory: hm_output/run_20260603_204357
INFO:historymatching.engine: Run log: hm_output/run_20260603_204357/log.txt
INFO:historymatching.engine: Parameters: ['beta', 'gamma']
INFO:historymatching.engine: Targets: ['peak_incidence', 'total_cases']
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Starting automated run with 2 max iterations
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 1 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 1] Phase 1/5 SAMPLING: generated 200 samples (acceptance rate: 1.000) [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: running 200 simulations...
Testing LHS + GPR (baseline)...
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=1171.38, total_cases=207.49
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=1171.38, mean_z=16.11, std_z=30.27)
INFO:historymatching.engine:[Wave 1] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0098, lengthscale range=[0.562, 1.279] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 33/200 (16%) | 220 tested | rate=15.0000% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 234/200 (117%) | 1,444 tested | rate=16.2050% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.2s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.2s]
DEBUG:historymatching.engine:NROY fraction: 0.138504
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: found 200 candidates [0.2s]
INFO:historymatching.engine:[Wave 1] ALL PHASES COMPLETE [0.9s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204357/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204357/wave1
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204357
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 200 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=30.10, total_cases=15.73
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=15.73, mean_z=-1.55, std_z=3.66)
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0094, lengthscale range=[0.815, 1.007] [0.7s]
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [0.7s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [0.7s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 14/200 (7%) | 220 tested | rate=6.3636% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 312/200 (156%) | 3,435 tested | rate=9.0830% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.3s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.3s]
DEBUG:historymatching.engine:NROY fraction: 0.058224
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 200 candidates [0.3s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [1.2s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204357/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204357/wave2
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204357
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Automated run completed. 2 iterations executed.
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:HISTORY MATCHING — CONFIGURATION
INFO:historymatching.engine:============================================================
INFO:historymatching.engine: Emulator type: linear
INFO:historymatching.engine: Parameters: 2
INFO:historymatching.engine: Observation targets: 2
INFO:historymatching.engine: Samples per wave: 200
INFO:historymatching.engine: Max iterations: 2
INFO:historymatching.engine: Implausibility threshold: 3.0
INFO:historymatching.engine: Auto space reduction: disabled
INFO:historymatching.engine: Random seed: 999
INFO:historymatching.engine: Oversample factor: 1.1
INFO:historymatching.engine: Max batch size: 10000
INFO:historymatching.engine: Output directory: hm_output/run_20260603_204405
INFO:historymatching.engine: Run log: hm_output/run_20260603_204405/log.txt
INFO:historymatching.engine: Parameters: ['beta', 'gamma']
INFO:historymatching.engine: Targets: ['peak_incidence', 'total_cases']
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Starting automated run with 2 max iterations
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 1 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 1] Phase 1/5 SAMPLING: generated 200 samples (acceptance rate: 1.000) [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=1179.60, total_cases=208.41
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=1179.60, mean_z=16.22, std_z=30.35)
INFO:historymatching.engine:[Wave 1] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (200 samples, 2 params, type=linear)...
DEBUG:historymatching.emulators.factory:Created linear emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: complete [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
Runtime: 7.3s Final samples: 200 beta error: 5.4%, gamma error: 16.2% Final acceptance rate: 0.118 Testing LHS + Linear...
INFO:historymatching.nroy_sampling: LHS rejection: 15/200 (7%) | 220 tested | rate=6.8182% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 252/200 (126%) | 3,204 tested | rate=7.8652% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.0s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.0s]
DEBUG:historymatching.engine:NROY fraction: 0.062422
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: found 200 candidates [0.0s]
INFO:historymatching.engine:[Wave 1] ALL PHASES COMPLETE [0.2s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204405/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204405/wave1
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204405
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 200 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): total_cases=18.53, peak_incidence=11.04
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=18.53, mean_z=-1.75, std_z=3.94)
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (200 samples, 2 params, type=linear)...
DEBUG:historymatching.emulators.factory:Created linear emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
WARNING:historymatching.nroy_sampling:Filter failed for 'peak_incidence': operands could not be broadcast together with shapes (40,) (72,) (40,)
INFO:historymatching.nroy_sampling: LHS rejection: 40/200 (20%) | 220 tested | rate=18.1818% [0s]
WARNING:historymatching.nroy_sampling:Filter failed for 'peak_incidence': operands could not be broadcast together with shapes (212,) (377,) (212,)
INFO:historymatching.nroy_sampling: LHS rejection: 252/200 (126%) | 1,188 tested | rate=21.2121% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.0s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.0s]
DEBUG:historymatching.engine:NROY fraction: 0.168350
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 200 candidates [0.0s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [0.2s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204405/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204405/wave2
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204405
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Automated run completed. 2 iterations executed.
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:HISTORY MATCHING — CONFIGURATION
INFO:historymatching.engine:============================================================
INFO:historymatching.engine: Emulator type: gpr
INFO:historymatching.engine: Parameters: 2
INFO:historymatching.engine: Observation targets: 2
INFO:historymatching.engine: Samples per wave: 200
INFO:historymatching.engine: Max iterations: 2
INFO:historymatching.engine: Implausibility threshold: 3.0
INFO:historymatching.engine: Auto space reduction: disabled
INFO:historymatching.engine: Random seed: 999
INFO:historymatching.engine: Oversample factor: 1.1
INFO:historymatching.engine: Max batch size: 10000
INFO:historymatching.engine: Output directory: hm_output/run_20260603_204410
INFO:historymatching.engine: Run log: hm_output/run_20260603_204410/log.txt
INFO:historymatching.engine: Parameters: ['beta', 'gamma']
INFO:historymatching.engine: Targets: ['peak_incidence', 'total_cases']
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Starting automated run with 2 max iterations
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 1 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 1] Phase 1/5 SAMPLING: generated 200 samples (acceptance rate: 1.000) [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: running 200 simulations...
Runtime: 5.0s Final samples: 200 beta error: 21.4%, gamma error: 15.4% Final acceptance rate: 0.131 Testing Random + GPR...
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=1154.06, total_cases=194.54
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=1154.06, mean_z=17.65, std_z=29.10)
INFO:historymatching.engine:[Wave 1] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0128, lengthscale range=[0.621, 1.037] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 36/200 (18%) | 220 tested | rate=16.3636% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 261/200 (130%) | 1,322 tested | rate=19.7428% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.2s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.2s]
DEBUG:historymatching.engine:NROY fraction: 0.151286
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: found 200 candidates [0.2s]
INFO:historymatching.engine:[Wave 1] ALL PHASES COMPLETE [0.9s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204410/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204410/wave1
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204410
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 200 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=40.09, total_cases=16.54
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=16.54, mean_z=-1.26, std_z=3.88)
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0078, lengthscale range=[0.737, 0.818] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [0.6s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 20/200 (10%) | 220 tested | rate=9.0909% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 254/200 (127%) | 2,398 tested | rate=10.5922% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.3s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.3s]
DEBUG:historymatching.engine:NROY fraction: 0.083403
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 200 candidates [0.4s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [1.2s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204410/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204410/wave2
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204410
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Automated run completed. 2 iterations executed.
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:HISTORY MATCHING — CONFIGURATION
INFO:historymatching.engine:============================================================
INFO:historymatching.engine: Emulator type: gpr
INFO:historymatching.engine: Parameters: 2
INFO:historymatching.engine: Observation targets: 2
INFO:historymatching.engine: Samples per wave: 200
INFO:historymatching.engine: Max iterations: 2
INFO:historymatching.engine: Implausibility threshold: 2.5
INFO:historymatching.engine: Auto space reduction: disabled
INFO:historymatching.engine: Random seed: 999
INFO:historymatching.engine: Oversample factor: 1.1
INFO:historymatching.engine: Max batch size: 10000
INFO:historymatching.engine: Output directory: hm_output/run_20260603_204417
INFO:historymatching.engine: Run log: hm_output/run_20260603_204417/log.txt
INFO:historymatching.engine: Parameters: ['beta', 'gamma']
INFO:historymatching.engine: Targets: ['peak_incidence', 'total_cases']
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Starting automated run with 2 max iterations
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 1 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 1] Phase 1/5 SAMPLING: generated 200 samples (acceptance rate: 1.000) [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
Runtime: 7.0s Final samples: 200 beta error: 8.6%, gamma error: 22.3% Final acceptance rate: 0.153 Testing LHS + GPR (conservative)...
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=1170.99, total_cases=208.49
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=1170.99, mean_z=16.09, std_z=30.28)
INFO:historymatching.engine:[Wave 1] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0094, lengthscale range=[0.644, 0.724] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 30/200 (15%) | 220 tested | rate=13.6364% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 233/200 (116%) | 1,591 tested | rate=14.6449% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.2s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.2s]
DEBUG:historymatching.engine:NROY fraction: 0.125707
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: found 200 candidates [0.2s]
INFO:historymatching.engine:[Wave 1] ALL PHASES COMPLETE [0.9s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204417/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204417/wave1
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204417
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 200 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 200 simulations...
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=24.68, total_cases=13.43
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=13.43, mean_z=-1.00, std_z=3.53)
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (200 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.gpr: Training GPR: 150 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0078, lengthscale range=[1.023, 1.331] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [0.6s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 200 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/200 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 10/200 (5%) | 220 tested | rate=4.5455% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 287/200 (143%) | 4,818 tested | rate=5.9568% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 200/200 [0.4s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 200/200 [OK] — sources: LHS=200 [0.4s]
DEBUG:historymatching.engine:NROY fraction: 0.041511
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 200 candidates [0.5s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [1.3s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204417/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204417/wave2
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204417
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Automated run completed. 2 iterations executed.
Runtime: 7.1s Final samples: 200 beta error: 2.4%, gamma error: 9.0% Final acceptance rate: 0.091
In [14]:
Copied!
# Summarize comparison results
comparison_df = pd.DataFrame(comparisons)
print("\n Configuration Comparison Summary:")
print("=" * 80)
print(f"{'Configuration':<25} {'Runtime(s)':<12} {'Samples':<10} {'β Error':<10} {'γ Error':<10} {'Acc. Rate':<10}")
print("=" * 80)
for _, row in comparison_df.iterrows():
print(f"{row['config']:<25} {row['runtime']:<12.1f} {row['samples']:<10} "
f"{row['beta_error']:<10.1%} {row['gamma_error']:<10.1%} {row['acceptance_rate']:<10.3f}")
print("\n Best performing configurations:")
print(f" Fastest: {comparison_df.loc[comparison_df['runtime'].idxmin(), 'config']}")
print(f" Most samples: {comparison_df.loc[comparison_df['samples'].idxmax(), 'config']}")
print(f" Lowest β error: {comparison_df.loc[comparison_df['beta_error'].idxmin(), 'config']}")
print(f" Lowest γ error: {comparison_df.loc[comparison_df['gamma_error'].idxmin(), 'config']}")
# Plot comparison
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Runtime comparison
ax = axes[0]
comparison_df.plot(x='config', y='runtime', kind='bar', ax=ax, color='skyblue')
ax.set_title('Runtime Comparison')
ax.set_ylabel('Time (seconds)')
ax.tick_params(axis='x', rotation=45)
# Error comparison
ax = axes[1]
x = range(len(comparison_df))
width = 0.35
ax.bar([i - width/2 for i in x], comparison_df['beta_error'], width, label='β error', alpha=0.7)
ax.bar([i + width/2 for i in x], comparison_df['gamma_error'], width, label='γ error', alpha=0.7)
ax.set_xlabel('Configuration')
ax.set_ylabel('Relative Error')
ax.set_title('Parameter Estimation Error')
ax.set_xticks(x)
ax.set_xticklabels([c[:15] + '...' if len(c) > 15 else c for c in comparison_df['config']], rotation=45)
ax.legend()
# Sample count vs acceptance rate
ax = axes[2]
scatter = ax.scatter(comparison_df['acceptance_rate'], comparison_df['samples'],
s=100, alpha=0.7, c=comparison_df['runtime'], cmap='viridis')
ax.set_xlabel('Final Acceptance Rate')
ax.set_ylabel('Final Sample Count')
ax.set_title('Acceptance Rate vs Sample Count')
plt.colorbar(scatter, ax=ax, label='Runtime (s)')
# Add config labels
for i, row in comparison_df.iterrows():
ax.annotate(row['config'][:10] + '...',
(row['acceptance_rate'], row['samples']),
xytext=(5, 5), textcoords='offset points', fontsize=8)
plt.tight_layout()
plt.show()
# Summarize comparison results
comparison_df = pd.DataFrame(comparisons)
print("\n Configuration Comparison Summary:")
print("=" * 80)
print(f"{'Configuration':<25} {'Runtime(s)':<12} {'Samples':<10} {'β Error':<10} {'γ Error':<10} {'Acc. Rate':<10}")
print("=" * 80)
for _, row in comparison_df.iterrows():
print(f"{row['config']:<25} {row['runtime']:<12.1f} {row['samples']:<10} "
f"{row['beta_error']:<10.1%} {row['gamma_error']:<10.1%} {row['acceptance_rate']:<10.3f}")
print("\n Best performing configurations:")
print(f" Fastest: {comparison_df.loc[comparison_df['runtime'].idxmin(), 'config']}")
print(f" Most samples: {comparison_df.loc[comparison_df['samples'].idxmax(), 'config']}")
print(f" Lowest β error: {comparison_df.loc[comparison_df['beta_error'].idxmin(), 'config']}")
print(f" Lowest γ error: {comparison_df.loc[comparison_df['gamma_error'].idxmin(), 'config']}")
# Plot comparison
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Runtime comparison
ax = axes[0]
comparison_df.plot(x='config', y='runtime', kind='bar', ax=ax, color='skyblue')
ax.set_title('Runtime Comparison')
ax.set_ylabel('Time (seconds)')
ax.tick_params(axis='x', rotation=45)
# Error comparison
ax = axes[1]
x = range(len(comparison_df))
width = 0.35
ax.bar([i - width/2 for i in x], comparison_df['beta_error'], width, label='β error', alpha=0.7)
ax.bar([i + width/2 for i in x], comparison_df['gamma_error'], width, label='γ error', alpha=0.7)
ax.set_xlabel('Configuration')
ax.set_ylabel('Relative Error')
ax.set_title('Parameter Estimation Error')
ax.set_xticks(x)
ax.set_xticklabels([c[:15] + '...' if len(c) > 15 else c for c in comparison_df['config']], rotation=45)
ax.legend()
# Sample count vs acceptance rate
ax = axes[2]
scatter = ax.scatter(comparison_df['acceptance_rate'], comparison_df['samples'],
s=100, alpha=0.7, c=comparison_df['runtime'], cmap='viridis')
ax.set_xlabel('Final Acceptance Rate')
ax.set_ylabel('Final Sample Count')
ax.set_title('Acceptance Rate vs Sample Count')
plt.colorbar(scatter, ax=ax, label='Runtime (s)')
# Add config labels
for i, row in comparison_df.iterrows():
ax.annotate(row['config'][:10] + '...',
(row['acceptance_rate'], row['samples']),
xytext=(5, 5), textcoords='offset points', fontsize=8)
plt.tight_layout()
plt.show()
Configuration Comparison Summary: ================================================================================ Configuration Runtime(s) Samples β Error γ Error Acc. Rate ================================================================================ LHS + GPR (baseline) 7.3 200 5.4% 16.2% 0.118 LHS + Linear 5.0 200 21.4% 15.4% 0.131 Random + GPR 7.0 200 8.6% 22.3% 0.153 LHS + GPR (conservative) 7.1 200 2.4% 9.0% 0.091 Best performing configurations: Fastest: LHS + Linear Most samples: LHS + GPR (baseline) Lowest β error: LHS + GPR (conservative) Lowest γ error: LHS + GPR (conservative)
Checkpoint and Resume Functionality¶
For long-running workflows, the engine supports saving and loading checkpoints.
In [15]:
Copied!
# Build a small engine to demonstrate monitoring
simple_obs = {
'peak_incidence': (incidence_obs.max(), 50),
'total_cases': (incidence_obs.sum(), 200),
}
progress_log = []
feature_history = []
def monitor_callback(progress):
progress_log.append({
'iteration': progress.current_iteration,
'acceptance': progress.acceptance_rate,
'total': progress.total_samples_generated,
'accepted': progress.total_samples_accepted,
})
monitor_engine = hm.HistoryMatching(
bounds={'beta': (0.5, 3.0), 'gamma': (0.1, 1.0)},
observations=simple_obs,
function=enhanced_sir_simulation,
emulator_type='gpr',
n_samples=300,
max_iterations=4,
random_seed=7,
)
monitor_engine.add_progress_callback(monitor_callback)
mon_results = monitor_engine.run()
for r in mon_results:
feature_history.append(r.emulated_outputs)
print("Iteration | Acceptance | Features selected")
print("-" * 55)
for i, (p, feats) in enumerate(zip(progress_log, feature_history), 1):
print(f" {i} | {p['acceptance']:.3f} | {feats}")
# Build a small engine to demonstrate monitoring
simple_obs = {
'peak_incidence': (incidence_obs.max(), 50),
'total_cases': (incidence_obs.sum(), 200),
}
progress_log = []
feature_history = []
def monitor_callback(progress):
progress_log.append({
'iteration': progress.current_iteration,
'acceptance': progress.acceptance_rate,
'total': progress.total_samples_generated,
'accepted': progress.total_samples_accepted,
})
monitor_engine = hm.HistoryMatching(
bounds={'beta': (0.5, 3.0), 'gamma': (0.1, 1.0)},
observations=simple_obs,
function=enhanced_sir_simulation,
emulator_type='gpr',
n_samples=300,
max_iterations=4,
random_seed=7,
)
monitor_engine.add_progress_callback(monitor_callback)
mon_results = monitor_engine.run()
for r in mon_results:
feature_history.append(r.emulated_outputs)
print("Iteration | Acceptance | Features selected")
print("-" * 55)
for i, (p, feats) in enumerate(zip(progress_log, feature_history), 1):
print(f" {i} | {p['acceptance']:.3f} | {feats}")
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:HISTORY MATCHING — CONFIGURATION
INFO:historymatching.engine:============================================================
INFO:historymatching.engine: Emulator type: gpr
INFO:historymatching.engine: Parameters: 2
INFO:historymatching.engine: Observation targets: 2
INFO:historymatching.engine: Samples per wave: 300
INFO:historymatching.engine: Max iterations: 4
INFO:historymatching.engine: Implausibility threshold: 3.0
INFO:historymatching.engine: Auto space reduction: disabled
INFO:historymatching.engine: Random seed: 7
INFO:historymatching.engine: Oversample factor: 1.1
INFO:historymatching.engine: Max batch size: 10000
INFO:historymatching.engine: Output directory: hm_output/run_20260603_204426
INFO:historymatching.engine: Run log: hm_output/run_20260603_204426/log.txt
INFO:historymatching.engine: Parameters: ['beta', 'gamma']
INFO:historymatching.engine: Targets: ['peak_incidence', 'total_cases']
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Starting automated run with 4 max iterations
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 1 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 1] Phase 1/5 SAMPLING: generated 300 samples (acceptance rate: 1.000) [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: running 300 simulations...
INFO:historymatching.engine:[Wave 1] Phase 2/5 SIMULATION: complete — 31 outputs [0.2s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=1154.09, total_cases=216.58
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=1154.09, mean_z=16.26, std_z=29.88)
INFO:historymatching.engine:[Wave 1] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (300 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 225 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0079, lengthscale range=[0.532, 1.172] [0.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 4/5 EMULATOR TRAINING: complete [0.6s]
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: finding 300 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/300 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 43/300 (14%) | 330 tested | rate=13.0303% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 391/300 (130%) | 2,499 tested | rate=15.6463% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 300/300 [0.3s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 300/300 [OK] — sources: LHS=300 [0.3s]
DEBUG:historymatching.engine:NROY fraction: 0.120048
INFO:historymatching.engine:[Wave 1] Phase 5/5 NROY SAMPLING: found 300 candidates [0.3s]
INFO:historymatching.engine:[Wave 1] ALL PHASES COMPLETE [1.1s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204426/checkpoint.pkl
INFO:historymatching.engine:Wave 1 output saved to hm_output/run_20260603_204426/wave1
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_samples_generated'
INFO:historymatching.engine:[Wave 1] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204426
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 2 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 2] Phase 1/5 SAMPLING: using 300 pre-computed samples from wave 1 [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: running 300 simulations...
INFO:historymatching.engine:[Wave 2] Phase 2/5 SIMULATION: complete — 31 outputs [0.3s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=23.88, total_cases=15.24
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=15.24, mean_z=-1.58, std_z=3.58)
INFO:historymatching.engine:[Wave 2] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (300 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.gpr: Training GPR: 225 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0081, lengthscale range=[0.891, 1.087] [0.8s]
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [0.8s]
INFO:historymatching.engine:[Wave 2] Phase 4/5 EMULATOR TRAINING: complete [0.8s]
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: finding 300 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/300 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 24/300 (8%) | 330 tested | rate=7.2727% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 353/300 (117%) | 4,504 tested | rate=7.8375% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 300/300 [0.5s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 300/300 [OK] — sources: LHS=300 [0.5s]
DEBUG:historymatching.engine:NROY fraction: 0.066607
INFO:historymatching.engine:[Wave 2] Phase 5/5 NROY SAMPLING: found 300 candidates [0.5s]
INFO:historymatching.engine:[Wave 2] ALL PHASES COMPLETE [1.6s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204426/checkpoint.pkl
INFO:historymatching.engine:Wave 2 output saved to hm_output/run_20260603_204426/wave2
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_samples_generated'
INFO:historymatching.engine:[Wave 2] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204426
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 3 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 3] Phase 1/5 SAMPLING: using 300 pre-computed samples from wave 2 [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: running 300 simulations...
INFO:historymatching.engine:[Wave 3] Phase 2/5 SIMULATION: complete — 31 outputs [0.5s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=28.54, total_cases=3.23
INFO:historymatching.feature_selection: Cooldown (skip): ['total_cases']
INFO:historymatching.feature_selection: SELECTED: peak_incidence (score=28.54, mean_z=1.08, std_z=5.24)
INFO:historymatching.engine:[Wave 3] Phase 3/5 FEATURE SELECTION: ['peak_incidence'] [0.0s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: creating (300 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: peak_incidence
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: training...
INFO:historymatching.emulators.gpr: Training GPR: 225 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0534, lengthscale range=[9.473, 12.025] [1.6s]
INFO:historymatching.emulators.factory: Emulator 1/1 [peak_incidence]: trained [1.6s]
INFO:historymatching.engine:[Wave 3] Phase 4/5 EMULATOR TRAINING: complete [1.6s]
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: finding 300 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/300 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 12/300 (4%) | 330 tested | rate=3.6364% [0s]
INFO:historymatching.nroy_sampling: LHS rejection: 384/300 (128%) | 9,042 tested | rate=4.2468% [1s]
INFO:historymatching.nroy_sampling: LHS rejection: 300/300 [1.0s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 300/300 [OK] — sources: LHS=300 [1.0s]
DEBUG:historymatching.engine:NROY fraction: 0.033179
INFO:historymatching.engine:[Wave 3] Phase 5/5 NROY SAMPLING: found 300 candidates [1.1s]
INFO:historymatching.engine:[Wave 3] ALL PHASES COMPLETE [3.2s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204426/checkpoint.pkl
INFO:historymatching.engine:Wave 3 output saved to hm_output/run_20260603_204426/wave3
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_samples_generated'
INFO:historymatching.engine:[Wave 3] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204426
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 4 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 4] Phase 1/5 SAMPLING: using 300 pre-computed samples from wave 3 [0.0s]
INFO:historymatching.engine:[Wave 4] Phase 2/5 SIMULATION: running 300 simulations...
INFO:historymatching.engine:[Wave 4] Phase 2/5 SIMULATION: complete — 31 outputs [0.4s]
INFO:historymatching.feature_selection: Feature ranking (mean_sq_z): peak_incidence=8.48, total_cases=3.94
INFO:historymatching.feature_selection: Cooldown (skip): ['peak_incidence']
INFO:historymatching.feature_selection: SELECTED: total_cases (score=3.94, mean_z=0.12, std_z=1.98)
INFO:historymatching.engine:[Wave 4] Phase 3/5 FEATURE SELECTION: ['total_cases'] [0.0s]
INFO:historymatching.engine:[Wave 4] Phase 4/5 EMULATOR TRAINING: training 1 emulators...
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: creating (300 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: total_cases
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: training...
INFO:historymatching.emulators.gpr: Training GPR: 225 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0130, lengthscale range=[1.086, 1.319] [1.0s]
INFO:historymatching.emulators.factory: Emulator 1/1 [total_cases]: trained [1.0s]
INFO:historymatching.engine:[Wave 4] Phase 4/5 EMULATOR TRAINING: complete [1.0s]
INFO:historymatching.engine:[Wave 4] Phase 5/5 NROY SAMPLING: finding 300 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/300 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 12/300 (4%) | 330 tested | rate=3.6364% [1s]
INFO:historymatching.nroy_sampling: LHS rejection: 369/300 (123%) | 9,042 tested | rate=4.0810% [1s]
INFO:historymatching.nroy_sampling: LHS rejection: 300/300 [1.1s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 300/300 [OK] — sources: LHS=300 [1.1s]
DEBUG:historymatching.engine:NROY fraction: 0.033179
INFO:historymatching.engine:[Wave 4] Phase 5/5 NROY SAMPLING: found 300 candidates [1.1s]
INFO:historymatching.engine:[Wave 4] ALL PHASES COMPLETE [2.6s total]. Committing...
INFO:historymatching.engine:Checkpoint saved to hm_output/run_20260603_204426/checkpoint.pkl
INFO:historymatching.engine:Wave 4 output saved to hm_output/run_20260603_204426/wave4
WARNING:historymatching.engine:Progress callback failed: 'HistoryMatching' object has no attribute 'total_samples_generated'
INFO:historymatching.engine:[Wave 4] COMMITTED — diagnostics and checkpoint saved to hm_output/run_20260603_204426
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:Automated run completed. 4 iterations executed.
Iteration | Acceptance | Features selected -------------------------------------------------------
In [16]:
Copied!
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# Acceptance rate over iterations
ax = axes[0]
iters = [p['iteration'] for p in progress_log]
rates = [p['acceptance'] for p in progress_log]
ax.plot(iters, rates, 'bo-', linewidth=2, markersize=8)
ax.set_xlabel("Iteration")
ax.set_ylabel("Acceptance Rate")
ax.set_title("Convergence: Acceptance Rate")
ax.grid(True, alpha=0.3)
# Cumulative samples
ax = axes[1]
ax.plot(iters, [p['total'] for p in progress_log], 'b-', label='Generated', linewidth=2)
ax.plot(iters, [p['accepted'] for p in progress_log], 'g-', label='Accepted', linewidth=2)
ax.set_xlabel("Iteration")
ax.set_ylabel("Cumulative Samples")
ax.set_title("Sample Generation Progress")
ax.legend()
ax.grid(True, alpha=0.3)
# Feature selection frequency
ax = axes[2]
from collections import Counter
feat_counts = Counter(f for feats in feature_history for f in feats)
ax.bar(feat_counts.keys(), feat_counts.values(), color='teal', alpha=0.7)
ax.set_xlabel("Feature")
ax.set_ylabel("Times Selected")
ax.set_title("Feature Selection Frequency")
ax.tick_params(axis='x', rotation=30)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# Acceptance rate over iterations
ax = axes[0]
iters = [p['iteration'] for p in progress_log]
rates = [p['acceptance'] for p in progress_log]
ax.plot(iters, rates, 'bo-', linewidth=2, markersize=8)
ax.set_xlabel("Iteration")
ax.set_ylabel("Acceptance Rate")
ax.set_title("Convergence: Acceptance Rate")
ax.grid(True, alpha=0.3)
# Cumulative samples
ax = axes[1]
ax.plot(iters, [p['total'] for p in progress_log], 'b-', label='Generated', linewidth=2)
ax.plot(iters, [p['accepted'] for p in progress_log], 'g-', label='Accepted', linewidth=2)
ax.set_xlabel("Iteration")
ax.set_ylabel("Cumulative Samples")
ax.set_title("Sample Generation Progress")
ax.legend()
ax.grid(True, alpha=0.3)
# Feature selection frequency
ax = axes[2]
from collections import Counter
feat_counts = Counter(f for feats in feature_history for f in feats)
ax.bar(feat_counts.keys(), feat_counts.values(), color='teal', alpha=0.7)
ax.set_xlabel("Feature")
ax.set_ylabel("Times Selected")
ax.set_title("Feature Selection Frequency")
ax.tick_params(axis='x', rotation=30)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Monitoring Convergence with Progress Callbacks¶
Progress callbacks let you track acceptance rate, feature selection, and space reduction across iterations — useful for deciding when the workflow has converged.
In [17]:
Copied!
# Save current state as checkpoint
checkpoint_path = Path('./advanced_engine_checkpoint.pkl')
advanced_engine.save_checkpoint(checkpoint_path)
print(f"Checkpoint saved to {checkpoint_path}")
# Demonstrate loading from checkpoint
print(f"\nLoading engine from checkpoint...")
# Create new strategies for loading (required parameters)
sampling_strategy = hm.SamplingStrategyFactory.create('lhs')
feature_strategy = hm.AutoFeatureSelection(method='fano', max_features=3)
emulator_factory = hm.EmulatorFactory('gpr')
loaded_engine = hm.HistoryMatching.load_checkpoint(
checkpoint_path,
sampling_strategy=sampling_strategy,
feature_selection=feature_strategy,
emulator_factory=emulator_factory
)
loaded_engine.function = enhanced_sir_simulation
print(f"Engine loaded successfully!")
print(f" Current iteration: {loaded_engine.current_iteration}")
print(f" Total samples accepted: {loaded_engine.samples_accepted}")
print(f" Emulators trained: {loaded_engine.emulators_trained}")
# Can continue from where we left off
if loaded_engine.current_iteration < loaded_engine.max_iterations:
print(f"\nContinuing workflow from checkpoint...")
# Run one more iteration to demonstrate
next_result = loaded_engine.step()
print(f" Additional iteration completed: {len(next_result.samples)} samples")
loaded_engine.commit_step()
print(f" Updated progress: Iteration {loaded_engine.current_iteration}, "
f"Rate: {loaded_engine.acceptance_rate:.3f}")
# Clean up checkpoint file
checkpoint_path.unlink()
print(f"\nCheckpoint file cleaned up")
# Save current state as checkpoint
checkpoint_path = Path('./advanced_engine_checkpoint.pkl')
advanced_engine.save_checkpoint(checkpoint_path)
print(f"Checkpoint saved to {checkpoint_path}")
# Demonstrate loading from checkpoint
print(f"\nLoading engine from checkpoint...")
# Create new strategies for loading (required parameters)
sampling_strategy = hm.SamplingStrategyFactory.create('lhs')
feature_strategy = hm.AutoFeatureSelection(method='fano', max_features=3)
emulator_factory = hm.EmulatorFactory('gpr')
loaded_engine = hm.HistoryMatching.load_checkpoint(
checkpoint_path,
sampling_strategy=sampling_strategy,
feature_selection=feature_strategy,
emulator_factory=emulator_factory
)
loaded_engine.function = enhanced_sir_simulation
print(f"Engine loaded successfully!")
print(f" Current iteration: {loaded_engine.current_iteration}")
print(f" Total samples accepted: {loaded_engine.samples_accepted}")
print(f" Emulators trained: {loaded_engine.emulators_trained}")
# Can continue from where we left off
if loaded_engine.current_iteration < loaded_engine.max_iterations:
print(f"\nContinuing workflow from checkpoint...")
# Run one more iteration to demonstrate
next_result = loaded_engine.step()
print(f" Additional iteration completed: {len(next_result.samples)} samples")
loaded_engine.commit_step()
print(f" Updated progress: Iteration {loaded_engine.current_iteration}, "
f"Rate: {loaded_engine.acceptance_rate:.3f}")
# Clean up checkpoint file
checkpoint_path.unlink()
print(f"\nCheckpoint file cleaned up")
INFO:historymatching.engine:Checkpoint saved to advanced_engine_checkpoint.pkl
INFO:historymatching.engine:Engine loaded from checkpoint advanced_engine_checkpoint.pkl
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:WAVE 4 STARTING
INFO:historymatching.engine:============================================================
INFO:historymatching.engine:[Wave 4] Phase 1/5 SAMPLING: using 800 pre-computed samples from wave 3 [0.0s]
INFO:historymatching.engine:[Wave 4] Phase 2/5 SIMULATION: running 800 simulations...
Checkpoint saved to advanced_engine_checkpoint.pkl Loading engine from checkpoint... Engine loaded successfully! Current iteration: 3 Total samples accepted: 4000 Emulators trained: 7 Continuing workflow from checkpoint...
INFO:historymatching.engine:[Wave 4] Phase 2/5 SIMULATION: complete — 31 outputs [0.6s]
INFO:historymatching.feature_selection: Feature ranking (fano): cases_week_1=520.67, cases_week_2=159.13, total_cases=35.21, attack_rate=35.21, incidence_day_8=33.67
INFO:historymatching.feature_selection: SELECTED: cases_week_1 (score=520.67, mean_z=-0.04, std_z=4.29)
INFO:historymatching.feature_selection: SELECTED: cases_week_2 (score=159.13, mean_z=0.03, std_z=2.10)
INFO:historymatching.feature_selection: SELECTED: incidence_day_8 (score=33.67, mean_z=0.36, std_z=3.49)
INFO:historymatching.engine:[Wave 4] Phase 3/5 FEATURE SELECTION: ['cases_week_1', 'cases_week_2', 'incidence_day_8'] [0.0s]
INFO:historymatching.engine:[Wave 4] Phase 4/5 EMULATOR TRAINING: training 3 emulators...
INFO:historymatching.emulators.factory: Emulator 1/3 [cases_week_1]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: cases_week_1
INFO:historymatching.emulators.factory: Emulator 1/3 [cases_week_1]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.4664, lengthscale range=[2.873, 3.569] [2.4s]
INFO:historymatching.emulators.factory: Emulator 1/3 [cases_week_1]: trained [2.4s]
INFO:historymatching.emulators.factory: Emulator 2/3 [cases_week_2]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: cases_week_2
INFO:historymatching.emulators.factory: Emulator 2/3 [cases_week_2]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.0289, lengthscale range=[1.159, 1.191] [2.1s]
INFO:historymatching.emulators.factory: Emulator 2/3 [cases_week_2]: trained [2.1s]
INFO:historymatching.emulators.factory: Emulator 3/3 [incidence_day_8]: creating (800 samples, 2 params, type=gpr)...
DEBUG:historymatching.emulators.factory:Created gpr emulator for feature: incidence_day_8
INFO:historymatching.emulators.factory: Emulator 3/3 [incidence_day_8]: training...
INFO:historymatching.emulators.gpr: Training GPR: 600 points, 2 dims
INFO:historymatching.emulators.gpr: Building GPflow model and optimizing hyperparameters...
INFO:historymatching.emulators.gpr: GPR training complete — noise=0.3203, lengthscale range=[1.648, 2.400] [1.8s]
INFO:historymatching.emulators.factory: Emulator 3/3 [incidence_day_8]: trained [1.8s]
INFO:historymatching.engine:[Wave 4] Phase 4/5 EMULATOR TRAINING: complete [6.3s]
INFO:historymatching.engine:[Wave 4] Phase 5/5 NROY SAMPLING: finding 800 plausible candidates for next wave...
INFO:historymatching.nroy_sampling: LHS rejection: 0/800 (0%) — starting
INFO:historymatching.nroy_sampling: LHS rejection: 8/800 (1%) | 880 tested | rate=0.9091% [1s]
INFO:historymatching.nroy_sampling: LHS rejection: 1178/800 (147%) | 96,712 tested | rate=1.2180% [7s]
INFO:historymatching.nroy_sampling: LHS rejection: 800/800 [7.2s]
INFO:historymatching.nroy_sampling: NROY SUMMARY: 800/800 [OK] — sources: LHS=800 [7.2s]
DEBUG:historymatching.engine:NROY fraction: 0.008272
INFO:historymatching.engine:[Wave 4] Phase 5/5 NROY SAMPLING: found 800 candidates [7.2s]
INFO:historymatching.engine:[Wave 4] ALL PHASES COMPLETE [14.2s total]. Committing...
INFO:historymatching.engine:[Wave 4] COMMITTED — diagnostics and checkpoint saved to None
INFO:historymatching.engine:============================================================
Additional iteration completed: 800 samples Updated progress: Iteration 4, Rate: 0.022 Checkpoint file cleaned up
Summary¶
This notebook demonstrated advanced features of the history matching library:
Advanced Configuration¶
- Single-constructor API: Fine-grained control over all components
- Custom Strategies: Sampling, feature selection, and emulator configuration
- DataFrame Inputs: Structured parameter and observation definitions
- Performance Tuning: Threshold, oversampling, and optimization settings
Interactive Workflow Control¶
- Step-by-Step Execution: Manual control over each iteration
- Strategy Switching: Dynamic reconfiguration during workflow
- Rollback Capability: Revert unsatisfactory iterations
- Progress Monitoring: Real-time callbacks and diagnostics
Advanced Analysis¶
- Comprehensive Diagnostics: Parameter recovery assessment
- Performance Comparison: Different configuration strategies
- Evolution Tracking: Parameter space reduction over iterations
- Feature Selection Analysis: Understanding which features matter most
Enterprise Features¶
- Checkpoint/Resume: Save and restore workflow state
- Extensible Design: Custom emulators and strategies
- Scalable Architecture: Handle complex, high-dimensional problems
Key Takeaways¶
- Emulator Choice Matters: GPR generally outperforms linear models for complex relationships
- Feature Selection Strategy: Automatic selection often finds good features, manual gives control
- Sampling Strategy: LHS typically provides better space coverage than random sampling
- Threshold Tuning: Conservative thresholds (≤3.0) balance precision and coverage
- Interactive Workflow: Powerful for research and model development scenarios
The single-constructor API provides the flexibility needed for advanced applications while maintaining ease of use for standard workflows.