Skip to content

Code walkthrough lp_run

Jip Claassens edited this page Sep 15, 2026 · 1 revision

Code walkthrough lp_run

In one sentence: lp_run.jl is the optimisation engine of the Julia side — it builds the soft-coverage facility-location LP once per study area, re-solves it at every price λ from the previous basis (the solver's record of the corner it finished at, so the next solve starts there instead of from scratch; Glossary), turns the fractional answer into a real set of open locations, and still carries the school-era LP that nothing on the pharmacy path uses any more.

Where this sits: GeoDMS exports the OD matrix (GeoDMS OD matrix and choice set) → settings.jl turns it into a CountryData (Code walkthrough settings) → lp_run.jl solves and rounds at one λ → lambda_sweep_simplex.jl walks the λ-grid and pins S1/S2 (Code walkthrough lambda_sweep_simplex) → logs and Arrow files (Log lines and output files) → deck (Deck and charts pipeline). The why of the model: The facility location model explained; of the rounding: From LP fractions to real pharmacies; of the warm start: The lambda sweep and the Pareto frontier.

Line numbers are for commit e25416b (10 September 2026, the branch head as of 15 September 2026) and drift; the permalinks in In the code are pinned.

One file, three eras

The file has 560 lines in three layers: the rounding helpers (June 2026) first, the legacy school-era LP run_scenario (March–May) in the middle, the live LP functions build_lp_warmstart, resolve_for_basis! and solve_at_w! (May, reworked July and September) at the bottom. Everything except run_scenario is on the pharmacy results path; the table in In the code dates each function. The file starts with include("settings.jl"), so c(t), big_cost(), FACILITY_MIN_COSTS, ROUNDING, the MS_* knobs and CountryData all come from there.

flowchart TD
  D["lambda_sweep_simplex.jl · analyze_country<br/>(the live driver)"] -->|once per study area| B[build_lp_warmstart]
  D -->|every grid point and every bisection step| S["solve_at_w!<br/>(rounding=true or false)"]
  D -->|once after a bisection detour| R[resolve_for_basis!]
  S --> G[greedy_round]
  S --> M[multistart_round]
  S --> T["travel_of (three times)"]
  S --> A[assign_nearest]
  M --> RS[randomized_seed]
  M --> CN[client_nearest]
  M --> SW["swap_round!"]
  M --> T
  T --> CN
  L["lp.jl (school-era runner)"] --> RSC["run_scenario (legacy)"]
  LS["lambda_sweep.jl (superseded sweep)"] --> RSC
  X["archive/test_greedy.jl (June one-off)"] -.-> B
  X -.-> S
Loading

The live driver calls build_lp_warmstart once and solve_at_w! many times on the same WarmStartState: consecutive λ share a basis, so every solve after the first is a warm start (Glossary). run_scenario is an island with its own model.

The structs

CountryData (settings.jl around line 208) is the input every function here takes: the OD matrix in column form — per OD row k (a client–candidate pair inside the choice set) the client clients_col[k], the candidate facilities_col[k], the travel time t_ij_col[k] in minutes (GeoDMS exports seconds; Log lines and output files) and the client's population wpop[k] — plus the look-ups locations[i] (the OD rows of client i, its choice set; Glossary), facility_rows[j] and client_pop[i]. N is the number of OD rows, M the number of candidates. Field by field: Code walkthrough settings.

WarmStartState{TX, TY} (around line 3) is what build_lp_warmstart returns and what solve_at_w! and resolve_for_basis! take:

struct WarmStartState{TX, TY}
    model::Model      # the JuMP model wrapping one live HiGHS instance
    y::TX             # openness variables, a DenseAxisArray keyed by facility id
    x::TY             # assignment variables, a Vector keyed by OD row 1..N
    data::CountryData
end

It exists because HiGHS keeps the optimal basis of the last solve inside the optimizer object: keep the Model alive, change only objective coefficients, and the next optimize! starts from that basis. The type parameters TX (for y) and TY (for x) are a leftover of the 24 July 2026 rename to the p-median convention (commit 4cee238; see Pitfalls): y is a location, x is a client share (The facility location model explained).

build_lp_warmstart(data) — build the LP once

Signature: build_lp_warmstart(data::CountryData) -> WarmStartState. Around lines 397–430.

HiGHS options set here: presolve = "on" (helps the cold first solve; warm solves skip it anyway; why it is left on is not recorded), output_flag = true (so the sweep log shows simplex iterations live; commit 011cafb: "IPM was silent during stalls") and solver = "simplex" (overwritten on every solve_at_w! from the SOLVER env var).

The variables

@variable(model, 0 <= x[1:N] <= 1)             # one per OD row: share of client i served by candidate j
@variable(model, 0 <= y[j in facilities] <= 1) # one per candidate: how "open" it is

x[k] is the fraction of OD row k's client sent to that row's candidate; y[j] is the openness of candidate j, relaxed from yes/no to a number in [0, 1]. That relaxation (Glossary) makes the problem an LP instead of a MIP, which is out of reach at this size; why that is acceptable: The facility location model explained, depth in Background theory.

The objective, exactly as coded

BIG = big_cost()
@objective(model, Min, sum(x[k] * (c(t_ij_col[k]) - BIG) * wpop[k] for k in 1:N))

Two things are deliberately absent: there is no λ and no y term, because the price of an open cell is the only thing that changes between sweep points and solve_at_w! adds it. And (c(t) − BIG) is negative for every row: under soft coverage (Glossary) a client's unserved share is charged BIG per resident (big_cost() = 120 under LINEAR, 1.0 under LOGISTIC, settings.jl around line 97; Travel cost functions), and multiplying that out leaves a per-client constant BIG·p_i that the minimiser ignores ("constant dropped from argmin", comment at line 413) — the derivation is on The facility location model explained. Folding the stranding price into the x-coefficient keeps all x-coefficients fixed across λ, so a sweep step touches only the M y-coefficients.

Tiny check. One client of 10 residents, two candidates at 20 and 50 minutes, LINEAR, BIG = 120. The coefficients are (20 − 120)·10 = −1 000 and (50 − 120)·10 = −700. Fully served by the near candidate the objective contributes −1 000; the dropped constant is 120·10 = 1 200; true travel = −1 000 + 1 200 = 200 person-minutes ✔. Left unserved the contribution is 0 and the true cost is the constant, 1 200 = 10 residents at BIG ✔.

Consequently the HiGHS "Objective value" line is Σ x(c−BIG)pop + λ·Σy, a large negative number under LINEAR — not a cost you can quote; quote travel_relax, sum_y and fac_€ from the sweep row instead (Log lines and output files).

The constraints

for (_, rows) in locations            # one per client:  Σ_{k in choice set} x[k] <= 1
    @constraint(model, sum(x[k] for k in rows) <= 1)
end
for k in 1:N                          # one per OD row:  x[k] <= y[candidate of k]
    @constraint(model, x[k] <= y[facilities_col[k]])
end

The first says a client is served at most once; the slack is its unserved share. It was == 1 (hard coverage) until commit 4ec47a9 (16 July 2026); why the change was needed (the full-coverage floor above today's count) is on The facility location model explained.

The second is the strong (disaggregated) formulation (Glossary): one x ≤ y per OD row, N of them, rather than one summed constraint per candidate. It forces y_j to be at least the largest share it serves, which pushes y to 0 or 1 in practice — the relaxation is nearly integral (old NL log: 0 to 22 fractional candidates out of 12 000–22 000 open; The facility location model explained). The price is size: Netherlands (old log) N = 1 001 117 OD rows → 1 024 025 rows and 1 030 910 columns, built in 4.2 s.

The function returns WarmStartState(model, y, x, data), unsolved. The explicit bound x ≤ 1 is redundant given the coverage constraint (rationale not recorded; probably inherited from run_scenario).

solve_at_w!(state, w; rounding=true) — one price, one point

Signature: solve_at_w!(state::WarmStartState, w::Real; rounding::Bool=true) -> NamedTuple. Around lines 460–560. The !: it mutates the model's objective and leaves the solver's basis behind.

What it rewrites

  1. λ = w * FACILITY_MIN_COSTS (line 464): w is the dimensionless knob, λ the price of one open cell in placeholder euro (€100 000 · w by default; The lambda sweep and the Pareto frontier).
  2. set_objective_coefficient(model, y[j], λ) for every candidate (lines 467–469). This is the only change to the model between sweep points. No constraint changed, so the previous answer is still feasible and the dual simplex (the HiGHS default) only re-checks which cells are still worth their new price — usually far cheaper than the cold solve (NL, old hard-coverage log: 45 876 iterations cold, then 1 369 and 7 716), but not at a large absolute step in λ: the high-λ "cliff" and the near-/far-basis mechanism are explained on The lambda sweep and the Pareto frontier.
  3. Solver options re-applied every call (lines 484–487): solver = env SOLVER (default "simplex"), run_crossover = "on" only when SOLVER=ipm, time_limit = LP_TIME_LIMIT (env, default 3 600 s). Crossover turns the interior solution of an interior-point method (IPM; Glossary) — thousands of tiny non-zero y, the "phantom fractions" of commit 51157a0 — into a corner solution with most y exactly 0 or 1.
  4. optimize!(model).
  5. If termination_status ∉ (OPTIMAL, LOCALLY_SOLVED, ALMOST_OPTIMAL) the function throws; the driver catches it, prints LP at w=… FAILED and skips the point. A TIME_LIMIT is the usual cause. Why these three statuses is not recorded.

The comment block at lines 471–483 records the solver history — the hard-coverage cliff (NL w = 5: 16.7 h), the July 2026 IPM detour and its retreat once it failed on soft-coverage LOGISTIC while the high-λ soft-coverage points became cheap anyway (few open cells ⇒ small basis). SOLVER=ipm remains as an escape hatch. Mechanism: The lambda sweep and the Pareto frontier; dated sequence: Decision log.

What it reads off

  • fractional = candidates with 1e-6 < y_relaxed < 1 − 1e-6; n_fractional_y = their count (log column frac_y).
  • sum_y = Σ_j y_relaxed[j] — the LP's fractional count of open cells. S1 is bracketed on it because it is monotone in λ while the rounded count is not (commit 9ee7ac6).
  • travel_relax (lines 506–513), the coverage-honest lower bound on travel:
served_relax   = Σ_k x_relaxed[k] · c(t_ij_col[k]) · wpop[k]
stranded_relax = Σ_i BIG · pop_i · (1 − Σ_{k ∈ locations[i]} x_relaxed[k])
travel_relax   = served_relax + stranded_relax

In words: the LP's travel with every unserved share charged at BIG — the HiGHS objective with the dropped constant added back and the λ·Σy part taken out. It certifies, on total cost, travel_relax + λ·sum_y ≤ integer optimum ≤ cost_c + λ·n_open; what that chain means, its worked example and what the reported "multi vs relax" column does not certify: From LP fractions to real pharmacies (and Pitfalls below).

The rounding=false path

Added in e25416b for issue #52. In refine-only mode (SWEEP_STOP_AFTER_REFINE=1) the driver walks the grid LP-only until a bracket can close, because the three roundings would be thrown away (Code walkthrough lambda_sweep_simplex owns the light walk; the concept is on Scenarios S1 S2 S3). The function returns right after travel_relax with placeholder rounding fields (see the table); the driver prints such light rows as "(LP only)" and writes no Arrow files for them.

The rounding=true path

All three roundings are always computed, whatever ROUNDING says (lines 531–535):

p_open      = clamp(round(Int, sum_y), 1, length(facilities))   # the count: round half to even
sorted_by_y = sort(collect(facilities), by=j -> y_relaxed[j], rev=true)
topp_set    = Set(sorted_by_y[1:p_open])                        # top-p
greedy_set  = greedy_round(y_relaxed, data, p_open)             # lazy greedy
multi_set   = multistart_round(y_relaxed, data, p_open, topp_set, greedy_set)

The count is fixed at the LP's own total so that the rounded point sits at the same place on the count axis as the LP point (issue #47; comment at lines 17–18; history since 51157a0 on From LP fractions to real pharmacies); the rounding only decides which p_open cells. Round half to even is Julia's default, not a choice.

Each set is scored with travel_of(set, data, BIG) → (travel with stranded residents at BIG, number of stranded client cells). ROUNDING (env, default "multistart", settings.jl around line 55) picks which set is the result: "greedy", "multistart", and anything else — including a typo — silently gives top-p (line 546 is the else branch). Finally assign_nearest(open_set, data) gives the per-client assignment for the Arrow output, and mean_t (line 551) divides assigned travel time by the total client population, so stranded residents count as zero minutes here but 120 min in baseline_metrics (Log lines and output files owns this asymmetry). Compare cost_c, not mean_t, with the baseline ★.

The result, field by field

Field Meaning Units Value when rounding=false
w grid value — as given
λ w · FACILITY_MIN_COSTS € (placeholder) as computed
n_open size of the selected open set cells clamp(round(Int, sum_y), 1, M)
cost_c travel_of of the selected set (stranded at BIG) person-minutes (LINEAR) or residents × c(t) ∈ [0, 1] (LOGISTIC) Inf
mean_t Σ assigned t·pop ÷ total pop minutes NaN
n_frac candidates with 1e-6 < y < 1 − 1e-6 candidates real value
sum_y Σ y_relaxed cells (fractional) real value
travel_relax LP travel + BIG × unserved share as cost_c real value
travel_c_topp, travel_c_greedy, travel_c_multi travel_of of each set as cost_c NaN
uncov_topp, uncov_greedy, uncov_multi stranded client cells per set (cells, not residents; zero-population cells included) cells -1
rounding value of ROUNDING — "none"
open_set selected candidate ids, Set{Int} — empty set
assigned_k client id → OD row of its nearest open candidate, Dict{Int,Int}; stranded clients absent — an empty Tuple{Int,Int}[] (a different type; nothing reads it for light rows)

The driver turns one result into one sweep-log row (columns: Log lines and output files; beware that fac_€ prices the fractional count sum_y) and, for rounded rows, into assignment.arrow and traveltime.arrow under <area>/lambda_sweep/<FUNC>/w=<w>/. The LP fractions themselves are never written to disk.

resolve_for_basis!(state, w) — put the basis back on the grid

Signature: resolve_for_basis!(state::WarmStartState, w::Real) -> termination status. Around lines 445–453.

It sets the y-coefficients to w · FACILITY_MIN_COSTS, calls optimize! and returns termination_status(model). Nothing else: no values read, no bound, no rounding; solver, run_crossover and time_limit stay whatever the last solve_at_w! set. The status is returned but not acted on.

Why it exists (issue #52, September 2026): S1 and S2 are pinned by bisection inside the grid loop (Glossary), which leaves the HiGHS basis at some w inside the bracket; carrying straight on to the next grid point would start from a far basis — exactly what produced the time-outs the fix was meant to cure. So the driver first re-solves the bracket's upper endpoint (resolve_for_basis!(state, r.w), printing "basis restored at w=…") and only then moves on, one near step away again. The driver side: Code walkthrough lambda_sweep_simplex; the FRC example (w = 0.1 in 75 s as a grid step, a 1 h time-out from a far basis): Scenarios S1 S2 S3.

The rounding functions, briefly

What each does and why, with a worked example, is on From LP fractions to real pharmacies. Here is only the interface; all take data::CountryData, measure travel with the same c(t) and use tol = 1e-6.

Function (around line) Signature Returns
greedy_round (23) (y_relaxed, data, p_open; tol=1e-6) Set{Int} of exactly p_open candidate ids
assign_nearest (91) (open_set, data) (Dict{Int,Int} client → OD row, travel of assigned clients only)
client_nearest (111) (open_set, data, BIG) (d1, d2, phi1): per client the cheapest and second-cheapest c(t) to an open candidate and the id of the nearest one (-1 if none; d1 = d2 = BIG then)
travel_of (136) (open_set, data, BIG) (Σ_i d1[i]·client_pop[i], number of clients with phi1 = −1)
randomized_seed (148) (must_open, frac_js, yvals, k, rng) Set{Int} of k fractional candidates
swap_round! (161) (chosen, frac_set, data, d1, d2, phi1) Bool — true if it performed one swap (mutates chosen)
multistart_round (211) (y_relaxed, data, p_open, topp_set, greedy_set; restarts=MS_RESTARTS, rounds=MS_ROUNDS, seed=MS_SEED, tol=1e-6) Set{Int} of p_open candidate ids

Code facts a caller must know: travel_of is the single metric every rounding is judged on, stranded residents charged BIG "so stranding the rural many-small-fraction clusters is never free" (comment at line 131) — except inside greedy_round, which prices an unreached client at Inf (Pitfalls). swap_round! applies one best-improving close/open swap per call, only if the exact change in travel_of is > 1e-6 (line 191); must-opens are never closed. multistart_round draws its random seeds from a MersenneTwister(MS_SEED) re-created at every call and keeps the best of all seeds, so travel_c_multi ≤ min(travel_c_topp, travel_c_greedy) always.

run_scenario — the legacy hard-coverage LP

Signature: run_scenario(data::CountryData, min_clients, w, apply_threshold, nearest) -> 16-tuple (open_set, fload, assigned_k, …, n_fractional_y, sum_y, travel_relax). Around lines 248–388. Only the school-era runners call it — lp.jl, and the superseded lambda_sweep.jl with min_clients = Inf; greedy.jl and greedy-merge.jl define their own functions of the same name. Nothing on the pharmacy results path calls it. It differs from build_lp_warmstart + solve_at_w! in four ways:

  1. A fresh model per call, IPM with crossover and output_flag = false (lines 253–257) — no warm start, no time limit.
  2. Hard coverage: sum(x[k] for k in rows) == 1 for every client (line 282); nobody can be stranded, so there is no BIG anywhere in the function — and a client whose candidates are all closed is silently absent from assigned_k, the blind spot that made travel_of necessary (commit e1c426b, 7 June 2026).
  3. A different cost of opening (lines 262–279): with a finite min_clients it is the school-era minimum-size term (deficit[j] ≥ min_clients·y[j] − load[j], priced at λ); with min_clients = Inf it is λ per open cell, i.e. the live model minus soft coverage.
  4. Top-p rounding only (lines 309–311), then nearest assignment (nearest=true, rows ≤ 60 min preferred when apply_threshold) or a re-solve of the fixed-y model (nearest=false; lines 322–360).

What the school-era term meant and why it was dropped: Alternatives not pursued and Service allocation procedure. cap_scenario.jl includes lp_run.jl only for settings.jl's helpers and builds its own LP (Alternatives not pursued).

In the code

Open lp_run.jl at build_lp_warmstart first, then solve_at_w!; read run_scenario last and only for contrast. Permalinks pinned to e25416b:

What Era Permalink
WarmStartState May 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L3
greedy_round June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L23
assign_nearest June 2026 (factored out of solve_at_w!) https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L91
client_nearest / travel_of June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L111
randomized_seed June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L148
swap_round! June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L161
multistart_round (the published rounding) June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L211
run_scenario (legacy) school era, March–May 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L248
build_lp_warmstart — objective and constraints May 2026, soft coverage July 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L397
LP_TIME_LIMIT and the solver-history comment July 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L432
resolve_for_basis! September 2026 (#52) https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L445
solve_at_w! May 2026, reworked July and September 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L460
travel_relax (coverage-honest bound) May 2026 (51157a0), soft-coverage form July 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L502
the rounding=false return September 2026 (#52) https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L515
the ROUNDING switch June 2026 https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lp_run.jl#L543
c(t) and big_cost() (settings.jl) — https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/settings.jl#L73
CountryData (settings.jl) — https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/settings.jl#L208
the call sites in the driver — https://github.com/ObjectVision/NetworkModel_EU/blob/e25416b/lambda_sweep_simplex.jl#L278

Knobs

Only the knobs read inside lp_run.jl; the travel function, BIG, the population column and the MS_* rounding effort live in settings.jl (Code walkthrough settings, From LP fractions to real pharmacies).

Name Where read Default What it changes Why you would change it Interactions
SOLVER solve_at_w!, every call (line 484) simplex HiGHS solver; ipm also turns run_crossover on Only for the very largest LINEAR areas where a grid step falls off the simplex cliff (Belgium/FRI, Poland country-level) Never for LOGISTIC under soft coverage (Pitfalls). IPM ignores the warm start, so every point costs the full IPM price. resolve_for_basis! inherits it.
LP_TIME_LIMIT const at include time (line 439), applied in solve_at_w! (line 487) 3600 s HiGHS time_limit per solve Raise for Poland-sized LPs (issue #52 mentions 4 h); lower to fail fast. Set to 4 h in 86476bf (7 July 2026), 1 h in f93aea1 (9 July) so that "any pathological point is a skipped point, never a stall". A time-out throws → the point is skipped, but the aborted basis poisons the next warm start (Pitfalls); once S1/S2 are pinned the driver ends the sweep at the first failure in the tail (grid points above the brackets). Not applied by run_scenario.
rounding keyword solve_at_w! true false skips all three roundings Set by the driver's refine-only walk; not by hand LP-only rows have cost_c = Inf, no Arrow output, and are ignored by scenario selection.
ROUNDING settings.jl, used at line 543 multistart which of the three sets is reported Comparison experiments only Any other string silently means top-p. S2 is bisected on the selected set's cost_c.

Pitfalls and open questions

  • The certified bound is on total cost, not on travel. n_open = round(sum_y) differs from sum_y by up to 0.5, so the travel-only "multi vs relax" column can be off by up to 0.5·λ; it is negative in 147 of the 1 406 rows in doc/deck_data.json at e25416b (every one a point rounded up), 129 of them by less than 0.1 %, the largest by 3.65 % (Luxembourg LOGISTIC w = 0.2). Mechanism: From LP fractions to real pharmacies.
  • Pre-24-July material has x and y swapped — WarmStartState{TX,TY}, the ROUNDING comment in settings.jl, randomized_seed's comment, issue #47 and the sum_x / frac_x columns of old logs; the reading rule is on Log lines and output files.
  • A time-out poisons the warm start (issue #52): the next point starts from an aborted basis and tends to time out too — the reason for in-loop bisection, resolve_for_basis! and the tail stop. Before S1/S2 are pinned a cascade can still burn hours.
  • greedy_round prices an unreached client at Inf, not BIG (lines 31, 41–42), so coverage comes lexicographically first for greedy; intent unrecorded. A candidate reaching an unreached zero-population client (the exclusion rule zeroes population, lambda_sweep_simplex.jl around line 124) gets gain Inf × 0 = NaN, which never wins a comparison, so the greedy pass skips it. Greedy is worse than top-p in 129 of the 1 406 logged rows; whether this is the cause is not recorded (inference). Multistart repairs it either way.
  • SOLVER=ipm on LOGISTIC fails outright (OTHER_ERROR on every point, commit 4ec47a9); the commit's coefficient-range diagnosis could not be checked against a log (unverified).
  • Unused leftovers: randomized_seed's must_open; multistart_round's t_ij_col; assign_nearest's second return; the x ≤ 1 bound; nearest_facility in CountryData; the two assigned_k types.
  • Commits after e25416b cited in issues #52–#54 are not on GitHub as of 15 September 2026; solve_at_w! on Maarten's machine may differ from this page (Branches data and environment).
  • Roadmap: an exact soft-coverage MIP to pin S1/S2 instead of rounding (doc/todo.md B5); Poland country-level LINEAR needs the IPM path or the MIP (issue #52). Decision log.

See also

Clone this wiki locally