perf: batch the background inventory and build the model from precomputed data

[TimoDiepers]

Why

Profiling the methanol & pig iron case study (ecoinvent 3.12 + four premise REMIND-EU vintages, four impact categories, 20 intermediate flows) showed almost all of the setup time going into LCADataProcessor: 460 s of a 613 s end-to-end run.

Two hotspots, both of them avoidable work rather than unavoidable computation.

Hotspot 1 — background inventory (~290 s)

Per intermediate flow, per impact method, _calculate_inventory_of_db did:

step	cost	problem
lca.lci(demand=...)	0.057 s	builds B @ diag(x), the full 3341 x 43648 unaggregated inventory, of which only the row sums are used
to_dataframe(cutoff=1e4)	0.749 s	0.70 s of it is annotate=True, calling bd.get_node() on ~3775 background activities per call — columns that groupby("row_code") then discards
the method loop	x4	the inventory does not depend on the LCIA method, so it was recomputed once per category

Hotspot 2 — dynamic characterization (~150 s)

_construct_characterization_tensor called characterize() once per elementary flow with a one-row dataframe. The per-call cost is not the characterization itself but create_characterization_functions_from_method, rebuilt every time:

CRF per-flow:              0.120 s each  ->  ~1000 flows = 120 s
CRF batched (196 flows):   0.13 s total

What changed

Batched inventory. All intermediate flows of a database are solved against a single factorization, and only the aggregated elementary flow vector is built:

lca.build_demand_array({activity.id: 1})
aggregated = lca.biosphere_matrix @ lca.solve_linear_system()

That is the column of B A^-1 belonging to the flow — exactly what the optimization consumes. lci() is bypassed on purpose, to_dataframe() is gone entirely, and the method loop is removed.

One-pass lookups. bd.get_node() costs 0.3–1.4 s per call; iterating a full background database costs ~1.2 s. Activities are now resolved from a per-database name index, biosphere codes and names from a single pass over the biosphere database.

Caching. Inventories are cached at module level per (project, database, modified token, activity identity, cutoff), in the style of bw_timex's background LCI cache. The modified token means editing a database invalidates its entries instead of silently reusing stale ones. A rerun in the same session takes 0.2 s. clear_lca_caches() empties it.

Parallel by default. calculation_method="parallel" was a NotImplementedError stub; it now runs one worker process per background database via stdlib ProcessPoolExecutor (no new dependency), and is the default. Building and factorizing each database's technosphere matrix is the remaining bulk of the work and is fully independent between databases. A single pending database stays in-process, so there is no pool overhead in that case. Workers inherit both the project and the Brightway base directories.

Batched characterization. One characterize() call for all flows per method, with the characterization functions built once and cached. The GWP branch also had an O(n^2) df.loc[df["flow"] == ...] lookup per result row, replaced with a map.

Changed results — the old numbers were wrong

background_inventory.cutoff defaulted to 1e4, which truncated the unaggregated inventory to its 10 000 largest entries before summing per elementary flow. The resulting flow amounts were systematically low. Verified against bw2calc's own LCIA score for the case study's intermediate flows in ei312_REMIND-EU_SSP2_NDC_2020:

intermediate flow                          bw2calc          new          old   new err   old err
market for blast furnace                7.3473e+00   7.3473e+00   5.7377e+00     0.00%   -21.91%
direct air capture system, solvent      8.3404e+00   8.3404e+00   6.0208e+00    -0.00%   -27.81%
market group for electricity, low v     4.2563e-09   4.2563e-09   3.2161e-09     0.00%   -24.44%
electrolyzer production, PEM, BoP       3.6189e-03   3.6189e-03   2.5335e-03     0.00%   -29.99%
market for tap water                    1.4545e-11   1.4545e-11   1.2133e-11     0.00%   -16.58%

The default is now None (keep every non-zero flow); a numeric cutoff still works but applies after aggregation. On the case study the unconstrained climate objective moves from 1.6890e-06 to 1.7327e-06 (+2.6%).

