ITG baseline: ~85% of its travel cost is the unreachable-client penalty — landbody-connectivity fix (79cb58b) + recalculation plan

[MaartenHilferink]

Problem: ITG's baseline travel cost is an extreme outlier — and it is mostly an artefact

ITG (Sicilia + Sardegna) has by far the worst baseline of all 42 study areas:

	ITG	next worst	typical
baseline mean travel time	22.30 min	9.46 min (SE3)	2–5 min
baseline travel cost / facility	€157,742	€80,549 (Denmark)	€10–50k
rect_rel (relative improvement potential, LINEAR)	0.384	0.320 (Portugal)	0.005–0.15

Taken at face value this says Sicily and Sardinia have catastrophic pharmacy access and
the largest optimisation headroom in the study. They do not. The number is dominated
by clients that the model cannot route to any pharmacy, which soft coverage prices at
BIG = c(120 min).

Decomposition (from logs/sweep_ITG_LINEAR.log)

total client population:  6,253,866
unreachable clients:      1,859 cells / 983,169 residents (15.7%) priced at BIG=119.9638
total travel cost(c):     1.39443726e8
mean travel time (min):   22.2972

Splitting that total:

component	value	share
BIG penalty on unreachable residents (983,169 × 119.96)	117,944,689	84.6 %
actual travel of the 5,270,697 served residents	21,499,037	15.4 %

⇒ the served population of ITG travels 4.08 min on average — entirely normal, right
next to ITI (4.18), Czechia (4.23) and Austria (4.82). Roughly 85 % of ITG's baseline
"travel cost" is not travel; it is the penalty for residents the network cannot reach.

This is self-consistent with the optimisation: at S1 (same facility count, optimally
relocated) the reported mean drops from 22.30 to 4.62 min. The optimiser is not
achieving a 5× improvement in accessibility — it is mostly reshuffling facilities so that
fewer people fall into the unreachable bucket.

ITG is the worst case, but not the only one

Same decomposition across all areas (LINEAR):

area	population	unreachable	%	baseline mean_t	BIG share of cost	served mean_t
ITG	6,253,866	983,169	15.7 %	22.30	84.6 %	4.08
Portugal	10,300,382	459,683	4.5 %	8.56	62.6 %	3.35
PL8	4,674,751	153,669	3.3 %	8.89	44.4 %	5.11
SE2	4,398,063	32,226	0.7 %	8.05	10.9 %	7.23
Norway	5,110,060	3,696	0.1 %	7.99	1.1 %	7.91
(all other areas)		≈0	0 %	1.9–5.6	≈0 %	= mean_t

Portugal's 459,683 unreachable residents match the combined population of the Azores and
Madeira almost exactly — separate landbodies, same mechanism. PL8's 153,669 needs a
separate look, since Poland has no comparable islands; it may be a genuinely disconnected
road component rather than an island.

Consequence: the ranked "improvement potential" lists (deck p64–p69) and the
rect_rel ordering are currently ranking network-coverage artefacts alongside real
policy headroom. ITG tops the LINEAR ranking for the wrong reason.

---

Proposed solution: the landbody-connectivity fix Chris added in 79cb58b

Commit 79cb58b ("Additional reporting elements") on branch ServiceAccess already
contains what looks like the right fix. In cfg/main/Templates.dms the network-selection
criterion changed:

- attribute<bool> Roads_isConnected (…) := …/Connectiveness/IsSterkVerbonden;
+ attribute<bool> Roads_isConnected (…) := …/Connectiveness_seldomain/
+                                          Strong_connectivity/IsConnected_seldomain;

with a new template Check_Connectiveness_T_seldomain, parameterised on
seldomain = /SourceData/RegionalUnits/Country_Split, which is

Country_Split := geos_split_union_polygon(Country/geometry, id(Country))

i.e. each country's multipolygon split into its separate landbodies.

Before: only the single largest strongly-connected component in the whole study-area
road extract survived; every island subnetwork was pruned. For ITG that is fatal — Sicilia
and Sardegna are two separate landbodies, so at most one of them could ever be kept, and
the smaller islands never.

After: the largest SCC is selected per landbody (IsConnected_seldomain via
Networks_per_seldomain / Main_in_seldomain), so each island keeps its own road network,
and the ferry links inside those components come back into play.

