fluxopt.math

fluxopt's math, and the engine that builds it.

Three layers: Elements -> ModelData -> this. The math is :data:~fluxopt.math.sources.PROGRAM, a directory of YAML fragments that mathspec composes into one spec, built and solved by specsolve <https://github.com/fluxopt/lpspec>_; :mod:~fluxopt.math.sources binds a ModelData to it, and :mod:~fluxopt.math.results reads the answer back as a :class:~fluxopt.results.Result.

There is no second implementation. model.py built the same math a second time in linopy calls, and deleting it is what docs/design/lpspec-direction.md was written to argue for: the math is a file now, so it is reviewed, diffed, typeset and extended as one.

specsolve and mathspec are pinned to commits — their language surface is pre-1.0 and still moves, so an unpinned ref would let a uv sync change what the program means.

One refusal is left, and it is a difference of formulation rather than a gap: PiecewiseConversion.method='lp' is linopy's tangent-line relaxation, and this program has only the exact formulation, so it would answer a different question. :class:UnsupportedFeatureError says so rather than quietly answering it.

Modules:

Name	Description
parameters	The numbers a system binds, as an artifact.
results	Turn an lpspec result into a :class:~fluxopt.results.Result.
solve	Run a :class:~fluxopt.model_data.ModelData through the streaming engine.
sources	Bind a :class:~fluxopt.model_data.ModelData to fluxopt's math program.

Classes:

Name	Description
Parameters	Every table a solve binds, in the three kinds the program declares.
UnsupportedFeatureError	The ModelData uses a feature the program does not express yet.

Functions:

Name	Description
to_result	Build a :class:~fluxopt.results.Result from a solved lpspec model.
build_sources	Emit the parameter tables and coordinates for :data:PROGRAM.
program	fluxopt's math, composed from its fragments, loaded and checked.

Parameters dataclass

Parameters(
    dimensions: dict[str, DataFrame],
    lookups: dict[str, DataFrame],
    parameters: dict[str, DataFrame],
)

Every table a solve binds, in the three kinds the program declares.

dimensions:, lookups: and parameters: are the language's own three blocks, and they are three different things: a dimension is a label set, a lookup is a map from one to another, and a parameter is a function of the dims it names. The binder hands lookups over as columns on the index table of the dimension they run over, which is what lpspec reads — but that is a transport detail, and an artifact that flattened the three into two would be keeping the transport instead of the meaning.

Methods:

Name	Description
__getitem__	One table by name, of whichever kind.
of	Bind data to a program and keep what was bound, split by kind.
save	Write the set as a directory of parquet, one file per table.
load	Read a set written by :meth:save.

Attributes:

Name	Type	Description
rows	int	Total rows across every table — the size of the bound problem.

rows property

rows: int

Total rows across every table — the size of the bound problem.

__getitem__

__getitem__(name: str) -> DataFrame

One table by name, of whichever kind.

Source code in src/fluxopt/math/parameters.py

def __getitem__(self, name: str) -> pl.DataFrame:
    """One table by name, of whichever kind."""
    for group in (self.parameters, self.lookups, self.dimensions):
        if name in group:
            return group[name]
    raise KeyError(name)

of classmethod

of(
    data: ModelData, objective: str | dict[str, float], math: Spec | None = None
) -> Parameters

Bind data to a program and keep what was bound, split by kind.

The program is what tells a lookup's source key from a parameter's, so the split is read off the declarations rather than guessed from the names.

Parameters:

Name	Type	Description	Default
data	ModelData	The model data to bind.	required
objective	str | dict[str, float]	Effect id, or effect ids mapped to objective weights — what :func:~fluxopt.math.solve.solve takes. The weights are a bound parameter like any other, so the set is not complete without knowing what is being minimised.	required
math	Spec | None	The program to bind against. Defaults to the shipped one.	None

Source code in src/fluxopt/math/parameters.py

@classmethod
def of(cls, data: ModelData, objective: str | dict[str, float], math: Model | None = None) -> Parameters:
    """Bind *data* to a program and keep what was bound, split by kind.

    The program is what tells a lookup's source key from a parameter's,
    so the split is read off the declarations rather than guessed from
    the names.

    Args:
        data: The model data to bind.
        objective: Effect id, or effect ids mapped to objective weights —
            what :func:`~fluxopt.math.solve.solve` takes. The weights are
            a bound parameter like any other, so the set is not complete
            without knowing what is being minimised.
        math: The program to bind against. Defaults to the shipped one.
    """
    from fluxopt.math.results import objective_weights
    from fluxopt.math.sources import build_sources, program

    model = math if math is not None else program()
    sources, coords = build_sources(data, objective_weights(data, objective))

    # A lookup either names the dimension its values are labels of, or
    # owns an inline label space and targets nothing. The first is a
    # relation between two dimensions and reads best keyed by both; the
    # second has only its own name to be keyed by.
    # A lookup is a source key of its own, so the three kinds arrive
    # already apart and nothing has to be taken out of an index table.
    names = set(model.relations)
    return cls(
        dimensions={k: _as_frame(v, k) for k, v in coords.items()},
        lookups={k: _as_frame(v, k) for k, v in sources.items() if k in names},
        parameters={k: _as_frame(v, k) for k, v in sources.items() if k not in names},
    )

save

save(path: str | Path) -> None

Write the set as a directory of parquet, one file per table.

Parquet rather than YAML for the reason the math is YAML: each format is asked what it is good at. A million-row table has no useful line diff and no dtype of its own in YAML, and it is the dtypes that the binder checks against the program's declarations.

Parameters:

Name	Type	Description	Default
path	str | Path	Directory to write into. Created if absent.	required

Source code in src/fluxopt/math/parameters.py

def save(self, path: str | Path) -> None:
    """Write the set as a directory of parquet, one file per table.

    Parquet rather than YAML for the reason the math is YAML: each format
    is asked what it is good at. A million-row table has no useful line
    diff and no dtype of its own in YAML, and it is the dtypes that the
    binder checks against the program's declarations.

    Args:
        path: Directory to write into. Created if absent.
    """
    root = Path(path)
    for group, frames in (
        ('dimensions', self.dimensions),
        ('lookups', self.lookups),
        ('parameters', self.parameters),
    ):
        (root / group).mkdir(parents=True, exist_ok=True)
        for name, frame in frames.items():
            frame.write_parquet(root / group / f'{name}.parquet')

load classmethod

load(path: str | Path) -> Parameters