Because nothing is dropped by magnitude any more, the number of elementary flows per database column rises sharply (445 -> 2794 for one flow), which by itself made model building 3x slower. So flows that have no characterization factor in any configured category are now dropped: they contribute exactly zero to every impact, and the inventory is expressed per (process, elementary flow, year) in the model. Flows constrained directly — iridium in the case study — are kept with background_inventory.retain_flows, or restrict_to_characterized_flows=False keeps everything.

Results

Case study, Apple M2 Pro, 10 cores:

stage	before	after (sequential)	after (parallel, default)
LCADataProcessor	460.1 s	124.4 s	35.6 s
create_model	98.8 s	32.0 s	32.0 s
solve (Pyomo wall clock)	54.0 s	17.4 s	17.4 s
total	613 s	174 s	85 s

Rerunning the processor in the same session: 0.2 s.

Docs and notebooks

• docs/content/optimization_setup.md gains a Background Inventory section covering every setting, with warnings about dropped flows and about the if __name__ == "__main__": guard that a script needs for the parallel default (notebooks do not).
• docs/content/constraints.md warns at the flow limit section that constraining an uncharacterized flow requires retain_flows.
• notebooks/methanol_and_iron.ipynb retains iridium explicitly and explains why.
• notebooks/methanol_and_iron_v0.5.0.ipynb preserves the case study exactly as it was run for the paper with optimex 0.5.0, since the results there were produced with the truncating inventory.

Tests

tests/test_background_inventory_performance.py, six new tests: agreement with bw2calc's inventory, cache reuse across processors, parallel == sequential, cutoff behaviour, flow dropping, and both ways of keeping a flow. Full suite: 127 passed, 3 skipped.

Notes for review