Worth stressing: no source data changed. The ferry and road links were already in the
TomTom NW2021_SP_streets_subset extract — they were being discarded together with the
components that got pruned. This is purely the selection rule, which is why the fix is
cheap and why no re-import is needed.

---

What this invalidates

The network is the first link in the chain, so the change propagates all the way through:

road network → OD matrices (existing + candidate) → baseline metrics
             → λ-sweep → S1/S2 + frontier → metrics/crossings → deck + exports

Concretely: the descriptives (deck p5–8), every sweep, S1/S2, the frontier and its
projections (p62/63), the ranked improvement lists (p64–69), the aggregate, the S1/S2
open-location Arrow files delivered in #45, and the CSVs delivered in #48.

Pre-fix and post-fix numbers will not be comparable and should not be mixed in one
chart or table.

---

Recalculation plan (steps of ≤ 24 h)

Full version committed as doc/replan_connectivity_fix.md. Sizing is anchored on
measured runtimes from the previous full recalculation (logs/recalc_q*.log):

stage	measured, all 42 areas
network + OD rebuild (both sides)	~1,400 s ≈ 0.4 h
λ-sweeps (LINEAR + LOGISTIC)	~81.5 h sequential

Rebuilding networks is nearly free; the sweeps are the entire bill. So the plan spends one
cheap step rebuilding all networks and diffing them, to decide which areas actually
moved, before spending any sweep time.

step	content	budget
0	Config decisions (see below) + snapshot current network/OD stats + end-to-end smoke test on ITG (island case) and a single-landbody control (expect no change)	≤ 4 h
1	Purge FinalSet_*, rebuild networks + ODs for all 42 areas, diff vs. snapshot → the affected-area table (Δlinks, Δnodes, ΔOD rows, Δunreachable). Areas with zero delta are excluded from re-sweeping. This is the decision point that sizes everything below.	≤ 12 h
2…k	Re-sweep affected areas via run_resweep_batch.ps1, 7 parallel workers, bin-packed on measured cost × 1.5 (reconnected islands add OD rows, so expect +20–50 %)	24 h × ~2–4
P	Poland country-level: 32.3 h measured, does not fit one step. Recommend dropping — it is already excluded from the aggregate (the 7 PL NUTS-1 are used, for disjointness) and only feeds its own deck page	optional
F1	Baseline + descriptives (p5–8). Expect the unreachable-resident figures to collapse — that is the visible payoff and should be called out explicitly	≤ 6 h
F2	build_deck_data.py → build_charts.py → frontier_metrics.py → build_deck.mjs → merge_deck.ps1, then doc/export_issue48.py; re-export the S1/S2 Arrow files	≤ 4 h
F3	Re-communicate: supersede #45 and #48, and note on deck p60 that the frontier now runs on landbody-complete networks	≤ 2 h

⚠️ In step 1, verify that run_recalc_batch.ps1's purge pattern (*set_O-1km*) really
catches the FinalSet_*.mmd directories on both sides. doc/todo.md §A1 records that
GeoDMS silently skips rewriting a stale .mmd, which would make the whole step a no-op.

---

Decisions needed before any batch run

1. Blocking — baseline and LP must run over the SAME client set. 79cb58b switched
   AllocateClientsToExistingPharmacies (the baseline) to Client_pharmacy_coverage
   (Analyses.dms:115), while AllocateClientsToNewPharmacies (the LP) still uses
   Client (:123). That re-opens the comparability defect fixed in doc/todo.md §B3.

   It is also more than a comparability issue. The same commit correctly generalised the
   template's attribute<Client> client_rel → attribute<Org> client_rel, so client_rel
   now follows whichever Org was passed. But the descriptives still hard-code Client:

   Analyses.dms:42  attribute<…/Allocation> best_od                     (Client) :=
                        min_index(…/Allocation/impedance, …/Allocation/client_rel);
   Analyses.dms:43  attribute<cells>        nearest_cell_by_road_client (Client) := …
   Analyses.dms:44  attribute<Client>       grid_client_rel             (grid)   := …
   

   These aggregate over the existing allocation, whose client_rel is now in
   Client_pharmacy_coverage, but they are declared over Client. Note this is not
   automatically a no-op for single-country areas: a client cell whose Country_rel is
   null (coastal/border cells that fall outside the country polygon) drops out of
   Client_pharmacy_coverage, so the two units can differ in count even for ITG — and
   where the counts happen to match, they are still distinct units.

   Proposed fix — one client unit, used everywhere. Introduce a single alias and route
   both allocations and the descriptives through it, so the two can never drift apart again:

   unit<uint32> ModelClient := Client_pharmacy_coverage;   // single source of truth
   

   then Analyses.dms:115 and :123 both take ModelClient, and :42–44 are retyped
   from Client to ModelClient. The s1/s2 reporting containers (:146, :153) should
   follow too if their travel times are to be compared with the baseline.

   The alternative — reverting :115 to Client — restores the pre-commit behaviour with
   zero risk, but drops the intent of excluding countries that have no pharmacy data.
   Recommendation: ModelClient = Client_pharmacy_coverage on both sides.
   @cjacobscrisioni — agreed, or was the coverage filter meant for multi-country areas only?

