Address people by an ancestor geographical unit, and date symptom policies

CMakeLists.txt

@@ -147,6 +147,7 @@ set(JUNE_CORE_SOURCES
     src/utils/filtered_csv.cpp
     src/utils/event_logging/event_writer.cpp
     src/utils/event_logging/event_writer_lookups.cpp
+    src/utils/mpi_logging.cpp
     src/utils/event_logging/event_merger.cpp
     src/utils/event_logging/event_merger_lookups.cpp
     src/epidemiology/epidemiology.cpp
@@ -294,6 +295,7 @@ if(BUILD_TESTING)
   add_june_test(test_domain_communicator_detail "tests/test_domain_communicator_detail.cpp")
   add_june_test(test_run_dir "tests/test_run_dir.cpp")
   add_june_test(test_world_state "tests/test_world_state.cpp")
+  add_june_test(test_selection_criterion_geography "tests/test_selection_criterion_geography.cpp")
   add_june_test(test_domain_loader_internals "tests/test_domain_loader_internals.cpp")
   add_june_test(test_infection "tests/test_infection.cpp;src/loaders/disease_loader.cpp")
   add_june_test(test_transmission_lookups "tests/test_transmission_lookups.cpp")

analysis_tools/symptom_policy_window_check.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+"""Check that a windowed symptom policy does nothing before its start day.
+
+Compares two runs of the same config that differ only in their symptom
+policies: one with the windowed self-isolation policies, one with them removed.
+The engine records no policy-relocation event, so the observable signature of a
+symptom policy is the infections it prevents. Identical infection streams up to
+the window's start day, differing streams from it, is the demonstration.
+
+Producing the two runs (config_2021 starts 2020-02-01, so 12 Mar 2020 = day 40):
+
+    ./build/disease_sim --config configs/config_2021/simulation.yaml \
+        --world world_state.h5 --days 50 --runs-dir runs --run-id windowed
+
+    # control: same config with every isolation policy deleted, leaving only
+    # hospitalize_severe_cases. Keep the copy where its data/ references still
+    # resolve (the loader reads them relative to the config file).
+    ./build/disease_sim --config <copy>/simulation.yaml \
+        --world world_state.h5 --days 50 --runs-dir runs --run-id no_policy
+
+Usage:
+    python analysis_tools/symptom_policy_window_check.py \
+        runs/windowed runs/no_policy --window-start-day 40
+
+Result on 2026-08-19: identical infection streams on days 0-39, first
+difference on day 40 (15584 infections windowed vs 17101 without).
+"""
+
+import argparse
+import sys
+from collections import Counter
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+
+def infections_by_day(run_dir):
+    """Multiset of infection events, keyed by whole simulation day."""
+    events_path = Path(run_dir) / "simulation_events.h5"
+    with h5py.File(events_path, "r") as events_file:
+        infections = events_file["events/infections"][:]
+
+    by_day = {}
+    for day in np.unique(np.floor(infections["time"]).astype(int)):
+        on_day = infections[np.floor(infections["time"]).astype(int) == day]
+        by_day[int(day)] = Counter(
+            zip(
+                on_day["person_id"].tolist(),
+                on_day["infector_id"].tolist(),
+                on_day["venue_id"].tolist(),
+                on_day["time"].tolist(),
+            )
+        )
+    return by_day
+
+
+def first_differing_day(windowed, no_policy):
+    for day in sorted(set(windowed) | set(no_policy)):
+        if windowed.get(day, Counter()) != no_policy.get(day, Counter()):
+            return day
+    return None
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("windowed_run_dir")
+    parser.add_argument("no_policy_run_dir")
+    parser.add_argument("--window-start-day", type=int, required=True)
+    arguments = parser.parse_args()
+
+    windowed = infections_by_day(arguments.windowed_run_dir)
+    no_policy = infections_by_day(arguments.no_policy_run_dir)
+
+    last_day = max(max(windowed), max(no_policy))
+    divergence_day = first_differing_day(windowed, no_policy)
+
+    print(f"days compared:      0 - {last_day}")
+    print(f"window starts:      day {arguments.window_start_day}")
+    print(f"first differing day: {divergence_day}")
+
+    for day in range(arguments.window_start_day - 2,
+                     min(arguments.window_start_day + 3, last_day + 1)):
+        print(f"  day {day:>3}: windowed={sum(windowed.get(day, Counter()).values()):>6}"
+              f"  no_policy={sum(no_policy.get(day, Counter()).values()):>6}")
+
+    if divergence_day is None:
+        print("FAIL: the policies changed nothing at all")
+        return 1
+    if divergence_day < arguments.window_start_day:
+        print(f"FAIL: the policies acted on day {divergence_day}, "
+              f"before their window opened")
+        return 1
+    if divergence_day > arguments.window_start_day:
+        print(f"WARN: no effect until day {divergence_day}; the window opened on "
+              f"day {arguments.window_start_day}. Expected if nobody was "
+              f"symptomatic and out of the house on the opening day.")
+    print("PASS: no effect before the window, effect from it")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

configs/config_2021/infection_seeds.yaml

@@ -9,7 +9,7 @@ infection_seeds:
   # Single uniform seeding event
   - name: "basic_seeding"
     type: "uniform"
-    date: "2024-01-01 08:00"
+    date: "2020-02-01 08:00"
     parameters:
       cases_per_capita_multiplier: 15.0
 

configs/config_2021/policies.yaml

@@ -20,11 +20,23 @@ policies:
   symptom_policies:
 
     # -------------------------------------------------------------------------
-    # Policy: Mild/Severe symptoms → Stay at residence
+    # 12 Mar 2020 (day 40) - "contain" becomes "delay". Anyone with a new
+    # continuous cough or fever is told to self-isolate for 7 days. There is
+    # deliberately no policy before this date: symptomatic people behaved
+    # normally through February, which is what the timeline says happened.
+    #
+    # The window closes on 16 Mar, when household isolation supersedes it. The
+    # follow-up chain mild -> severe lives entirely inside this window: a
+    # closing window is a change of government instruction, not symptom
+    # progression, so nothing is inherited across the 16 Mar edge. Compliance
+    # is re-rolled against the superseding policy's own rate.
     # -------------------------------------------------------------------------
 
     - name: "isolate_when_mild"
 
+      start_date: "2020-03-12"
+      end_date: "2020-03-16"
+
       # Symptoms that trigger this policy
       symptoms: ["mild"]
 
@@ -42,35 +54,70 @@ policies:
       inherit_compliance: true
       inherit_refusal: false # Let them change their mind for severe symptoms
 
-      # Optional: Only apply to certain people (using selection criteria)
-      # applies_to:
-      #   - property: "age"
-      #     operator: "<"
-      #     value: 65
-
     - name: "isolate_when_severe"
 
-      # Symptoms that trigger this policy
+      start_date: "2020-03-12"
+      end_date: "2020-03-16"
+
       symptoms: ["severe"]
 
-      # Override these activities when symptoms trigger
-      # Options: List of activity names, or "*" for all activities
       override_activities: ["primary_activity", "leisure"]
 
-      # Replace with this activity
       replacement: "residence"
 
       compliance_rate: 0.8
 
-      # Optional: Only apply to certain people (using selection criteria)
-      # applies_to:
-      #   - property: "age"
-      #     operator: "<"
-      #     value: 65
+    # -------------------------------------------------------------------------
+    # 16 Mar 2020 (day 44) - household isolation. The whole household of a
+    # symptomatic person is told to stay at home for 14 days.
+    #
+    # A SymptomPolicy acts on the person carrying the symptom, so the
+    # household members are not themselves confined here - what the stronger
+    # instruction buys is a higher compliance rate for the symptomatic person.
+    # Confining well household members would need a policy trigger keyed on a
+    # housemate's symptoms, which the engine does not have.
+    #
+    # No end_date: the measure runs to the end of the simulated period. The
+    # mild -> severe follow-up chain is declared again inside this window, so
+    # symptom progression on or after 16 Mar inherits within the window.
+    # -------------------------------------------------------------------------
+
+    - name: "isolate_household_when_mild"
+
+      start_date: "2020-03-16"
+
+      symptoms: ["mild"]
+
+      override_activities: ["primary_activity", "leisure"]
+
+      replacement: "residence"
+
+      # Placeholder rates: household isolation was a firmer instruction than
+      # the 12 Mar advice, so compliance is higher. Fit against data before
+      # drawing conclusions from a run.
+      compliance_rate: 0.7
+
+      follow_up_policy: "isolate_household_when_severe"
+      inherit_compliance: true
+      inherit_refusal: false
+
+    - name: "isolate_household_when_severe"
+
+      start_date: "2020-03-16"
+
+      symptoms: ["severe"]
+
+      override_activities: ["primary_activity", "leisure"]
+
+      replacement: "residence"
+
+      compliance_rate: 0.9
 
     # -------------------------------------------------------------------------
-    # Policy: Hospitalized → Go to medical facility
+    # Policy: Hospitalized -> Go to medical facility
     # -------------------------------------------------------------------------
+    # Undated: hospitalisation is clinical care, not an intervention of the
+    # period, so it is in force for the whole run.
     - name: "hospitalize_severe_cases"
 
       symptoms: ["hospitalised", "intensive_care"]

configs/config_2021/schedules.yaml

@@ -1,7 +1,7 @@
 # Schedule Configuration - With Aligned Boundaries
 # All schedule types share the SAME time slot boundaries
 
-day_type_cycle: [weekday, weekday, weekday, weekday, weekday, weekend, weekend]  # Mon-Fri, Sat-Sun
+day_type_cycle: [weekend, weekend, weekday, weekday, weekday, weekday, weekday]  # Sat-Sun, Mon-Fri; start_date 2020-02-01 is a Saturday
 
 default_schedule_type: "has_primary_activity"
 

configs/config_2021/simulation.yaml

@@ -1,6 +1,10 @@
 time:
-    start_date : "2024-01-01"
-    end_date : "2024-07-01"
+    # The COVID-19 UK timeline this config reproduces. End date is exclusive,
+    # so the run covers 1 Feb 2020 through 31 Jul 2020. 12 Mar 2020 (the
+    # "contain" -> "delay" switch) is day 40; 16 Mar 2020 (household
+    # isolation) is day 44.
+    start_date : "2020-02-01"
+    end_date : "2020-08-01"
 
 config_paths:
 # Configuration paths
@@ -60,5 +64,5 @@ checkpoint:
     enabled: true
     output_dir: checkpoints/
     every_n_days: null          # checkpoint weekly; ignored if on_dates is set
-    on_dates: ["2024-02-05"]           # or a list, e.g. ["2024-01-10", "2024-01-20"]
+    on_dates: ["2020-03-12"]           # or a list, e.g. ["2020-03-12", "2020-03-16"]
     keep_last: 3             # 0 = keep all

include/core/config.h

@@ -45,7 +45,25 @@ struct SelectionCriterion {
   // block in the error message.
   void resolveOrThrow(const WorldState& world, const std::string& context);
 
+  // True for the property paths whose `value` is a geographical unit name (or
+  // a list of them) rather than a number: `geo_unit.<LEVEL>` and nothing else.
+  // Lexical, so config loaders can ask before a world exists. Kept beside the
+  // path dispatch in resolve() so the two cannot drift.
+  static bool comparesAgainstUnitNames(const std::string& property_path);
+
+  // Reference tables can name geographical units a given world does not have,
+  // e.g. a national table run on a regional world. With this set, a
+  // geo_unit.<LEVEL> name found at no level at all is recorded instead of being
+  // a resolve error, and matches nobody. A name that exists at a different
+  // level is still an error: that is a misspelt level, not a missing place.
+  bool allow_absent_geo_units = false;
+  const std::vector<std::string>& absentGeoUnitNames() const {
+    return absent_geo_unit_names;
+  }
+
  private:
+  mutable std::vector<std::string> absent_geo_unit_names;
+
   enum class PropertyType {
     UNKNOWN,
     AGE,
@@ -59,7 +77,30 @@ struct SelectionCriterion {
     PARTNER_IN_NETWORK,
     // is_alive: convenience for `!person.is_dead`. Path: "is_alive".
     IS_ALIVE,
+    // geo_unit.<LEVEL>: the person's ancestor geographical unit at a named
+    // level, compared by unit name. Path: "geo_unit.XLGU".
+    GEO_ANCESTOR,
   };
+
+  // Ancestor-geography membership, one entry per geographical unit, built
+  // once at resolve time so evaluate is a single array read: 0 = under some
+  // other unit, 1 = under a target unit, 2 = no ancestor at this level
+  // (absent — false whichever way the criterion is written).
+  mutable std::vector<uint8_t> geo_ancestor_mask;
+  // Empty when the mask is keyed by geo_unit_id directly (ids dense enough to
+  // make that cheap). Otherwise holds the world's unit ids in sorted order and
+  // the mask is keyed by position in it, found by binary search. Both forms are
+  // self-contained: evaluate answers from the criterion alone, never from a
+  // WorldState it may not have been handed.
+  mutable std::vector<GeoUnitId> geo_mask_unit_ids;
+  // Recorded rather than thrown: evaluate resolves lazily and must not throw
+  // from the hot path, so resolveOrThrow is what turns this into an error.
+  mutable std::string geo_resolve_error;
+
+  void buildGeoAncestorMask(const WorldState& world) const;
+  // Position of `id` in geo_ancestor_mask, or geo_ancestor_mask.size() when the
+  // mask has no entry for it.
+  size_t geoMaskSlot(GeoUnitId id) const;
   mutable PropertyType cached_type = PropertyType::UNKNOWN;
   mutable std::string cached_activity_name;  // (also reused for facet name)
   mutable std::string cached_sub_property;   // (also reused for facet field)
@@ -224,7 +265,8 @@ struct ScheduleType {
 
   void resolve(const WorldState& world) {
     for (auto& criterion : selection_criteria) {
-      criterion.resolve(world);
+      criterion.resolveOrThrow(world,
+                               "schedule type '" + name + "' selection");
     }
     // force_hybrid_mask is resolved in ScheduleConfig::resolveSlots (defined
     // in config.cpp where WorldState is complete).
@@ -352,7 +394,8 @@ struct VaccinationCampaignConfig {
 
   void resolve(const WorldState& world) {
     for (auto& crit : selection_criteria) {
-      crit.resolve(world);
+      crit.resolveOrThrow(world,
+                          "vaccination campaign '" + name + "' selection");
     }
   }
 };

include/core/variant.h

@@ -13,7 +13,8 @@ using PropertyValue =
                  int32_t,              // Integer
                  double,               // Double
                  std::string,          // String
-                 std::vector<int32_t>  // List of ints (e.g., social_contacts)
+                 std::vector<int32_t>,     // List of ints (e.g. social_contacts)
+                 std::vector<std::string>  // List of strings (e.g. unit names)
                  >;
 
 // Helper functions for PropertyValue

include/core/world_state.h

@@ -9,6 +9,7 @@
 #include <sstream>
 #include <string>
 #include <unordered_map>
+#include <unordered_set>
 #include <vector>
 
 #include "epidemiology/disease.h"
@@ -229,6 +230,12 @@ class WorldState {
   // descendants
   std::unordered_map<GeoUnitId, std::vector<uint32_t>> people_by_geo_unit;
 
+  // The units people are assigned to directly, i.e. the set of values of
+  // Person::geo_unit_id. A strict subset of people_by_geo_unit's keys, which
+  // also carry every ancestor of those units. Diagnostics that ask "would this
+  // exclude anybody?" want this one.
+  std::unordered_set<GeoUnitId> directly_inhabited_geo_units;
+
   // Build lookup indices (call after loading)
   void buildIndices();