• black (from the repo's pre-commit config) reformatted parts of lca_processor.py that this PR does not otherwise touch; the file was not previously black-clean. flake8 still reports pre-existing E501s in untouched docstrings (37 on main, 26 here).
• scikit-umfpack is in the working tree's pyproject.toml but not committed here. It only affects which solver backs the factorization; the batching wins do not depend on it.

---

Part 2 — building and solving the model

With the LCA processing fixed, the remaining time was in Pyomo. Profiling create_model on the paper inputs (7 processes, 20 intermediate flows, 751 elementary flows, 6 background databases, 31 years):

28,841,851 calls   61.1 s cumulative   pyomo/core/base/param.py:1029(__getitem__)
34,138,164 calls   56.0 s cumulative   pyomo/core/base/indexed_component.py:613(__getitem__)

Nearly 29 million Param lookups, each running full index validation, to read tensors that contain no decision variables at all.

What changed

Data is combined as data. The background inventory is collapsed over background databases once — sum_bkg G[bkg, i, e] * M[bkg, t] — into a plain dict keyed by (intermediate flow, elementary flow, year), instead of being re-evaluated through Params inside every expression rule.

The sum over elementary flows is folded into the characterization factors. The impact of a process was written as

sum(Q[c, e, t] * scaled_inventory[p, e, t] for e in ELEMENTARY_FLOW)

where each scaled_inventory[p, e, t] is itself a sum over intermediate flows. Since Q and the background inventory are both data, the sum over e can be done up front, giving a characterized impact per unit of intermediate flow. The expression then runs over 20 intermediate flows instead of 751 elementary flows — the same number, from a tree that is orders of magnitude smaller. That matters twice: once when Pyomo builds it, and again when the LP writer walks it.

Sparse, combined flow expressions. scaled_*_dependent_on_installation and scaled_*_dependent_on_operation were always used as a pair, and both were indexed densely over (process, flow, year) — for the elementary flows that is 163 k entries of which a handful are non-zero, since a process emits only a few flows directly. They are now a single scaled_technosphere_flow / scaled_biosphere_flow / scaled_internal_demand_flow expression indexed only over the (process, flow) pairs that carry data or a vintage override. Nothing outside optimizer.py referenced the old names.

Caching in the rules. Production rates per (process, product, vintage), operation flow rates per (process, flow, vintage) and the active vintages per (process, year) were recomputed inside rules that run tens of thousands of times, each scanning the whole ACTIVE_VINTAGE_TIME set or the whole process-time range. All three are now precomputed or memoized.

Post-processing. get_dynamic_inventory() evaluated scaled_inventory[p, e, t] — a deep expression tree — 163 k times, then called biosphere_db.get(code=...) once per row. It now multiplies solved flow values by the same background inventory data, and maps codes to ids from a single pass over the biosphere database.

Results

Paper case study (notebooks/data/paper/model_inputs_2050_paper.json plus the notebook's constraints), Apple M2 Pro:

stage	before	after
create_model	103.6 s	1.7 s
solve_model wall clock	55.7 s	0.1 s
get_dynamic_inventory()	19.8 s	0.2 s

Gurobi itself was always ~0.04 s here; the 55 s was Pyomo writing the LP file, which shrinks with the expression tree.

Smaller scenario (unconstrained, 230 flows, 4 databases): build 24.4 s → 0.9 s, solve 16.8 s → 0.1 s.

Equivalence

Both scenarios were solved before and after and compared:

• objective identical to all 12 digits (1.732682875061e-06, 2.880103392725e-06)
• all 217 var_installation and up to 2 983 var_operation values identical
• identical variable and constraint counts (3 136 / 4 410 and 3 200 / 4 501)
• get_dynamic_inventory() returns the same 45 651 non-zero rows, max relative deviation 1.5e-15

The only intentional output change: rows whose amount is exactly zero are no longer emitted by get_dynamic_inventory() (162 967 → 45 651 rows). characterize() discards them anyway.

solve_model gained no new behaviour; product_demand_fulfillment_rule now returns Constraint.Feasible when a sparse expression collapses to a constant, which previously could not happen because every expression existed densely.

Two things found on the way, not changed here

• Results are not reproducible across processes at ~1e-5 relative. The same model solved in different processes returned 1.7326829e-06 and 1.7326549e-06. Within one process it is deterministic, and Method=1, Threads=1 also pins it. The LP has alternative optima and Gurobi's default concurrent method picks between them. Worth pinning for published runs.
• appsi_gurobi cannot load the model: ProcessDeploymentLimitMax[...] has neither a lower nor an upper bound, because the limit defaults to float("inf") and the constraint is built unconditionally. Constraint.Skip for infinite limits would remove those rows (Gurobi reports 4 056 rows for a model Pyomo counts as 4 410), but it changes the reported model size, so I left it alone.

Addendum — on-disk inventory cache

With everything above in place, the LCA processing was ~97% of the runtime, essentially all of it bw2calc building and factorizing one technosphere matrix per background database. That work was cached for the running session but repeated by every new one.

Calculated inventories are now also written to optimex-inventory-cache/ inside the current Brightway project, one file per background database, keyed by the database's modified token, the activity identity and the cutoff — a stale entry misses rather than being served. Writes are atomic (temp file plus os.replace) and happen in the parent process, so parallel workers never contend for the same file. Superseded files of the same database are pruned on write.

Methanol & iron case study, whole pipeline:

run	processor	build	solve	dynamic inventory	total
first (cold cache)	37.2 s	0.9 s	0.2 s	0.1 s	38.4 s
next session (warm cache)	0.3 s	0.9 s	0.1 s	0.1 s	1.5 s

The tensor built from the cache is bit-identical to a freshly calculated one (27 600 entries, max relative deviation 0.0). Cache size is ~1.3 MB per background database.

Settings: background_inventory.use_disk_cache (default True) and disk_cache_dir. lca_processor.clear_lca_caches(include_disk=True) removes the files. Three tests cover reuse across a cleared session, invalidation after a database edit, and the disabled case.

[codecov]

Codecov Report

❌ Patch coverage is 86.95652% with 63 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/optimex/lca_processor.py	86.24%	41 Missing ⚠️
src/optimex/postprocessing.py	0.00%	17 Missing ⚠️
src/optimex/converter.py	50.00%	3 Missing ⚠️
src/optimex/optimizer.py	98.76%	2 Missing ⚠️

📢 Thoughts on this report? Let us know!