2. Circular dependency. The new services_allocated container reads the S1/S2 Arrow
   files (service_access_202608/<c>_s1_logistic_open.arrow) that this re-sweep will
   invalidate. It also lists 21 areas with a lowercase naming scheme, whereas the
   delivered set is 41 areas named <region>_<S1|S2>_<LINEAR|LOGISTIC>_open.arrow — as
   written it will silently read nothing. Naming and coverage need reconciling in step F2.

3. Minor. ModelParameters.dms StudyArea_default is now 'ITG'. Harmless for
   batches (the pipelines export STUDY_AREA), but it changes what an interactive GUI open
   loads — revert, or keep deliberately?

4. Reporting. Once stranding is largely gone, should the ranked improvement lists still
   report the BIG-penalty share per area, so a coverage artefact can never again masquerade
   as policy headroom?

cc @cjacobscrisioni

[MaartenHilferink]

Client vs ModelClient — intended difference and purpose

Following up on decision 1. Agreed: baseline and LP must run over the same client set.
Here is what each unit is for, so the distinction stays deliberate rather than accidental.

The two units answer different questions

	Client	ModelClient (= Client_pharmacy_coverage)
Question	Where do people live?	Where do people live and do we have the supply data to model them?
Definition	cells with ModelParameters/Client > 0 (i.e. population > 0 for pharmacy runs)	the same, and covered_country[Country_rel] — the cell's country contains ≥1 pharmacy in the study area
Depends on	population grid only	population grid + pharmacy source coverage
Role	the demand universe — a property of the geography	the modelled demand — a property of this run's data availability

ModelClient ⊆ Client by construction.

Why the distinction is worth keeping

The filter is not a modelling assumption about accessibility; it is an honest statement
about data coverage. In a multi-country study area (e.g. EU) some countries have no
pharmacy records at all. Those residents are not "badly served" — we simply do not know
how they are served. Including them would make every one of them unreachable and price
them at BIG, which would:

• inflate the baseline travel cost with pure data-absence,
• create the exact artefact this issue is about (see the ITG decomposition above — 84.6 %
  of ITG's baseline cost is a BIG penalty), and
• let a data gap masquerade as a policy finding in the ranked improvement lists.

So the separation is right. What was wrong is only that the two sides used different
units.

Why they must be the SAME unit on both sides

ModelClient is the correct choice for both, because the λ-sweep compares a baseline and
an optimum on the same demand:

• The baseline (AllocateClientsToExistingPharmacies) measures what the current
  network delivers.
• The LP (AllocateClientsToNewPharmacies) measures what an optimal network could
  deliver.

S1 ("same facility count") and S2 ("same travel cost") are both defined relative to the
baseline point. If the two sides carry different populations, S1/S2, the frontier, the
rectangle and its diagonal crossing are all comparing incomparable quantities — the defect
doc/todo.md §B3 fixed once already.

It is also a live correctness bug, not only a comparability concern

79cb58b correctly generalised the template's attribute<Client> client_rel →
attribute<Org> client_rel, so client_rel now follows whichever Org was passed. But
the descriptives still hard-code Client while reading the existing allocation, whose
client_rel is now in Client_pharmacy_coverage:

Analyses.dms:42  attribute<…/Allocation> best_od                     (Client) :=
                     min_index(…/Allocation/impedance, …/Allocation/client_rel);
Analyses.dms:43  attribute<cells>        nearest_cell_by_road_client (Client) := …
Analyses.dms:44  attribute<Client>       grid_client_rel             (grid)   := …

And this is not automatically a no-op for single-country areas: a client cell whose
Country_rel is null — coastal or border cells that fall just outside the country polygon —
drops out of Client_pharmacy_coverage. So the two units can differ in row count even for
ITG, and where the counts happen to coincide they are still distinct units.

Proposed shape

One alias as the single source of truth, so the two can never drift apart again:

// Client      : the demand universe — every inhabited cell in the study area.
// ModelClient : the demand we can legitimately MODEL — inhabited cells in countries
//               for which pharmacy supply data exists. Baseline and LP must both use
//               THIS unit, or S1/S2 and the frontier compare different populations.
unit<uint32> ModelClient := Client_pharmacy_coverage;

• Analyses.dms:115 (baseline) and :123 (LP) → both ModelClient
• Analyses.dms:42–44 (descriptives over the existing allocation) → retype to ModelClient
• :146 / :153 (the s1/s2 reporting containers) → follow, if their travel times are to be
  compared with the baseline

Client stays, and stays meaningful: it is the denominator for "share of residents we can
model", which is worth reporting per area precisely so a data gap can never again be read
as a policy result.

Open question for @cjacobscrisioni: should the candidate set (NewPharmacyLocations)
be restricted to covered countries as well? With λ > 0 the LP will not open a facility that
serves nobody, so it is harmless for the optimum — but leaving it unrestricted means the
reported candidate count M includes cells that can never be selected, which is confusing
in the logs.

---

Unrelated, but worth recording here

The pipelines defaulted GEODMS_EXE to C:\dev\GeoDMS_2026\bin\Release\x64\GeoDmsRun.exe
— the live build tree of the GeoDMS engine project. That couples two unrelated
workstreams in both directions: a batch run can pick up binaries that are mid-relink, and
the running GeoDmsRun holds a handle on Dm*.dll that makes the next engine link
silently skip. Repointed to the installed GeoDms 20.19.1.m in c368a01
(GEODMS_EXE still overrides). Verified: cfg/main.dms loads and /Analyses/Client
computes under 20.19.1.m.

[MaartenHilferink]

Step 1 done: 42/42 rebuilt — the fix works, and only 3 areas actually move. But it breaks Norway, SE2 and SE3.

Purged the cached FinalSet_*/Linkset_* artefacts on both sides and rebuilt network +
OD for every area under the landbody-connectivity fix, on the installed GeoDms 20.19.1.m.
Runner and results are committed (run_rebuild_all.ps1, scratch/rebuild_all_all.csv,
scratch/affected_areas.csv).

The fix delivers exactly what was predicted — ITG measured

ITG	before	after	change
unreachable residents	983,169 (15.7 %)	591 (0.009 %)	−99.94 %
baseline mean travel time	22.30 min	4.878 min	−78 %
baseline travel cost	1.394e8	3.063e7	−78 %
Existing OD rows	57,375	90,299	+57 %
candidate OD rows	505,944	809,249	+60 %

The decomposition in the issue body predicted the served population travels ≈4.08 min.
The new overall mean is 4.878 min — slightly higher, which is exactly right: the
newly-connected island residents travel further than the mainland average. ITG is no
longer an outlier.

Only 3 areas move at all

Candidate-OD growth after the rebuild, all 42 areas:

area	growth	pre-fix unreachable	pre-fix %	pre-fix mean_t
PL8	+27.9 %	153,669	3.3 %	8.89
Portugal	+10.7 %	459,683	4.5 %	8.56
ITG	+60 % (OD rows)	983,169	15.7 %	22.30
ITF	+0.1 %	0	0.0 %	6.67
(35 further areas)	+0.0 % … −0.1 %	0	0.0 %	1.03–6.52
FRH	−0.5 %	0	0.0 %	4.34
SE1	−2.9 %	6	0.0 %	5.41

ITG's byte-size row reads +0.0 % only because its snapshot was taken after it was
rebuilt first as the smoke test; its real growth comes from the sweep-log OD row counts.

This lines up precisely with the BIG-share table in the issue body: the areas that move
are exactly the ones whose baseline cost was dominated by the unreachable-client penalty
(ITG 84.6 %, Portugal 62.6 %, PL8 44.4 %). Every area with zero unreachable residents is
untouched. SE1's −2.9 % is the new candidate coverage filter, not connectivity.

So steps 2…k shrink from 2–4 × 24 h to a single short step: 3 areas, ~2.5 h of
sequential sweep at the previously measured cost.

Also worth noting for the mechanism: the Netherlands, with its Wadden islands and
ferries, shows 0.0 % growth. A ferry link joins the island to the mainland component,
so those islands were never pruned. The pruning only bit where a populated landbody was
connected to the main component by no modelled link at all — which is why ITG (whose
study area is only Sicilia + Sardegna, two separate landbodies) was the extreme case.

---

Defect: Norway, SE2 and SE3 now fail

All three fail at alloc, identically:

[[/NetworkSetup/NewPharmacy_Analysis/NetwerkSpec/CreateMoreEfficientNetwork/
   FinalSet_stored/FinalLinkSet/F2]] Link_Node2_rel: out of range or undefined

network1 and network2 succeed, so the .mmd is written; it breaks when the link set
is read back. These three swept fine before, so this is new.

Attribution established by experiment, not assumed. I reverted only the
connectivity criterion in Templates.dms — back to Connectiveness/IsSterkVerbonden and
Connectiveness/NodesSterkVerbonden — leaving the MMD read-holder split and every other
change in place, purged Norway and re-ran:

Norway, new side	result
per-landbody criterion (Strong_connectivity/IsConnected_seldomain)	fails
old single-largest-SCC criterion	exit 0, 51 s

So the cause is the connectivity change itself, not the read-holder split. The criterion
has been restored to the new one; the three areas stay failing pending your call.

Hypothesis, not verified — @cjacobscrisioni this is your algorithm, so I have not
touched it. In CreateInitialWorkingNetwork:

Roads_isConnected := …/Strong_connectivity/IsConnected_seldomain    // link: both endpoints kept
Nodes_isConnected := …/Strong_connectivity/NodesConnected_seldomain // node: kept
unit<uint32> Roads := select(Roads_isConnected) {
    unit<uint32> UqGeomPointSet := select(Nodes_isConnected) { attribute<uint64> values := JnctIds; }
    attribute<UqGeomPointSet> F1 := rlookup(F_JnctId, UqGeomPointSet/values);
    attribute<UqGeomPointSet> F2 := rlookup(T_JnctId, UqGeomPointSet/values);
}

The per-landbody rule retains nodes from several disconnected components at once, where
the old rule retained one. F1/F2 resolve through rlookup on JnctIds, which needs
those ids to be unique across the retained set — if two retained landbodies carry a
duplicate junction id, rlookup yields undefined and F2 goes out of range. Norway, SE2
and SE3 being the most island-rich areas fits that, but I have not confirmed it.

One more thing worth a look while you are in there: Strong_connectivity/Networks mixes
pcount(part_rel) with PartNr == Main and rlookup(PartNr, …). If
strongly_connected_components exposes part_rel rather than PartNr, those are not the
same attribute.

Where this leaves the plan

• Re-sweeping PL8, Portugal and ITG can start now — they rebuilt clean and are the
  only areas that moved.
• Norway, SE2 and SE3 are blocked on the defect above. They showed near-zero unreachable
  residents before (3,696 / 32,226 / 475), so they are not expected to move much; the
  blocker is correctness, not results.
• The other 36 areas need no re-sweep at all.

[MaartenHilferink]

Correction: Portugal is NOT fixed by the connectivity change — and the cause is not the road network

Re-swept the three areas that moved. Post-fix baselines (LINEAR):

area	unreachable before → after	mean_t before → after
ITG	983,169 (15.7 %) → 591	22.30 → 4.878
PL8	153,669 (3.3 %) → 0	8.89 → 4.534
Portugal	459,683 (4.5 %) → 467,704	8.56 → 8.644

ITG and PL8 are resolved, PL8 completely. Portugal is not, and my earlier framing —
that its stranding was "separate landbodies, same mechanism" — implied it would be. That
was wrong, so here is the actual cause.

The Azores and Madeira have no pharmacies in the source data

Probing the rebuilt Portugal artefacts directly (coordinates in EPSG:3035, km):

all clients          x:  944 .. 2976    y: 1506 .. 2788
unreachable clients  x:  944 .. 1884    y: 1506 .. 2788     (1,360 cells, 467,704 residents)
facilities (_j)      x: 2640 .. 2966    y: 1736 .. 2296     (1,893 pharmacies)

Every one of the 1,893 pharmacies sits at x ≥ 2640 km — mainland Portugal. Every
unreachable client sits at x ≤ 1884 km — the Azores and Madeira. The two sets do not
overlap at all.

So this was never a connectivity problem. No road-network fix can help: even a
perfectly retained island network has no pharmacy to route to. The islands are in the
population grid and in the study area, but absent from the pharmacy source data.

For contrast, ITG's residual 591 unreachable residents (12 cells) fall at x 4276–4530,
inside the range its 906 pharmacies span (4168–4818) — genuinely remote cells, not a data
gap. That is the expected residual.

Consequence: Portugal's baseline is inflated by pure data absence

467,704 residents × BIG(120 min) ≈ 5.6e7 of Portugal's 8.9e7 baseline travel cost — about
63 % — is the penalty for demand we have no supply data for. Exactly the failure mode
this issue was opened about, one level down: covered_country filters at country
granularity, and Portugal-the-country does have pharmacies, so its islands pass the filter
and are then priced as if catastrophically underserved.

Proposed fix: move the coverage filter from country to landbody

ModelClient currently uses

covered_country (Country/subset) := pcount(Pharmacies/within_StudyArea/Country_rel) > 0

The natural refinement is to evaluate coverage per Country_Split — the same
per-landbody unit Chris's connectivity fix already introduced — instead of per country:
a landbody with no pharmacy in the data is excluded from ModelClient and from
NewPharmacyLocations, exactly as an uncovered country is today.

That keeps the existing principle ("model demand only where we have supply data") and
just applies it at the granularity the data gap actually has. It would drop the Azores and
Madeira from Portugal's baseline and remove the ~63 % inflation.

This is a modelling-semantics decision, so I have not applied it. Two things to weigh:

1. It changes what "Portugal" means in the deck — 467k residents leave the denominator.
   That must be stated on the page, not silently applied.
2. The alternative reading is that the islands are genuinely underserved and should stay
   in. But we cannot tell that from this data: absence of pharmacy records is not evidence
   of absence of pharmacies. Pricing an unknown at BIG asserts the strong version.

@cjacobscrisioni — is the pharmacy source known to exclude the Azores and Madeira, or is
this an unexpected gap worth chasing upstream first?

Status of the plan

• ITG, PL8 — re-swept, resolved.
• Portugal — pipeline fine, result blocked on the decision above.
• Norway, SE2, SE3 — still blocked on the FinalLinkSet/F2 defect reported above.
• The other 36 areas — unchanged, no re-sweep needed.

F1 (descriptives) and F2 (deck + exports) should wait for Norway/SE2/SE3: the aggregate
sums over 41 disjoint areas, so it cannot be rebuilt while three of them are failing.

One measurement caveat

Client populations rose in all three areas (ITG +25,895; PL8 +382,273; Portugal +8,021).
That is expected rather than alarming — total client population counts clients present
in the OD, and a larger retained network admits more cells as OD clients. It does mean the
before/after unreachable counts are not exactly like-for-like; the direction is
unambiguous, the precise deltas are not. It also explains the ITG population rise I had
flagged earlier as needing verification: it is the network, not the Italian data update.

[MaartenHilferink]

Spun the Norway/SE2/SE3 failure off into its own issue, with the analysis.

Two corrections to what I proposed earlier in this thread, both ruled out:

1. Duplicate JnctIds — no. NodeSet_src := unique(PointSet/JnctIds) (Templates.dms:79) makes junction ids unique by construction, and Roads/UqGeomPointSet is a select() subset of it.
2. PartNr vs part_rel — no. The engine's ConnectedParts.cpp creates both (s_PartRel and the deprecated s_PartNr, GeoDMS #1177). Same data, two names — a style issue, not this bug.

The actual mechanism is node-id arithmetic: LinkSet_Calc packs three disjoint id ranges into one node space by hand (#Roads/UqGeomPointSet, + #UniqueLocations, + id(.)), and the per-landbody criterion grows the retained node set, shifting every offset. Details and the proposed per-connected-group offset fix are in the new issue.