Read a set written by :meth:save.

Parameters:

Name	Type	Description	Default
path	str | Path	The directory :meth:save wrote.	required

Raises:

Type	Description
OSError	If the directory holds no parameter set.

Source code in src/fluxopt/math/parameters.py

@classmethod
def load(cls, path: str | Path) -> Parameters:
    """Read a set written by :meth:`save`.

    Args:
        path: The directory :meth:`save` wrote.

    Raises:
        OSError: If the directory holds no parameter set.
    """
    root = Path(path)
    if not (root / 'parameters').is_dir():
        raise OSError(f'No fluxopt parameter set found in {root} (missing parameters/)')

    def read(group: str) -> dict[str, pl.DataFrame]:
        return {f.stem: pl.read_parquet(f) for f in sorted((root / group).glob('*.parquet'))}

    return cls(dimensions=read('dimensions'), lookups=read('lookups'), parameters=read('parameters'))

UnsupportedFeatureError

Bases: RuntimeError

The ModelData uses a feature the program does not express yet.

Raised rather than silently dropping the feature: a missing constraint would still solve, just to the wrong answer.

to_result

to_result(
    result: Result,
    data: ModelData,
    weights: dict[str, float],
    model: Spec | None = None,
) -> Result

Build a :class:~fluxopt.results.Result from a solved lpspec model.

Parameters:

Name	Type	Description	Default
result	Result	The lpspec result — primal frames and named expressions.	required
data	ModelData	The model data that was bound. Carries the labels the program indexes by ordinal, and the tables Result reads alongside the solution.	required
weights	dict[str, float]	Effect ids mapped to their objective weight, for provenance.	required
model	Spec | None	The program that was solved, whose expressions: are read back. Defaults to the shipped one.	None

Returns:

Type	Description
Result	The same object the linopy lane returns.

Raises:

Type	Description
NoSolutionError	If the solve did not produce a primal solution — there is no result to read, and an empty one would read as zeros. lpspec's own, so a caller sweeping scenarios catches one class for every lane and records the outcome.

Source code in src/fluxopt/math/results.py

def to_result(result: Solved, data: ModelData, weights: dict[str, float], model: Model | None = None) -> Result:
    """Build a :class:`~fluxopt.results.Result` from a solved lpspec model.

    Args:
        result: The lpspec result — primal frames and named expressions.
        data: The model data that was bound. Carries the labels the program
            indexes by ordinal, and the tables ``Result`` reads alongside the
            solution.
        weights: Effect ids mapped to their objective weight, for provenance.
        model: The program that was solved, whose ``expressions:`` are read
            back. Defaults to the shipped one.

    Returns:
        The same object the linopy lane returns.

    Raises:
        NoSolutionError: If the solve did not produce a primal solution —
            there is no result to read, and an empty one would read as zeros.
            lpspec's own, so a caller sweeping scenarios catches one class for
            every lane and records the outcome.
    """
    if model is None:
        from fluxopt.math.sources import program

        model = program()
    if not result.has_primal:
        raise NoSolutionError(f'no primal solution to read: the solver terminated {result.termination_condition!r}')

    solution: dict[str, xr.DataArray] = {}
    for program_name, var_name in _DIRECT.items():
        if (arr := _read(result, program_name, data)) is not None:
            solution[var_name] = arr
    for program_name, (flow_var, component_var) in _STATUS.items():
        if (arr := _read(result, program_name, data)) is None:
            continue
        on_flows, on_components = _split_status(arr, data)
        if on_flows is not None:
            solution[flow_var] = on_flows
        if on_components is not None:
            solution[component_var] = on_components

    # The lump half of the ledger is a named expression rather than a variable:
    # nothing in the model decides it on its own, so it is read back through
    # the same compiler that built the accounting row.
    expressions = _expressions(result, model, data)
    solution[Var.EFFECT_LUMP] = expressions['effect_lump']

    dataset = xr.Dataset(
        solution,
        attrs={'objective': float(result.objective), 'objective_weights': json.dumps(weights)},
    )
    return Result(solution=dataset, data=data, expressions=expressions)

build_sources

build_sources(
    data: ModelData, objective: dict[str, float]
) -> tuple[dict[str, Any], dict[str, Any]]

Emit the parameter tables and coordinates for :data:PROGRAM.

Parameters:

Name	Type	Description	Default
data	ModelData	The model data to bind. Both backends read the same object.	required
objective	dict[str, float]	Effect ids mapped to their objective weight, as :func:~fluxopt.math.solve.solve takes it.	required

Returns:

Type	Description
tuple[dict[str, Any], dict[str, Any]]	(sources, coords) ready to pass to lpspec.solve.

Raises:

Type	Description
UnsupportedFeatureError	If data uses a feature the program does not express yet, rather than dropping it silently.

Source code in src/fluxopt/math/sources.py

def build_sources(data: ModelData, objective: dict[str, float]) -> tuple[dict[str, Any], dict[str, Any]]:
    """Emit the parameter tables and coordinates for :data:`PROGRAM`.

    Args:
        data: The model data to bind. Both backends read the same object.
        objective: Effect ids mapped to their objective weight, as
            :func:`~fluxopt.math.solve.solve` takes it.

    Returns:
        ``(sources, coords)`` ready to pass to ``lpspec.solve``.

    Raises:
        UnsupportedFeatureError: If *data* uses a feature the program does
            not express yet, rather than dropping it silently.
    """
    _reject_unsupported(data)
    reject_varying_contribution_into_lump(data)
    fds, dims = data.flows, data.dims

    # Keyed by the numpy labels an array coordinate carries, not by the frame's
    # python ones: `tidy` looks these up with whatever xarray hands it.
    time_ord = {v: i for i, v in enumerate(dims.time.values)}
    ordinals = dims.timesteps['time'].to_list()
    dt_by_time = dims.timesteps.select(['time', 'dt'])
    size_upper_of = _size_upper(data)

    def tidy(da: xr.DataArray, *, drop_zero: bool) -> pd.DataFrame:
        """`_tidy` with this model's time-ordinal mapping bound in."""
        return _tidy(da, drop_zero=drop_zero, time_ord=time_ord)

    flow_ids = fds.ids
    sources: dict[str, Any] = {}
    # Bound up front: the feature blocks below fill these only when the
    # corresponding container is present, but later blocks read them.
    sizing_ids: list[str] = []
    status_ids: list[str] = []
    # Lump-domain accumulators, filled by the flow- and storage-sizing blocks.
    # (parameter name, entity dim, coefficients) — the entity dim is carried so
    # an absent term still emits a correctly keyed empty table.
    lump_terms: list[tuple[str, str, xr.DataArray | None]] = []
    #: (parameter, entity dim, frame, value column) — the frame-backed ones
    lump_frames: list[tuple[str, str, pl.DataFrame, str]] = []
    #: (parameter, frame, value column) — charged where the build happened
    at_build_frames: list[tuple[str, pl.DataFrame, str]] = []

    # --- flow rate bounds: dense, they sit at the variable's own grid -----
    # The envelope carries the relative pair, the profile and the fixed size
    # join onto it, and the bound is then one chain of overrides read top to
    # bottom — the same order the array version applied them in.
    envelope = fds.envelope
    profile = fds.fixed_profile.rename({'value': 'fixed'})
    grid = envelope.join(fds.sizes, on='flow', how='left').join(profile, on=['flow', 'time', 'period'], how='left')

    sz, inv, st = fds.sizing, fds.invest, data.status
    sizing_ids = sz.ids if sz is not None else []
    invest_ids = inv.ids if inv is not None else []
    # The status axis holds both kinds; these are the flows whose own rate
    # bounds the binary takes over, so components have no business here.
    status_ids = [e for e in st.ids if e in set(flow_ids)] if st is not None else []
    is_bounded, is_profile = fds.bounded_ids, fds.profiled_ids

    def among(ids: list[str]) -> pl.Expr:
        """Rows whose flow is one of *ids*."""
        return pl.col('flow').is_in(pl.Series(ids, dtype=pl.String).implode())

    scaled_max = pl.col('size') * pl.when(among(is_profile)).then(pl.col('fixed')).otherwise(
        pl.col('relative_rate_max')
    )
    bounds = grid.with_columns(
        pl.when(among(is_bounded)).then(pl.col('size') * pl.col('relative_rate_min')).otherwise(0.0).alias('rate_min'),
        pl.when(among(is_bounded))
        .then(pl.col('size') * pl.col('relative_rate_max'))
        .otherwise(np.inf)
        .alias('rate_max'),
    )
    # A fixed profile pins the rate: both bounds land on the same value.
    pinned = among(is_profile) & pl.col('fixed').is_not_null()
    bounds = bounds.with_columns(
        pl.when(pinned).then(pl.col('size') * pl.col('fixed')).otherwise(pl.col('rate_min')).alias('rate_min'),
        pl.when(pinned).then(pl.col('size') * pl.col('fixed')).otherwise(pl.col('rate_max')).alias('rate_max'),
    )
    # An optimized size is a variable, so the rate is free here and the
    # envelope is applied against the size variable instead.
    free = among([*sizing_ids, *invest_ids])
    bounds = bounds.with_columns(
        pl.when(free).then(0.0).otherwise(pl.col('rate_min')).alias('rate_min'),
        pl.when(free).then(np.inf).otherwise(pl.col('rate_max')).alias('rate_max'),
    )
    # `on` carries the envelope for status flows; the variable itself is free
    # above 0, except when it is sized too and the size variable takes over.
    gated = among(status_ids)
    bounds = bounds.with_columns(
        pl.when(gated).then(0.0).otherwise(pl.col('rate_min')).alias('rate_min'),
        pl.when(gated & among(sizing_ids))
        .then(np.inf)
        .when(gated)
        .then(scaled_max)
        .otherwise(pl.col('rate_max'))
        .alias('rate_max'),
    )
    for key in ('rate_min', 'rate_max'):
        sources[key] = bounds.select(['flow', 'time', 'period', pl.col(key).alias('value')])
    # The big-M for every binary that has to release a rate — a ramp across a
    # start-up, a storage's charge/discharge exclusion. Stated once on `flow`;
    # the storage side reads it through `charge_storage` rather than keeping a
    # copy on its own axis.
    sources['size_bound'] = pd.DataFrame({'flow': flow_ids, 'value': [size_upper_of[f] for f in flow_ids]})

    # --- carrier balance --------------------------------------------------
    membership = data.carriers.membership
    flow_index = pd.DataFrame({'flow': flow_ids, 'carrier_of': membership['carrier'].to_list()})
    sources['carrier_sign'] = pd.DataFrame({'flow': flow_ids, 'value': membership['sign'].to_numpy()})

    # --- converters -------------------------------------------------------
    if data.converters is not None:
        cds = data.converters
        coeffs = cds.coefficients
        conv_of = dict(zip(coeffs['flow'], coeffs['converter'], strict=True))
        flow_index['converter_of'] = [conv_of.get(f) for f in flow_ids]
        # Already the table the parameter wants — only the time labels are
        # the element layer's and the program indexes them by ordinal.
        sources['conversion_factor'] = _with_time_ordinals(coeffs.filter(pl.col('value') != 0), dims).select(
            ['flow', 'eq_idx', 'time', 'value']
        )
        # One row per equation each converter states — the counts, expanded.
        sources['conversion_active'] = pl.DataFrame(
            {
                'converter': [c for c, n in zip(cds.ids, cds.equations['n_equations'], strict=True) for _ in range(n)],
                'eq_idx': [i for n in cds.equations['n_equations'] for i in range(n)],
            }
        ).with_columns(pl.lit(True).alias('value'))
    else:
        flow_index['converter_of'] = None
        sources['conversion_factor'] = _empty('conversion_factor', 'flow', 'eq_idx', 'time')
        sources['conversion_active'] = _empty('conversion_active', 'converter', 'eq_idx')

    # --- storage ----------------------------------------------------------
    storage_ids: list[str] = []
    if data.storages is not None:
        sds = data.storages
        storage_ids = sds.ids
        for lookup, column in (('charge_storage', 'charge_flow'), ('discharge_storage', 'discharge_flow')):
            of_flow = dict(zip(sds.storages[column], storage_ids, strict=True))
            flow_index[lookup] = [of_flow.get(f) for f in flow_ids]

        # One join carries every per-timestep storage parameter, since they
        # all live on the same (storage, time) rows.
        profiles = _with_time_ordinals(sds.profiles, dims).join(dt_by_time, on='time')

        gains = profiles.with_columns(
            (pl.col('eta_charge') * pl.col('dt')).alias('charge_gain'),
            (pl.col('dt') / pl.col('eta_discharge')).alias('discharge_draw'),
            ((1 - pl.col('loss')) ** pl.col('dt')).alias('retention'),
        )
        for key in ('charge_gain', 'discharge_draw', 'retention'):
            sources[key] = gains.select(['storage', 'time', pl.col(key).alias('value')])

        # An absent capacity row is a storage whose capacity is a variable, so
        # its absolute level bounds are not knowable here: 0 and infinity are
        # what the program reads while the relative pair does the bounding.
        absolute = profiles.join(sds.capacity, on='storage', how='left')
        sources['level_min'] = absolute.select(
            ['storage', 'time', (pl.col('relative_level_min') * pl.col('capacity')).fill_null(0.0).alias('value')]
        )
        sources['level_max'] = absolute.select(
            ['storage', 'time', (pl.col('relative_level_max') * pl.col('capacity')).fill_null(np.inf).alias('value')]
        )
        csz = sds.sizing
        if csz is not None:
            cap_ids = csz.ids
            sources['has_capacity_sizing'] = _flags('has_capacity_sizing', 'storage', cap_ids)
            sources['capacity_mandatory'] = _flags(
                'capacity_mandatory', 'storage', csz.bounds.filter('mandatory')['entity'].to_list()
            )
            sources['capacity_min'] = pd.DataFrame({'storage': cap_ids, 'value': csz.bounds['size_min'].to_numpy()})
            sources['capacity_max'] = pd.DataFrame({'storage': cap_ids, 'value': csz.bounds['size_max'].to_numpy()})
            sized = profiles.filter(pl.col('storage').is_in(pl.Series(cap_ids).implode()))
            for key, column in (
                ('relative_level_min', 'relative_level_min'),
                ('relative_level_max', 'relative_level_max'),
            ):
                sources[key] = sized.select(['storage', 'time', pl.col(column).alias('value')]).filter(
                    pl.col('value') != 0
                )
            lump_frames += [
                ('effects_per_capacity', 'storage', csz.effects, 'per_size'),
                ('effects_fixed_capacity', 'storage', csz.effects, 'fixed'),
            ]
        for key, column in (('is_cyclic', 'cyclic'), ('prevent_simultaneous', 'prevent_simultaneous')):
            sources[key] = _flags(key, 'storage', sds.storages.filter(pl.col(column))['storage'].to_list())
        # `prior_level` is dense: it sits on the constant side of the first
        # step's balance, where an absent row is a binding zero — which is
        # also what an unset prior level means.
        sources['prior_level'] = (
            sds.storages.select('storage')
            .join(sds.levels.select(['storage', pl.col('prior_level').alias('value')]), on='storage', how='left')
            .with_columns(pl.col('value').fill_null(0.0))
        )
        for key in ('final_level_min', 'final_level_max'):
            sources[key] = sds.levels.select(['storage', pl.col(key).alias('value')]).drop_nulls()
    else:
        flow_index['charge_storage'] = None
        flow_index['discharge_storage'] = None
        for name, dcols in (
            ('is_cyclic', ['storage']),
            ('prior_level', ['storage']),
            ('final_level_min', ['storage']),
            ('final_level_max', ['storage']),
            ('prevent_simultaneous', ['storage']),
        ):
            sources[name] = pd.DataFrame({c: [] for c in [*dcols, 'value']})
        for name, dcols in (
            ('charge_gain', ['storage', 'time']),
            ('discharge_draw', ['storage', 'time']),
            ('retention', ['storage', 'time']),
            ('level_min', ['storage', 'time']),
            ('level_max', ['storage', 'time']),
        ):
            sources[name] = pd.DataFrame({c: [] for c in [*dcols, 'value']})

    # --- status -----------------------------------------------------------
    # One table, one axis. An entity here is a flow carrying its own on/off
    # decision or a component whose decision governs several flows; they carry
    # the same fields and obey the same math, so the program states the family
    # once over `status_entity`. Which rows read which binary is `status_of`,
    # and nothing else distinguishes them.
    ec_extra: list[tuple[str, xr.DataArray]] = []
    #: (parameter, frame, value column) — status coefficients, per timestep
    status_effect_frames: list[tuple[str, pl.DataFrame, str]] = []
    #: flow id -> the entity whose binary gates it. A self-status flow maps to
    #: itself; a governed flow to its component; an ungated flow to nothing.
    status_of: dict[str, str] = {f: f for f in status_ids}
    # A piecewise curve's flows are gated by its convexity row, which already
    # pins every weight to zero when the binary is off. Gating them a second
    # time per flow would be redundant, and wrong for the links the curve
    # only bounds.
    pw_comps = set(data.piecewise.converter_ids()) if data.piecewise is not None else set()
    status_of.update(
        {
            fid: owner
            for fid, owner in zip(fds.governed_by['flow'], fds.governed_by['component'], strict=True)
            if owner not in pw_comps
        }
    )

    flow_index['status_of'] = [status_of.get(f) for f in flow_ids]
    entity_ids: list[str] = st.ids if st is not None else []
    gated_ids = [f for f in flow_ids if f in status_of]

    sources['is_gated'] = _flags('is_gated', 'flow', gated_ids)
    sources['is_bounded'] = _flags('is_bounded', 'flow', is_bounded)
    sources['is_profile'] = _flags('is_profile', 'flow', is_profile)

    if st is not None:
        durations, prior, status_effects = st.durations, st.prior, st.effects
        all_entities = pl.DataFrame({'status_entity': entity_ids}, schema={'status_entity': pl.String})

        def per_entity(frame: pl.DataFrame, column: str) -> Any:
            """One column of a frame, keyed on the entity axis, live rows only."""
            live = frame.filter(pl.col(column).is_not_null()).select(
                [pl.col('entity').alias('status_entity'), pl.col(column).alias('value')]
            )
            return live if len(live) else _empty(column, 'status_entity')

        horizon = float(dims.timesteps['dt'].sum())
        bounded_up = durations.filter(pl.col('uptime_min').is_not_null() | pl.col('uptime_max').is_not_null())[
            'entity'
        ].to_list()
        bounded_down = durations.filter(pl.col('downtime_min').is_not_null() | pl.col('downtime_max').is_not_null())[
            'entity'
        ].to_list()
        sources['has_uptime'] = _flags('has_uptime', 'status_entity', bounded_up)
        sources['has_downtime'] = _flags('has_downtime', 'status_entity', bounded_down)
        sources['uptime_min'] = per_entity(durations, 'uptime_min')
        sources['downtime_min'] = per_entity(durations, 'downtime_min')
        sources['initial_status'] = per_entity(prior, 'initial')
        sources['previous_uptime'] = per_entity(prior, 'previous_uptime')
        sources['previous_downtime'] = per_entity(prior, 'previous_downtime')

        # Big-M covers the whole horizon plus whatever ran before it. An
        # entity with no prior has no row in `prior`, and no prior is zero.
        big_m = all_entities.join(
            prior.select([pl.col('entity').alias('status_entity'), 'previous_uptime', 'previous_downtime']),
            on='status_entity',
            how='left',
        ).with_columns(
            (pl.col('previous_uptime').fill_null(0.0) + horizon).alias('up'),
            (pl.col('previous_downtime').fill_null(0.0) + horizon).alias('down'),
        )
        sources['uptime_big_m'] = big_m.select(['status_entity', pl.col('up').alias('value')])
        sources['downtime_big_m'] = big_m.select(['status_entity', pl.col('down').alias('value')])

        # The duration variables' ceiling: the declared maximum where there is
        # one, the big-M where there is not. Per entity — lpspec broadcasts it
        # over the (time, period) grid the variable is declared on.
        for kind, ids, declared, fallback in (
            ('uptime', bounded_up, 'uptime_max', 'up'),
            ('downtime', bounded_down, 'downtime_max', 'down'),
        ):
            ceiling = (
                pl.DataFrame({'status_entity': ids}, schema={'status_entity': pl.String})
                .join(durations.select([pl.col('entity').alias('status_entity'), declared]), on='status_entity')
                .join(big_m.select(['status_entity', fallback]), on='status_entity')
                .with_columns(pl.col(declared).fill_null(pl.col(fallback)).alias('value'))
                .select(['status_entity', 'value'])
            )
            sources[f'{kind}_upper'] = ceiling

        # A prior run shorter than the minimum forces continuation.
        for forced_key, prev_column, min_column in (
            ('forced_on_at_start', 'previous_uptime', 'uptime_min'),
            ('forced_off_at_start', 'previous_downtime', 'downtime_min'),
        ):
            forced = prior.join(durations.select(['entity', min_column]), on='entity', how='inner').filter(
                (pl.col(prev_column) > 0) & (pl.col(prev_column) < pl.col(min_column))
            )
            sources[forced_key] = _flags(forced_key, 'status_entity', forced['entity'].to_list())

        # `dt` turns a per-running-hour rate into the step's cost; a startup
        # happens once at the step, so it is not scaled.
        scaled = (
            _with_time_ordinals(status_effects, dims)
            .join(dt_by_time, on='time')
            .with_columns((pl.col('running') * pl.col('dt')).alias('running'))
        )
        status_effect_frames += [
            ('effects_per_running_hour', scaled, 'running'),
            ('effects_per_startup', scaled, 'startup'),
        ]
    else:
        for n in ('has_uptime', 'has_downtime'):
            sources[n] = _empty(n, 'status_entity')
        for n in (
            'uptime_min',
            'uptime_big_m',
            'uptime_upper',
            'downtime_big_m',
            'downtime_min',
            'downtime_upper',
            'initial_status',
            'previous_uptime',
            'previous_downtime',
            'forced_on_at_start',
            'forced_off_at_start',
        ):
            sources[n] = _empty(n, 'status_entity')

    sources['dt'] = dims.timesteps.select(['time', pl.col('dt').alias('value')])
    sources['is_last'] = pd.DataFrame({'time': ordinals, 'value': [i == len(ordinals) - 1 for i in ordinals]})

    # --- sizing -----------------------------------------------------------
    if sz is not None:
        bounds = sz.bounds
        sources['has_sizing'] = _flags('has_sizing', 'flow', sizing_ids)
        sources['mandatory'] = _flags('mandatory', 'flow', bounds.filter('mandatory')['entity'].to_list())
        sources['size_min'] = pd.DataFrame({'flow': sizing_ids, 'value': bounds['size_min'].to_numpy()})
        sources['size_max'] = pd.DataFrame({'flow': sizing_ids, 'value': bounds['size_max'].to_numpy()})
        lump_frames += [
            ('effects_per_size', 'flow', sz.effects, 'per_size'),
            ('effects_fixed', 'flow', sz.effects, 'fixed'),
        ]
    else:
        for n in ('has_sizing', 'mandatory', 'size_min', 'size_max'):
            sources[n] = _empty(n, 'flow')

    if 'has_capacity_sizing' not in sources:
        for n in ('has_capacity_sizing', 'capacity_mandatory', 'capacity_min', 'capacity_max'):
            sources[n] = _empty(n, 'storage')
        for n in ('relative_level_min', 'relative_level_max'):
            sources[n] = _empty(n, 'storage', 'time')

    # --- ramps ------------------------------------------------------------
    # A ramp limit is per hour, so the step's allowance is limit x dt. Where
    # the size is a number that allowance is absolute; where it is a variable
    # the per-unit coefficient travels instead and the program multiplies.
    ramps = fds.ramps.join(dt_by_time, on='time').join(fds.sizes, on='flow', how='left')
    for kind in ('up', 'down'):
        declared = ramps.filter(pl.col(f'ramp_{kind}').is_not_null())
        sources[f'has_ramp_{kind}'] = _flags(
            f'has_ramp_{kind}', 'flow', declared['flow'].unique(maintain_order=True).to_list()
        )
        allowance = pl.col(f'ramp_{kind}') * pl.col('dt')
        # The split follows the program's `has_sizing`, not the absence of a
        # fixed size: those are the two branches the constraints are written in.
        sized = pl.col('flow').is_in(pl.Series(sizing_ids, dtype=pl.String).implode())
        sources[f'ramp_{kind}_coeff'] = _live(declared.filter(sized), allowance)
        sources[f'ramp_{kind}_limit'] = _live(declared.filter(~sized), allowance * pl.col('size'))
    # --- investment -------------------------------------------------------
    if inv is not None:
        period_labels_inv: list[Any] = dims.periods['label'].to_list()
        n_p = len(period_labels_inv)
        # No row means forever, so the lookup's default is the absence.
        expires = dict(zip(inv.lifetime['entity'], inv.lifetime['periods'], strict=True))
        lifetime = [expires.get(f) for f in invest_ids]
        prior = inv.bounds['prior_size'].to_numpy()
        window = np.zeros((len(invest_ids), n_p, n_p))
        prior_active = np.zeros((len(invest_ids), n_p))
        for f_idx in range(len(invest_ids)):
            lt_int = lifetime[f_idx]
            for p_idx in range(n_p):
                for b_idx in range(n_p):
                    alive = b_idx <= p_idx if lt_int is None else b_idx <= p_idx < b_idx + lt_int
                    window[f_idx, p_idx, b_idx] = float(alive)
                if prior[f_idx] > 0 and (lt_int is None or p_idx < lt_int):
                    prior_active[f_idx, p_idx] = 1.0
        coords_w = {'flow': invest_ids, 'period': period_labels_inv, 'build_period': period_labels_inv}
        sources['lifetime_window'] = _tidy(
            xr.DataArray(window, dims=['flow', 'period', 'build_period'], coords=coords_w), drop_zero=True
        )
        sources['prior_capacity_active'] = _tidy(
            xr.DataArray(
                prior_active, dims=['flow', 'period'], coords={'flow': invest_ids, 'period': period_labels_inv}
            ),
            drop_zero=False,
        )
        sources['has_invest'] = _flags('has_invest', 'flow', invest_ids)
        # The same parameter the sizing block fills. A flow's size is a
        # `Sizing` or an `Investment`, never both, so the two never collide.
        sources['mandatory'] = pd.concat(
            [sources['mandatory'], pd.DataFrame({'flow': invest_ids, 'value': inv.bounds['mandatory'].to_numpy()})],
            ignore_index=True,
        )
        # The same two parameters the sizing block fills, over the other half
        # of the sized flows.
        for key in ('size_min', 'size_max'):
            sources[key] = pd.concat(
                [sources[key], pd.DataFrame({'flow': invest_ids, 'value': inv.bounds[key].to_numpy()})],
                ignore_index=True,
            )
        sources['prior_capacity'] = pd.DataFrame({'flow': invest_ids, 'value': prior})
        sources['has_prior_capacity'] = _flags(
            'has_prior_capacity', 'flow', [f for f, ps in zip(invest_ids, prior, strict=True) if ps > 0]
        )
        lump_frames += [
            ('effects_per_size_recurring', 'flow', inv.effects, 'per_size_recurring'),
            ('effects_fixed_recurring', 'flow', inv.effects, 'fixed_recurring'),
        ]
        # A one-time cost belongs to the period the build happened in, not to
        # every period the unit is alive — which as rows is the diagonal: the
        # same period twice, rather than a square matrix multiplied by an eye.
        at_build_frames += [
            ('effects_per_size_at_build', inv.effects, 'per_size_at_build'),
            ('effects_fixed_at_build', inv.effects, 'fixed_at_build'),
        ]
    else:
        for name in (
            'has_invest',
            'prior_capacity',
            'has_prior_capacity',
        ):
            sources[name] = _empty(name, 'flow')
        sources['lifetime_window'] = _empty('lifetime_window', 'flow', 'period', 'build_period')
        sources['prior_capacity_active'] = _empty('prior_capacity_active', 'flow', 'period')

    # The envelope is read by both mechanisms and by the status family, so it
    # must span every flow whose rate is scaled by a size: sized flows, whose
    # `size` is a variable, and gated flows, whose `size_bound` is a number.
    sized_ids = [*sizing_ids, *invest_ids]
    scaled = grid.filter(pl.col('flow').is_in(pl.Series([*sized_ids, *gated_ids], dtype=pl.String).implode()))
    sources['relative_rate_max'] = _live(scaled, pl.col('relative_rate_max'))
    sources['fixed_relative_profile'] = _live(scaled, pl.col('fixed'))
    # Dense: the lower bound also stands on the constant side of
    # `status_sizing_rate_min`, where a dropped zero is a bound rather than an
    # absent coefficient — and a flow whose lower bound is zero is the
    # ordinary case, so dropping it would break exactly the common one.
    sources['relative_rate_min'] = _live(scaled, pl.col('relative_rate_min'), drop_zero=False)

    # --- effects: the sparse one -----------------------------------------
    eds = data.effects
    effect_ids = eds.ids
    # dt stays: a per-flow-hour rate times a duration is the step's energy.
    # The aggregation weight does not — the program applies it in the sum,
    # so a named contribution reads as the physical per-step quantity.
    cf = eds.cf_matrix()
    leo = leontief(cf) if cf is not None else None
    pairs = fds.effect_pairs.join(dt_by_time, on='time').select(
        ['flow', 'effect', 'time', 'period', (pl.col('value') * pl.col('dt')).alias('value')]
    )
    sources['effects_per_flow_hour'] = _fold_effect_rows(pairs.rename({'flow': 'entity'}), 'flow', 'value', leo, tidy)
    for name, arr in ec_extra:
        sources[name] = tidy(apply_leontief(leo, arr) if leo is not None else arr, drop_zero=True)
    for name, frame, column in status_effect_frames:
        sources[name] = _fold_effect_rows(frame, 'status_entity', column, leo, tidy)
    for name in ('effects_per_running_hour', 'effects_per_startup'):
        sources.setdefault(name, _empty(name, 'status_entity', 'effect', 'time', 'period'))

    # Lump domain: effect_lump = (I - cf_lump)^-1 . lump_direct, folded into
    # the coefficients so no self-referential effect_lump variable is needed.
    leo_lump = leontief(cf.mean('time')) if cf is not None else None

    def fold(arr: xr.DataArray) -> xr.DataArray:
        """Apply the lump-domain Leontief inverse, if there are cross-effects."""
        return arr if leo_lump is None else apply_leontief(leo_lump, arr)

    for name, entity_dim, arr in lump_terms:
        sources[name] = tidy(fold(arr), drop_zero=True) if arr is not None else _empty(name, entity_dim, 'effect')
    for name, entity_dim, frame, column in lump_frames:
        sources[name] = _fold_lump_rows(frame, entity_dim, column, leo_lump)
    for name, frame, column in at_build_frames:
        folded = _fold_lump_rows(frame, 'flow', column, leo_lump)
        sources[name] = folded.with_columns(pl.col('period').alias('build_period'))
    for name in ('effects_per_size', 'effects_fixed'):
        sources.setdefault(name, _empty(name, 'flow', 'effect', 'period'))
    for name in ('effects_per_capacity', 'effects_fixed_capacity'):
        sources.setdefault(name, _empty(name, 'storage', 'effect', 'period'))
    for name in ('effects_per_size_at_build', 'effects_fixed_at_build'):
        sources.setdefault(name, _empty(name, 'flow', 'effect', 'period', 'build_period'))
    for name in ('effects_per_size_recurring', 'effects_fixed_recurring'):
        sources.setdefault(name, _empty(name, 'flow', 'effect', 'period'))
    # Objective weight x period weight, folded into one parameter. Both are
    # per (effect, period), so the fold is a join and the defaults are what a
    # missing row means: no override, then no global weight, then 1.
    global_weights = dims.periods.select(['period', pl.col('weight').alias('global_weight')])
    period_axis = list(range(dims.n_periods))
    grid = (
        pl.DataFrame({'effect': effect_ids}, schema={'effect': pl.String})
        .join(pl.DataFrame({'period': period_axis}, schema={'period': pl.Int64}), how='cross')
        .join(eds.period_weights, on=['effect', 'period'], how='left')
        .join(global_weights, on='period', how='left')
        .with_columns(pl.col('weight').fill_null(pl.col('global_weight')).fill_null(1.0).alias('weight'))
    )
    objective_by_effect = pl.DataFrame(
        {'effect': effect_ids, 'objective': [float(objective.get(e, 0.0)) for e in effect_ids]},
        schema={'effect': pl.String, 'objective': pl.Float64},
    )
    sources['objective_weight'] = (
        grid.join(objective_by_effect, on='effect')
        .with_columns((pl.col('objective') * pl.col('weight')).alias('value'))
        .filter(pl.col('value') != 0)
        .select(['effect', 'period', 'value'])
    )
    # Weights for the across-period sum: per-effect override, else global, else 1.
    sources['period_weight'] = grid.select(['effect', 'period', pl.col('weight').alias('value')])

    # --- effect limits ----------------------------------------------------
    for key, frame, axes in (
        ('periodic_min', eds.periodic, ['effect', 'period']),
        ('periodic_max', eds.periodic, ['effect', 'period']),
        ('total_min', eds.totals, ['effect']),
        ('total_max', eds.totals, ['effect']),
    ):
        sources[key] = frame.filter(pl.col(key).is_not_null()).select([*axes, pl.col(key).alias('value')])

    # --- temporal boundary mask ------------------------------------------
    sources['is_first'] = pd.DataFrame({'time': ordinals, 'value': [i == 0 for i in ordinals]})
    sources['time_weight'] = dims.timesteps.select(['time', pl.col('weight').alias('value')])

    # --- flow aggregates ------------------------------------------------
    # `size` here is the static one; a sized flow's is a variable, so its bound
    # travels as a coefficient instead of a product and the program multiplies.
    total_duration = float((dims.timesteps['dt'] * dims.timesteps['weight']).sum())
    # A load factor bounds the mean rate as a fraction of the size. Where the
    # size is a number the bound is one too; where it is a variable the bound
    # travels as a coefficient and the program multiplies.
    aggregates = fds.aggregates.join(fds.sizes, on='flow', how='left')
    for name in ('flow_hours_min', 'flow_hours_max'):
        sources[name] = aggregates.select(['flow', pl.col(name).alias('value')]).drop_nulls('value')
    for kind in ('min', 'max'):
        declared = aggregates.filter(pl.col(f'load_factor_{kind}').is_not_null())
        sources[f'load_factor_{kind}_bound'] = declared.select(
            ['flow', (pl.col(f'load_factor_{kind}') * pl.col('size') * total_duration).alias('value')]
        ).drop_nulls('value')
        sources[f'load_factor_{kind}_coeff'] = declared.filter(pl.col('size').is_null()).select(
            ['flow', (pl.col(f'load_factor_{kind}') * total_duration).alias('value')]
        )

    # --- piecewise conversion ------------------------------------------
    # The curve tables are already the shape the program wants: a link is a
    # row on `flow`, so nothing has to be reshaped into link slots.
    linear_convs = data.converters.ids if data.converters is not None else []
    bp_width = 0
    pw_status_of: dict[str, str | None] = {}
    pw = data.piecewise
    if pw is not None:
        pw_convs = pw.converter_ids()
        links, curves = pw.links, pw.curves
        # A link is a (converter, flow, bound) — the breakpoints it passes
        # through are its rows, so the identity is what remains after dropping
        # the axes a curve varies along.
        identity = links.select(['converter', 'flow', 'bound']).unique(maintain_order=True)
        pair_flow = identity['flow'].to_list()

        sources['pw_bp_value'] = _with_time_ordinals(links.filter(pl.col('value') != 0), dims).select(
            ['flow', 'bp', 'time', 'value']
        )

        # Which breakpoints a curve has. Curves of different width share one
        # `bp` axis, so the mask is what stops a weight existing past the end
        # of a narrower one.
        present = links.select(['converter', 'bp']).unique(maintain_order=True).sort(['converter', 'bp'])
        bp_width = int(present['bp'].max() or 0) + 1 if len(present) else 0  # type: ignore[arg-type]
        sources['pw_bp_present'] = present.with_columns(pl.lit(True).alias('value'))

        gated = curves.filter(pl.col('has_status'))['converter'].unique(maintain_order=True).to_list()
        sources['has_piecewise'] = _flags('has_piecewise', 'converter', pw_convs)
        sources['pw_gated'] = _flags('pw_gated', 'converter', gated)
        sources['is_piecewise'] = _flags('is_piecewise', 'flow', pair_flow)
        for name, sign in (('pw_equal', '=='), ('pw_upper', '<='), ('pw_lower', '>=')):
            sources[name] = _flags(name, 'flow', identity.filter(pl.col('bound') == sign)['flow'].to_list())

        # Availability scales the envelope of the reference link — a curve's
        # first, which is what the eager lane bounds too.
        reference = identity.group_by('converter', maintain_order=True).first()
        sources['pw_ref'] = reference.select(['flow', pl.lit(1.0).alias('value')])
        widest = (
            links.join(reference.select(['converter', 'flow']), on=['converter', 'flow'])
            .group_by(['converter', 'time'])
            .agg(pl.col('value').max().alias('widest'))
        )
        sources['pw_avail_bound'] = _with_time_ordinals(
            curves.join(widest, on=['converter', 'time']).with_columns(
                (pl.col('availability') * pl.col('widest')).alias('value')
            ),
            dims,
        ).select(['converter', 'time', 'value'])

        # A gated curve's Status is keyed by the converter's own id, so the
        # lookup maps it to itself — mapping it to None reads as 'no Status'
        # and leaves the curve ungated.
        pw_status_of = {c: c for c in gated}
        of = dict(zip(identity['flow'], identity['converter'], strict=True))
        flow_index['converter_of'] = [of.get(f) or c for f, c in zip(flow_ids, flow_index['converter_of'], strict=True)]
        converter_ids = linear_convs + [c for c in pw_convs if c not in set(linear_convs)]
    else:
        for name, dcols in (
            ('has_piecewise', ('converter',)),
            ('pw_gated', ('converter',)),
            ('is_piecewise', ('flow',)),
            ('pw_equal', ('flow',)),
            ('pw_upper', ('flow',)),
            ('pw_lower', ('flow',)),
            ('pw_ref', ('flow',)),
            ('pw_bp_present', ('converter', 'bp')),
            ('pw_bp_value', ('flow', 'bp', 'time')),
            ('pw_avail_bound', ('converter', 'time')),
        ):
            sources[name] = _empty(name, *dcols)
        converter_ids = linear_convs

    # A map is its own source key, keyed `(over, into)` and holding only the
    # labels it maps — a flow charging no storage has no row, rather than a
    # null saying so. An index carrying a column named after a lookup over it
    # is refused, so the two facts stay apart all the way down.
    def maps(index: pd.DataFrame, over: str, into: dict[str, str]) -> dict[str, pl.DataFrame]:
        """One table per lookup declared over *over*, from its index columns."""
        return {
            name: pl.DataFrame(
                {over: mapped[over].tolist(), target: mapped[name].tolist()},
                schema={over: pl.String, target: pl.String},
            )
            for name, target in into.items()
            if name in index.columns
            # pandas 3 reads a missing label as NaN, which no string column takes.
            for mapped in (index[index[name].notna()],)
        }

    flow_axis = pl.DataFrame({'flow': flow_ids}, schema={'flow': pl.String})
    lookup_tables = maps(
        flow_index,
        'flow',
        {
            'carrier_of': 'carrier',
            'converter_of': 'converter',
            'charge_storage': 'storage',
            'discharge_storage': 'storage',
            'status_of': 'status_entity',
        },
    )

    # Single-period models supply a length-1 period so one program serves both.
    period_labels = dims.period_labels
    period_ord = {v: i for i, v in enumerate(period_labels)}
    p_ordinals = list(range(len(period_labels)))

    def onto_periods(df: Any, axis: str) -> Any:
        """Put *df* on the program's period axis, whichever library holds it.

        A table already naming periods has its labels mapped to ordinals; one
        that does not is the same in every period, so it is crossed onto all
        of them. Both shapes arrive in pandas and in polars, and will keep
        doing so until the last container is converted.
        """
        if isinstance(df, pl.DataFrame):
            if axis in df.columns:
                return df.with_columns(pl.col(axis).replace_strict(period_ord, return_dtype=pl.Int64))
            return df.join(pl.DataFrame({axis: p_ordinals}, schema={axis: pl.Int64}), how='cross')
        if axis in df.columns:
            return df.assign(**{axis: [period_ord[v] for v in df[axis]]})
        return df.merge(pd.DataFrame({axis: p_ordinals}), how='cross')

    for name in PERIOD_PARAMS:
        if (df := sources.get(name)) is not None:
            sources[name] = onto_periods(df, 'period')

    for name in BUILD_PERIOD_PARAMS:
        df = sources.get(name)
        if df is not None and 'build_period' in df.columns:
            sources[name] = onto_periods(df, 'build_period')

    def labels(values: Any) -> np.ndarray:
        """A string dimension's labels, carrying their type even when empty.

        A system with no storages still declares the dimension, and its
        lookups are string maps into it. A bare ``[]`` has no dtype to infer,
        binds as a null label space and fails the join against those columns —
        so the labels travel as a string array, which says what the space
        would have held.
        """
        return np.array(list(values), dtype=str)

    def axis(dim: str, values: Any) -> pl.DataFrame:
        """A dimension's labels as a one-column frame.

        A frame rather than a sequence because a strategy slices its sources:
        `solve_over` filters every table that carries the axis it cuts on, and
        a bare list is not a table it can cut. Costing nothing to always do,
        so the shipped sources are sliceable whether or not this build is.
        """
        return pl.DataFrame({dim: list(values)}, schema={dim: pl.Int64})

    coords: dict[str, Any] = {
        'time': axis('time', ordinals),
        'period': axis('period', p_ordinals),
        'build_period': axis('build_period', p_ordinals),
        'flow': flow_axis,
        'carrier': labels(data.carriers.ids),
        # Both kinds: a converter states linear equations, a piecewise curve,
        # or one of each. The axis is the union, or a curve's own converter
        # would not be a coordinate of the dimension its rows are keyed on.
        'converter': pl.DataFrame({'converter': converter_ids}, schema={'converter': pl.String}),
        'eq_idx': axis('eq_idx', range(data.converters.width) if data.converters is not None else []),
        'storage': labels(storage_ids),
        'effect': labels(effect_ids),
        'status_entity': labels(entity_ids),
        # numpy, not a list: with no piecewise converter the width is 0 and a
        # bare `[]` has no integer type for the join to match.
        'bp': axis('bp', range(bp_width)),
    }
    lookup_tables['pw_status_of'] = pl.DataFrame(
        {'converter': list(pw_status_of), 'status_entity': [pw_status_of[c] for c in pw_status_of]},
        schema={'converter': pl.String, 'status_entity': pl.String},
    ).drop_nulls('status_entity')
    _stamp_empty_dtypes(sources)
    return {**sources, **lookup_tables}, coords

program

program() -> Any

fluxopt's math, composed from its fragments, loaded and checked.

A :class:mathspec.Spec. Each file under :data:PROGRAM states one feature and loads on its own; effects.yaml declares the two halves of the ledger as sums, and every feature adds its own term to them, so merge writes the ledger. The engine verbs come from specsolve.

Source code in src/fluxopt/math/sources.py

def program() -> Any:
    """fluxopt's math, composed from its fragments, loaded and checked.

    A :class:`mathspec.Spec`. Each file under :data:`PROGRAM` states one
    feature and loads on its own; ``effects.yaml`` declares the two halves of
    the ledger as sums, and every feature adds its own term to them, so
    ``merge`` writes the ledger. The engine verbs come from specsolve.
    """
    from mathspec import merge

    fragments = {path.stem: path for path in sorted(PROGRAM.glob('*.yaml'))}
    return merge(fragments, description='fluxopt: flows, converters and storages, and what they cost.')