# /for-agents

# For agents

The whole manual is at [`/llms-full.txt`](/llms-full.txt); the index is at
[`/llms.txt`](/llms.txt).

## The mental model

**A variable is a dimension.** `m.var("x", (P, W))` occupies a block of the
model's single column space. A member's column is computed from its
multi-index rather than stored, so a variable over millions of columns
costs its members and not its columns.

**An expression is symbolic.** `cost[P, W] * x[P, W]` holds references,
not arrays. Writing it costs nothing. It becomes matrix entries only when a
constraint is materialised.

**A constraint is an array.** It is a `nimblend` array over its free sets
crossed with the column space, so there is no assembly step: the array is
the matrix.

**A definition is a model without its data.** `Definition` mirrors the
vocabulary a model is written in, `set`, `param`, `var`, `eq` and
`set_objective`, over symbols declared with no members and no values.
`explain()` reports what it declares; `build(data)` binds a copy and returns
a `Model`, so one definition builds as many models as it is given datasets.

**A built model answers questions about itself.** `explain()` reports what
it built, `row(name, **coords)` reads one row back out of the assembled
matrix, and `absent(name)` reports which coordinates were dropped from a
constraint and by which rule. All three read what was built rather than
walking the expression a second time.

**A session keeps the solver open.** `model.session()` assembles once and
keeps the solver's model, so `diagnose()` asks the solved instance which
rows conflict, or which direction an unbounded model runs off in.
`available()` lists the adapters installed and `capabilities(name)` reports
what each does, including what it refuses: a model with integer columns has
no duals, because a mixed-integer model's duals are not its relaxation's.

**`nimblend` is the layer below.** It knows dimensions, labels, entries and
alignment, and nothing about optimization. Import from `nimblend` itself,
never from `nimblend.sparse` or another submodule, and never read an array's
`.index` or `.data` or a domain's `.codes`. Each has a reader above it:
`coordinates()`, `values()`, `positions_of_coordinates()` and `as_coord()`.
Nor build one: a domain returns the array over its own members through
`array(values)` and `identity(into, coord, start)`.

## The public surface

| From | Names |
| --- | --- |
| `nimopt` | `COLUMN`, `ROW`, `Absence`, `Alias`, `Assembled`, `Coefficient`, `Constraint`, `Definition`, `Diagnosis`, `Explanation`, `Expression`, `Model`, `Option`, `Param`, `Relation`, `Row`, `Session`, `Set`, `Solution`, `Sum`, `Term`, `Variable`, `available`, `capabilities`, `load`, `loads`, `options`, `product`, `save`, `subset`, `subset_of` |
| `nimblend` | `Array`, `DenseArray`, `Domain`, `EntryBuffer`, `SparseArray`, `combined_dims`, `from_long`, `from_dense`, `is_canonical`, `StoredCoord`, `ProductCoord`, `SubsetCoord` |

**A coefficient composes.** A coefficient is a parameter read at its sets
or an arithmetic combination of such readings: `price[G, T] / eta[G, T]` is
a coefficient written before any data exists, read at its sets like a
parameter, and evaluated once when the matrix is built. `+`, `-`, `*`, `/`
and a power by a number combine coefficients. An expression also carries a
constant, so `x + 1 <= 5` produces the row `x <= 4`.

## What goes wrong

**A chained comparison.** `0 <= expr <= 10` raises `TypeError`. Python
evaluates it as two comparisons joined by `and`, which keeps only the
second, so a relation has no truth value rather than letting the first
bound be dropped. Write each bound as its own constraint.

**A sum over a lag.** `Sum(T - 1, ...)` raises: a sum runs over a set's
members. Put the lag on the variable reference, `x[T - 1]`.

**The built-in `sum` over a set's members.** `sum(x[S, t] for t in members)`
gives the correct answer at a cost: it produces one term per member, where
`Sum(T, x[S, T])` produces one term and reduces a dimension. The terms
concatenate pairwise and each materialises its own block, so building a
model that way runs 24 times slower at 25 members and 275 times at 400, and
the gap widens. Use the built-in `sum` for a short list of distinct
expressions and `Sum` for a set's members.

**A right-hand side over the wrong dimensions.** A constraint's right-hand
side is a parameter over exactly its free dimensions. The error message
gives both.

**Reading values from a model that did not solve.** `objective`, `primal`
and `dual` raise unless `status` is `"optimal"`. Read `status` first.

**A domain over a definition's sets.** `product((B, T))` needs each set's
coordinate, and a declared set has none. In a definition, give `where=`,
`over=` and `subset=` as a tuple of its sets or as one of its parameters,
whose coefficients are the coordinates.

**A row that is not there.** `row()` raises for a coordinate at which the
constraint has no row. `absent()` reports which rule dropped it: a
coefficient absent inside a sum removes a **term** and leaves the row
standing; a term absent along a **free** dimension removes the **row**.

**Reading a MILP's duals.** A model with integer columns has none, and
`dual()` raises rather than returning the relaxation's. Read `primal`.

**A conflict HiGHS cannot prove.** HiGHS computes its conflict over the
linear relaxation, so a model infeasible only through its integrality
produces none and `diagnose()` raises. Gurobi's covers the integrality.

**Two operands that share no dimension.** Every binary operator combines two
dimensioned operands only where they share a dimension, and the rule
applies to a coefficient meeting a variable exactly as it applies to two
coefficients. Frames sharing nothing raise: their combination would be an
outer product no model asks for. A number has no dimension and scales.

**A division by zero.** A divisor that is zero raises `ZeroDivisionError`
with the coordinate, for a Python number, a NumPy scalar and a coefficient
with a zero at one coordinate alike. Handle the divisor before it reaches an
expression.

**A derived coefficient read at the wrong sets.** A combination is read at
its sets as a parameter is, and the reading is checked against the
dimensions it has: `unit_cost[T, G]` raises where it is written, naming
`('G', 'T')`.

**Reaching into `nimblend`.** A test fails on an import from a `nimblend`
submodule, on any read of an array's `.index` or `.data` or a domain's
`.codes`, and on a module of the package assembling an index matrix of its
own.

## Every refusal, and where it is shown

The prose above covers the mistakes worth explaining. This is every
refusal the documentation demonstrates, each executed to produce the
message beside it.

| Raises | Message | Shown at |
| --- | --- | --- |
| `ValueError` | the upper bound 'cap' carries no value for member ('b',) of variable 'x'; a bound covers every column of the variable it bounds | [/guides/bounds-from-parameters](/guides/bounds-from-parameters) |
| `ValueError` | variable 'x' is declared over ('G',) and does not carry ['W']; its upper bound 'cap' is declared over ('W',) | [/guides/bounds-from-parameters](/guides/bounds-from-parameters) |
| `TypeError` | a coefficient is a parameter; build one with `Param.from_dense` or `Param.from_long` and read it at its sets. A product of two expressions is not linear. | [/guides/coefficient-arithmetic](/guides/coefficient-arithmetic) |
| `ValueError` | coefficient (fuel_price / efficiency) is over ('G', 'T'); got ('T', 'G') | [/guides/coefficient-arithmetic](/guides/coefficient-arithmetic) |
| `ZeroDivisionError` | divisor holed carries a zero at 1 coordinate(s), the first at {'G': 'base', 'T': 1}; a quotient there states a coefficient no solver can read | [/guides/coefficient-arithmetic](/guides/coefficient-arithmetic) |
| `ValueError` | frames ('G',) and ('T',) share no dimension; there is nothing to align them on | [/guides/coefficient-arithmetic](/guides/coefficient-arithmetic) |
| `ValueError` | constraint 'capacity' has free dimensions ('P',); its condition is over ('W',) | [/guides/conditions](/guides/conditions) |
| `ValueError` | constraint 'capacity' states its rows with over= and narrows them with where=; state one | [/guides/conditions](/guides/conditions) |
| `ValueError` | variable 'x' is read at member 't9' of dimension 'T', which that set does not carry | [/guides/fixed-members](/guides/fixed-members) |
| `ValueError` | a lag is a whole number of members; got 1.7 | [/guides/lags](/guides/lags) |
| `ValueError` | a sum is over the members of ['T'], so it takes the set and not a lag of it; state the lag at the variable's reference | [/guides/lags](/guides/lags) [/reference/expression](/reference/expression) |
| `ValueError` | parameter 'rate' is read at a lag ['T']; state the lag at the variable's reference, where a coefficient multiplies the row it lands on | [/guides/lags](/guides/lags) |
| `ValueError` | 'max(gen[G, T]) <= 10': Sum is the one call the spelling carries | [/guides/saving-and-loading](/guides/saving-and-loading) |
| `ValueError` | capital does not fall from base to what follows it | [/models/expansion](/models/expansion) |
| `ValueError` | frames ('P',) and ('Q',) share no dimension; there is nothing to align them on | [/nimblend/arrays](/nimblend/arrays) |
| `ValueError` | label column 't' has length 2 and the value column has length 1; they name the same entries | [/nimblend/arrays](/nimblend/arrays) |
| `ValueError` | this array declares absence 'unknown' and does not carry every coordinate of its frame, so densifying must state fill=<value> to place at the rest | [/nimblend/arrays](/nimblend/arrays) [/tutorial/reading-the-answer](/tutorial/reading-the-answer) |
| `ValueError` | 3 member(s) numbered from 4 reach position 6, and dimension 'k' spans 6 | [/nimblend/domains](/nimblend/domains) |
| `ValueError` | a domain of 3 member(s) takes one value each, as a column of that length; got shape (2,) | [/nimblend/domains](/nimblend/domains) |
| `ValueError` | constraint 'supply' has free dimensions ('P',); its right-hand side 'demand' is over ('W',) | [/reference/constraint](/reference/constraint) [/tutorial/constraints](/tutorial/constraints) |
| `ValueError` | data does not cover ['S'] | [/reference/definition](/reference/definition) |
| `ValueError` | parameter 'S' is already declared as a set; a name means one symbol, in an expression and in the data | [/reference/definition](/reference/definition) |
| `TypeError` | a relation has no truth value; a chained comparison such as 0 <= expr <= 10 reads as two comparisons joined by `and` and keeps only the second, so state each bound separately | [/reference/expression](/reference/expression) [/tutorial/constraints](/tutorial/constraints) |
| `TypeError` | a relation is already an equation and states one bound; compare the expression a second time in its own equation rather than comparing the relation | [/reference/expression](/reference/expression) |
| `TypeError` | an LP has no row for a strict inequality; state `<=` or `>=`. `min` and `max` compare two expressions this way and are not linear either, so reduce with `Sum` over the sets instead | [/reference/expression](/reference/expression) |
| `TypeError` | an expression is reduced over the sets it is summed across; state them with `Sum(I, J, expression)` | [/reference/expression](/reference/expression) |
| `TypeError` | nimopt expresses a linear term, so a variable in a denominator is not one; state the reciprocal as a coefficient the variable multiplies | [/reference/expression](/reference/expression) |
| `TypeError` | nimopt expresses a linear term, so a variable raised to a power is not one; a coefficient takes the power instead, and a variable multiplies it | [/reference/expression](/reference/expression) |
| `TypeError` | nimopt expresses a linear term, so the absolute value of one is not linear; reduce with `Sum` over its sets, or state the magnitude with two rows bounding the expression | [/reference/expression](/reference/expression) |
| `ValueError` | term 'x' already sums over ['T']; a dimension is reduced once, and a second reduction has nothing left to reduce | [/reference/expression](/reference/expression) |
| `ValueError` | constraint 'cap' states where= with a domain that has no name; declare its members as a parameter and name that | [/reference/files](/reference/files) |
| `ValueError` | parameter 'c' is given columns ['value', 'S']; a table states the dimensions then value: ['S', 'value'] | [/reference/files](/reference/files) |
| `ValueError` | variable 'x' carries ['bound'], which the format does not; it takes ('sets', 'subset', 'lower', 'upper', 'integer') | [/reference/files](/reference/files) |
| `ValueError` | constraint 'cap' states no row at {'P': 'p3'}; `absent('cap')` names the rule that dropped it | [/reference/inspection](/reference/inspection) |
| `ValueError` | parameter 'cost': label column 'P' has length 1 and the value column has length 2; they name the same entries | [/reference/param](/reference/param) |
| `TypeError` | parameter 'price' carries ('G',) and states no coefficient until it is read; read it at its sets as price[G] | [/reference/param](/reference/param) |
| `ValueError` | the model's status is 'infeasible', so it carries no objective; read `status` before reading values | [/reference/solution](/reference/solution) [/tutorial/solving](/tutorial/solving) |
| `ValueError` | this model carries integer columns and 'highs' refuses duals for a model with integrality, so there is no dual for constraint 'cap' to read: a mixed-integer model's duals are not its relaxation's | [/reference/solvers](/reference/solvers) |
| `TypeError` | parameter 'supply' carries ('P',) and states no coefficient until it is read; read it at its sets as supply[P] | [/tutorial/constraints](/tutorial/constraints) |
| `ValueError` | parameter 'cost' is over sets of shape (2, 3); got values of shape (2, 2) | [/tutorial/sets-and-parameters](/tutorial/sets-and-parameters) |

---

# /

# nimopt

`nimopt` is a Python library for building linear and mixed-integer programs. A model is declared symbolically over named index sets — as parameters, variables and constraints — and is expanded into a coefficient matrix only when it is assembled or solved. Solutions are returned as arrays over those same index sets, so a primal value is read by label rather than by column position.

`nimopt` is built on `nimblend`, a labelled sparse N-dimensional array library with no knowledge of optimisation. The dependency runs in one direction, and `nimblend` is documented in [its own section](/nimblend).

## Design

**A variable is a dimension.** A constraint is an array indexed over its free sets crossed with the model's column space, with the coefficients as values. There is no assembly step converting the model into a matrix, because the array is the matrix.

**Absence is distinct from zero.** An entry is either stored or absent, and every array declares what absence means: `"empty"` for a coordinate that contributes nothing, `"unknown"` for one that was never modelled. A missing result is never counted as zero, and division by an absent value raises rather than producing an infinity.

**A subset stays a subset.** A variable declared over a subset of a set product has one column per member of the subset and none for the rest. The full product is never materialised, at declaration or at any point after it.

**Expressions are symbolic.** An expression holds references to variables and parameters rather than their values. `cost[P, W] * x[P, W]` costs the same to write over a million routes as over six; the values are read when the matrix is built.

**Dropped rows are reported.** A row whose terms have no value at some coordinate is dropped rather than written incompletely. `absent()` lists every dropped row with the rule that dropped it, and `row()` returns one row of the assembled matrix as the solver receives it.

## Install

Neither package is published yet, so both are installed from a checkout. `nimblend` is a dependency and is installed first, from wherever it is cloned; `nimopt` then installs from its own root:

```bash
pip install /path/to/nimblend
pip install ".[highs]"
```

HiGHS is the default solver, and `[highs]` installs it. `[gurobi]` and `[mosek]` add those adapters instead, `[bench]` adds the comparison suite and `[dev]` the test and lint tooling. `available()` reports the solvers whose backend can be imported in the current environment, and `capabilities(name)` answers for an adapter whether or not its backend is installed.

## A first model

A transport problem: two plants with limited supply ship to three warehouses with fixed demand, and the objective is total shipping cost.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))

m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

solution = m.solve()
print(solution.status, solution.objective)
print(solution.primal("x").to_dense())
```

Output:

```text
optimal 135.0
[[20.  0. 10.]
 [ 0. 15.  5.]]
```

The primal values are returned as a 2 by 3 array over plants and warehouses, in the order the sets declare their members.

## Declaring before the data exists

A `Definition` states the same model without binding any data. Its sets and parameters are declared by name, its constraints are written in the same expression syntax, and `explain()` reports the whole declaration before a single value has been read.

```python
import numpy as np
from nimopt import Definition, Sum

d = Definition("transport")
P, W = d.set("P"), d.set("W")
cost = d.param("cost", (P, W))
supply = d.param("supply", (P,))
demand = d.param("demand", (W,))
x = d.var("x", (P, W))

d.eq("supply", Sum(W, x[P, W]) <= supply[P])
d.eq("demand", Sum(P, x[P, W]) >= demand[W])
d.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

print(d.explain())
```

Output:

```text
transport  min  not built
  sets        P · W
  parameters  cost (P,W) · supply (P) · demand (W)
  variables   x (P×W) [0.0, inf]
  constraint  supply (P)  Sum(W, x[P, W]) <= supply[P]
  constraint  demand (W)  Sum(P, x[P, W]) >= demand[W]
  objective   min  Sum(P, W, cost[P, W] * x[P, W])
```

A definition is copied before it is bound, so one definition builds as many models as it is given datasets for and is unchanged by any of them. The built model is inspected the same way, and `row()` reads one row back out of the assembled matrix as the solver receives it.

```python
import numpy as np
from nimopt import Definition, Sum

d = Definition("transport")
P, W = d.set("P"), d.set("W")
cost = d.param("cost", (P, W))
supply = d.param("supply", (P,))
demand = d.param("demand", (W,))
x = d.var("x", (P, W))
d.eq("supply", Sum(W, x[P, W]) <= supply[P])
d.eq("demand", Sum(P, x[P, W]) >= demand[W])
d.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

data = {
    "P": np.array(["lisbon", "porto"]),
    "W": np.array(["berlin", "paris", "rome"]),
    "cost": np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]),
    "supply": np.array([30.0, 25.0]),
    "demand": np.array([20.0, 15.0, 15.0]),
}
m = d.build(data)
print(m)
print(m.row("demand", W="paris"))
print(m.absent("demand"))
```

Output:

```text
Model('transport', 1 variables, 6 columns, 5 rows)
demand[W='paris']  row 3
  1·x[lisbon,paris] + 1·x[porto,paris] >= 15
demand  3 of 3 rows  stated by terms
```

## A variable over a subset

Where a variable spans an arc list rather than a full product, it has a column per arc and the product is never built. A thousand plants each serving three warehouses is three thousand columns, not a million.

```python
import numpy as np
from nimopt import Model, Set, subset

P = Set("P", np.array([f"p{i}" for i in range(1000)]))
W = Set("W", np.array([f"w{i}" for i in range(1000)]))

served = np.array([f"w{(i * 7 + k) % 1000}" for i in range(1000) for k in range(3)])
arcs = subset((P, W), {"P": np.repeat(P.labels, 3), "W": served})

m = Model("transport")
x = m.var("x", (P, W), subset=arcs)
print(f"{m.n_columns} columns over a product of {len(P) * len(W)}")
```

Output:

```text
3000 columns over a product of 1000000
```

## Features

- Sets, aliases, subsets and set products as the index structure of every declaration
- Parameters from dense arrays or long-form columns, broadcast where a parameter is narrower than the variable it multiplies
- Composable coefficients: a parameter read at its sets, or an arithmetic of parameters written before data exists
- Conditions on a sum and on a constraint, lags that drop or wrap at the ends of a set, and members fixed at a label
- Per-column bounds from a parameter, and variables declared over a subset of a set product
- A `Definition` written before data exists and built against any number of datasets
- `explain()` on a definition or a built model, `row()` into the assembled matrix, and `absent()` reporting dropped rows and the rule that dropped each
- A `Session` that keeps the solved instance open, and `diagnose()` reporting the conflicting rows of an infeasible model or the ray of an unbounded one
- Primals and duals returned over their index sets, with absence distinct from zero
- Continuous and integer columns, solved through HiGHS, Gurobi or Mosek behind one adapter contract, with `capabilities()` stating what each adapter does and which capabilities it refuses together
- One option vocabulary translated into each solver's own spelling, so a time limit is stated the same way whatever solves the model
- A model written to and read back from YAML, with its data inline or in a sidecar
- A corpus of worked models under `nimopt.models`, each stating its formulation, inputs at any size, and an objective computed by arithmetic rather than by a solver

## Performance

Where a variable's columns are a subset of a set product, not materialising the product is worth a great deal. On a transport model of 400 000 arcs over a 20 000 000-cell product, `nimopt` builds the same matrix in 72.9 MB of resident memory against a dense rival's 1 682.9 MB, and in 292.9 ms against 844.5 ms.

Where nothing is sparse, the alignment work is a cost with no corresponding saving. On a fully dense temporally-coupled model at 2 111 080 rows, the same comparison reverses: the rival builds the matrix three times faster, for seven per cent more resident memory.

Both rows are in the suite for the same reason. The [benchmark page](/explanation/what-the-numbers-measure) gives each figure, the baseline it is measured against, and what it does not claim.

## Documentation

- [Get started](/get-started): installation, and the transport model above solved and read back.
- [Vocabulary](/vocabulary): the terms used throughout the documentation.
- [Tutorial](/tutorial/sets-and-parameters): the transport model built in six steps, one concept per page.
- [Playground](/playground): every example runs in the browser and can be edited.
- [For agents](/for-agents): the mental model, the public surface and the failure modes on one page.

Every Python example on this site is executed by the test suite and shows the output it produced.

## Licence

MIT. See `LICENSE`.

## Citing

The package carries a `CITATION.cff`. Cite it by author, name and version:

> Gaete-Morales, Carlos. *nimopt* (version 0.1.2). MIT.

## Contributing

Issues and patches are welcome once the repositories are published. Until then, the most useful contribution is a model that does not fit: the formulations that are awkward to write determine the next features.

---

# /get-started

# Get started

## Install

```bash
pip install nimopt
```

`nimblend` is installed as a dependency. HiGHS is the default solver.

## A transport model

Two plants, Lisbon and Porto, ship to three warehouses, Berlin, Paris and
Rome. Plant `p` has supply `s[p]`, warehouse `w` has demand `d[w]`, and one
unit shipped on route `(p, w)` costs `c[p, w]`. The decision variable
`x[p, w]` is the quantity shipped on each route.

```text
minimise    Σ_{p,w} c[p,w] · x[p,w]
subject to  Σ_w x[p,w] ≤ s[p]        for each plant p
            Σ_p x[p,w] ≥ d[w]        for each warehouse w
            x[p,w] ≥ 0
```

In `nimopt`, the sets index every declaration, the parameters hold the data,
`m.var` declares the decision variable, `m.eq` adds each constraint family
under a name, and `set_objective` sets the objective function.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))

m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

solution = m.solve()
print(solution.status, solution.objective)
```

Output:

```text
optimal 135.0
```

The solver reports an optimal solution with objective 135.

## Reading the solution

`primal("x")` returns the shipments as an array indexed over the sets `x`
was declared on. `dual("demand")` returns the dual value of each demand row:
the change in the objective per unit increase in that warehouse's demand.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

shipped = solution.primal("x")
print(shipped.dims)
print(shipped.to_dense())
print(solution.dual("demand").to_dense())
```

Output:

```text
('P', 'W')
[[20.  0. 10.]
 [ 0. 15.  5.]]
[3. 1. 6.]
```

Rows are plants and columns are warehouses. Lisbon ships 20 to Berlin and 10
to Rome; Porto ships 15 to Paris and 5 to Rome. The duals of the demand rows
are 3, 1 and 6: the marginal cost of one additional unit at each warehouse.

Every code block in this documentation is self-contained, which is why the
second block repeats the model. Each block can be pasted into a Python
session as it is, or opened in the playground with "Run this example".

## Next

- [Vocabulary](/vocabulary) defines the terms used throughout: set, member,
  frame, row, absence, and others.
- The [tutorial](/tutorial/sets-and-parameters) builds this model one concept
  per page.
- The [guides](/guides/subsets) cover sparse networks, time lags, conditions
  on rows, and bounds from data.
- [Explanation](/explanation/a-variable-is-a-dimension) covers the design and
  its costs.

---

# /vocabulary

# Vocabulary

Terms used throughout the documentation, each defined once. Examples refer
to the transport model: plants `P = {lisbon, porto}` ship to warehouses
`W = {berlin, paris, rome}`.

## Index sets

**Set.** A named index dimension with labels. `Set("P", np.array(["lisbon",
"porto"]))` is the set of plants. Parameters, variables and constraints are
indexed over sets, and solution values are returned over the same sets.

**Member.** One element of a set. `"lisbon"` is a member of `P`.

**Label.** The name of a member, a string or a number. Labels are the
caller-facing identifiers; integer positions are used internally.

**Set product.** The Cartesian product of several sets. `P × W` has six
members, `("lisbon", "berlin")`, `("lisbon", "paris")` and so on. Variables
and parameters are indexed over set products.

**Subset.** An explicit list of members of a set product. `subset((P, W),
{"P": ..., "W": ...})` lists the routes that exist. A variable over a subset
has a column per listed member and none for the rest.

**Alias.** A second name for a set, sharing its labels. It allows a
parameter or a constraint to relate a set to itself, such as a flow between
two nodes of one node set.

**Domain.** A set of coordinates over some dimensions. `product((P, W))` and
`subset(...)` return one. `subset=`, `where=` and `over=` take a domain.

## Data and decisions

**Parameter.** Data indexed over a set product: one value per member.
`cost` is indexed over `(P, W)`; `supply` over `P`. A parameter has no
column in the matrix.

**Coefficient.** The multiplier of a variable in a row. A parameter indexed
at its sets, `cost[P, W]`, is a coefficient, as is an arithmetic
combination of such readings, `price[G, T] / eta[G, T]`.

**Variable.** A decision variable. `m.var("x", (P, W))` declares one
decision per route. Values are read after a solve with `primal("x")`.

**Column.** One decision in the coefficient matrix. Each member of a
variable is one column. Column indices are computed from member positions
and are never assigned by the caller.

**Bound.** The interval a column may take values in. The default lower
bound is 0 and the default upper bound is infinity.

## Expressions

**Expression.** A linear combination of variables, such as `Sum(W, x[P,
W])`. Writing an expression records its structure and computes nothing.
Values are read when the matrix is built.

**Term.** One component of an expression: one variable, an optional
coefficient, and the sets summed over. `cost[P, W] * x[P, W]` is one term.

**Frame.** The dimensions an expression is still indexed over, also called
its free dimensions. `x[P, W]` has frame `(P, W)`; `Sum(W, x[P, W])` has
frame `(P,)`. An empty frame is a scalar.

**Sum.** Summation over the members of the named sets. The summed sets are
removed from the frame.

**Lag.** A reference to the previous or next member of a set. `x[T - 1]`
references the previous period. A lag either drops the row with no
predecessor or, with `T.cyclic`, wraps to the last member.

**Fixed member.** A label in place of a set in a reference, `x[G, "t0"]`.
It selects that member and removes the set from the frame.

## Constraints and the matrix

**Relation.** An expression compared with `<=`, `>=` or `==` to a
right-hand side. `Sum(W, x[P, W]) <= supply[P]` is a relation. It becomes
of the model when passed to `m.eq`.

**Constraint.** A relation added to the model under a name. It produces one
row per member of its expression's frame.

**Row.** One inequality or equality of the coefficient matrix. The supply
constraint over two plants produces two rows.

**Right-hand side.** The scalar or parameter on the other side of the
relation. A scalar applies to every row. A parameter must be indexed over
exactly the constraint's frame, so that each row has its own value.

**Objective.** A scalar expression, one with an empty frame, that the solver
minimises or maximises. `Sum(P, W, cost[P, W] * x[P, W])` is the total
shipping cost.

**Sense.** The optimisation direction, `"min"` or `"max"`, set once on the
`Model`.

**Materialise.** Evaluate a parameter or an expression into an array of
values. This happens when the matrix is built, not when the expression is
written.

**Assemble.** Build the coefficient matrix from every constraint. `solve()`
assembles before calling the solver; `assemble()` returns the matrix without
solving.

**Nonzero.** One stored coefficient of the matrix. `nnz` is the count.

## Solutions

**Solution.** The return value of `solve()`: a status, an objective value,
and primal and dual values.

**Status.** The outcome the solver reported: `optimal`, `infeasible`,
`unbounded`, or a limit reached. Values are defined only for `optimal`.

**Primal.** The value of a variable in the solution, returned over the sets
it was declared on.

**Dual.** The dual value of a constraint, also called the shadow price: the
change in the objective per unit change in that row's right-hand side.
Returned over the constraint's frame.

**Absence.** A coordinate at which an array has no value, as distinct from
a stored zero. Every array declares the meaning of absence: `"empty"` for a
coordinate that contributes nothing, used by parameters, or `"unknown"` for
one that was never modelled, used by solutions.

**Session.** A solver instance kept open on one assembled model, so that
questions can be asked of it after the solve, such as which rows conflict.

**Definition.** A model written before its data exists, in the same
vocabulary. `build(data)` produces a `Model` for one dataset.

---

# /reference/constraint

# Constraint

## `Constraint`

Returned by `Model.eq`. Rows over an expression's frame, bounded by a
right-hand side.

```
Model.eq(name, relation, where=None, over=None)
```

| Argument | Meaning |
| --- | --- |
| `name` | the name `Solution.dual` reads it back by |
| `relation` | an expression, a sense and a right-hand side |
| `where` | a domain intersecting the rows |
| `over` | the rows, given explicitly |

| Member | Returns |
| --- | --- |
| `n_rows` | the number of rows it produces |
| `nnz` | the number of coefficients they hold |
| `row_of(name)` on the `Assembled` | where those rows sit in the matrix |

A row derived from the terms exists where every term has a value and the
right-hand side has a value. A coefficient absent inside a sum removes a
term and leaves the row standing; a term absent along a free dimension
removes the row, because a row missing one of its terms would express a
constraint that was not written.

`over=` gives the rows explicitly instead, so a term covering some of them
contributes where it has values. A condition given with `where=`
intersects the row domain, so a row outside the condition is not produced.

The expression is symbolic, so the constraint holds the recipe rather than
a block: it is materialised once to compute its shape and once to write it,
and holds nothing in between.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))

m = Model("transport")
x = m.var("x", (P, W))

rows = m.eq("supply", Sum(W, x[P, W]) <= supply[P])
print(rows.n_rows, rows.nnz)
print(m.assemble().row_of("supply"))
```

Output:

```text
2 6
slice(0, 2, None)
```

The right-hand side is a number, applied to every row, or a parameter over
exactly the constraint's free dimensions, giving each row its own value. A
parameter over other dimensions raises `ValueError`; the message gives the
constraint's free dimensions and the parameter's.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))

m.eq("supply", Sum(W, x[P, W]) <= demand[W])
```

Raises ValueError:

```text
ValueError: constraint 'supply' has free dimensions ('P',); its right-hand side 'demand' is over ('W',)
```

---

# /reference/definition

# Definition

## `Definition`

```
Definition(name="definition", sense="min")
```

A definition declares the sets, parameters and variables a model is written
from, and its constraints, in the expression syntax a model uses. It holds
no data: a set declared here names a dimension and has no members, and a
parameter names a shape and has no values.

An expression holds references rather than arrays, so an equation's free
dimensions and its sense are read off the relation rather than declared
beside it. `sense` is `"min"` or `"max"`, set once here.

| Member | Returns |
| --- | --- |
| `set(name)` | a declared `Set`, whose members arrive with the data |
| `alias(name, base)` | a declared `Alias` over one of this definition's sets |
| `param(name, sets)` | a declared `Param`, whose values arrive with the data |
| `var(name, sets, subset=None, lower=0.0, upper=inf, integer=False)` | a declared `Variable` |
| `eq(name, relation, where=None, over=None)` | nothing; registers the equation |
| `build(data)` | a `Model` over the declarations, bound to `data` |
| `explain()` | an `Explanation` of what is declared |
| `to_yaml()` | the text of this definition's file, structure and no data |
| `set_objective(expression)` | nothing; sets the objective |
| `sense` | `"min"` or `"max"`, as declared |
| `sets`, `aliases`, `parameters`, `variables`, `constraints` | the registries, keyed by name |

```python
from nimopt import Definition, Sum

d = Definition("dispatch", sense="min")
snapshot = d.set("snapshot")
generator = d.set("generator")
p_max = d.param("p_max", (generator,))
load = d.param("load", (snapshot,))
cost = d.param("cost", (generator,))
p = d.var("p", (snapshot, generator), lower=0.0, upper=p_max)
d.eq("balance", Sum(generator, p[snapshot, generator]) == load[snapshot])
d.set_objective(Sum(snapshot, generator, cost[generator] * p[snapshot, generator]))

print(list(d.sets), list(d.parameters))
print(d.constraints["balance"][0].expression.frame)
print(list(d.variables), d)
```

Output:

```text
['snapshot', 'generator'] ['p_max', 'load', 'cost']
('snapshot',)
['p'] Definition('dispatch', 1 variables, 1 constraints)
```

## One namespace for sets and parameters

Sets and parameters share one key space, because the data a definition is
built from is keyed by declared name and one would otherwise shadow the
other. Declaring a parameter under a set's name raises `ValueError`.

```python raises=ValueError
from nimopt import Definition

d = Definition("d")
S = d.set("S")
d.param("S", (S,))
```

Raises ValueError:

```text
ValueError: parameter 'S' is already declared as a set; a name means one symbol, in an expression and in the data
```

Equations are in no data mapping, so a constraint may take the name of the
parameter that bounds it.

```python
from nimopt import Definition, Sum

d = Definition("d")
S = d.set("S")
supply = d.param("supply", (S,))
one = d.param("one", (S,))
x = d.var("x", (S,))
d.eq("supply", Sum(S, one[S] * x[S]) <= supply[S])

print(list(d.parameters), list(d.constraints))
```

Output:

```text
['supply', 'one'] ['supply']
```

## An alias in a definition

`alias(name, base)` declares a second name for one of the definition's sets,
which is how a model relates a set to itself. The alias carries no data of
its own: it reads the labels its base set binds, so `build` takes members for
the set and none for the alias, and naming the alias in `data` is refused.

```python
import numpy as np
from nimopt import Definition, Sum

d = Definition("network", sense="min")
N = d.set("N")
NP = d.alias("NP", N)
limit = d.param("limit", (N, NP))
flow = d.var("flow", (N, NP), lower=0.0)
d.eq("cap", flow[N, NP] <= limit[N, NP])
d.set_objective(Sum(N, NP, limit[N, NP] * flow[N, NP]))

m = d.build({"N": np.array(["a", "b"]), "limit": np.ones((2, 2))})
print(m.n_columns, m.n_rows)
```

Output:

```text
4 4
```

## Domains in a definition

A `Domain` resolves labels through each set's coordinate, and a declared
set has none. `where=` and `over=` on `eq`, and `subset=` on `var`,
therefore take a tuple of the definition's sets, meaning their full
product, or one of its parameters, whose coefficients are the coordinates.
Both forms resolve to the same domain, so a model and a definition declare
a sparse variable or an explicit row domain the same way.

## Building

`build(data)` copies the declaration graph, binds the copy, numbers the
columns and returns a `Model`. `data` maps a declared set's name to its
members and a declared parameter's name to its values. The definition is
unchanged, so it builds as many models as it is given datasets.

A parameter's values arrive dense over its product, as an array of one
value per cell, or long over its entries, as a pair of one mapping of label
columns and one value column. The long form is how a parameter with
coefficients at some coordinates and none at the rest is given, and it is
what a variable declared with `subset=` that parameter takes its members
from.

```python
import numpy as np
from nimopt import Definition, Sum

d = Definition("transport", sense="min")
P, W = d.set("P"), d.set("W")
cost = d.param("cost", (P, W))
supply = d.param("supply", (P,))
demand = d.param("demand", (W,))
flow = d.var("flow", (P, W), subset=cost, lower=0.0)
d.eq("supply", Sum(W, cost[P, W] * flow[P, W]) <= supply[P])
d.eq("demand", Sum(P, cost[P, W] * flow[P, W]) >= demand[W])
d.set_objective(Sum(P, W, cost[P, W] * flow[P, W]))

m = d.build(
    {
        "P": np.array(["p1", "p2"]),
        "W": np.array(["w1", "w2"]),
        "cost": (
            {"P": np.array(["p1", "p1", "p2"]), "W": np.array(["w1", "w2", "w1"])},
            np.array([1.0, 2.0, 3.0]),
        ),
        "supply": np.array([3.0, 3.0]),
        "demand": np.array([1.0, 1.0]),
    }
)

# three arcs, so three columns rather than the four the product would span
print(m.n_columns, m.n_rows)
print(m.solve().status)
```

Output:

```text
3 4
optimal
```

Data that misses a declaration, or names something the definition never
declared, raises `ValueError` before anything is bound.

```python raises=ValueError
from nimopt import Definition

d = Definition("d")
d.set("S")
d.build({})
```

Raises ValueError:

```text
ValueError: data does not cover ['S']
```

---

# /reference/explanation

# Explanation

## `Explanation`

Returned by `Definition.explain` and `Model.explain`. A frozen record of
every declaration and what it built. The shapes it is made of are frozen
too, so a reader takes a field rather than parsing a rendering.

| Field | Holds |
| --- | --- |
| `name`, `sense` | the declaration's name and the direction it optimises |
| `built` | whether counts are facts about data or absent |
| `sets` | one `SetShape` per dimension |
| `parameters` | one `ParamShape` per parameter |
| `variables` | one `VariableShape` per variable |
| `constraints` | one `ConstraintShape` per equation |
| `objective` | the objective's spelling, or `None` |
| `columns`, `rows`, `nonzeros` | the model's shape, or `None` |

A count is `None` where nothing is bound. It is never zero: a count of zero
is a fact a caller acts on, and reporting one for a declaration would be
false.

| Shape | Fields |
| --- | --- |
| `SetShape` | `name`, `size` |
| `ParamShape` | `name`, `dims`, `entries` |
| `VariableShape` | `name`, `dims`, `members`, `columns`, `lower`, `upper`, `integer` |
| `ConstraintShape` | `name`, `free`, `sense`, `rows`, `nonzeros`, `relation` |

`VariableShape.members` names the parameter a sparse variable took its
members from, and is `None` for one over the full product. Columns are
absent until data binds, so without it a sparse declaration and a dense one
would otherwise read identically.

`ConstraintShape.free` and `.sense` are read off the relation rather than
declared beside it, because an expression holds references and reports
both.

```python
from nimopt import Definition, Sum

d = Definition("transport", sense="min")
P, W = d.set("P"), d.set("W")
cost = d.param("cost", (P, W))
supply = d.param("supply", (P,))
flow = d.var("flow", (P, W), subset=cost, lower=0.0)
d.eq("supply", Sum(W, cost[P, W] * flow[P, W]) <= supply[P])
d.set_objective(Sum(P, W, cost[P, W] * flow[P, W]))

e = d.explain()
print(e.built, e.columns, e.variables[0].members)
print(e.constraints[0].free, e.constraints[0].sense)
print(e)
```

Output:

```text
False None cost
('P',) <=
transport  min  not built
  sets        P · W
  parameters  cost (P,W) · supply (P)
  variables   flow (P×W) over cost [0.0, inf]
  constraint  supply (P)  Sum(W, cost[P, W] * flow[P, W]) <= supply[P]
  objective   min  Sum(P, W, cost[P, W] * flow[P, W])
```

---

# /reference/expression

# Expressions

## `Term`

One variable, an optional coefficient, the dimensions summed over, and a
scale factor. A term is the recipe for a block of coefficients and holds
references rather than arrays, so writing it costs nothing: an expression
over a million columns costs the same as one over ten.

| Member | Returns |
| --- | --- |
| `free_dims` | the dimensions it is still indexed over |
| `carried_dims` | every dimension it has |
| `with_coefficient(coefficient)` | the term, scaled by a parameter |
| `summing(dims)` | the term, reduced over those dimensions |
| `scaled(by)` | the term, multiplied by a number |
| `restricted_to(domain)` | the term, over those members only |

A caller builds terms through the operators rather than these members:
`cost[P, W] * x[P, W]` gives a coefficient, `Sum` gives the reduction, and
`-` gives the scale.

## `Expression`

A list of terms and the frame they share. The frame is the union of the
terms' free dimensions, ordered by the term that introduces each. A term
narrower than the frame is broadcast over it when the expression is
materialised.

| Member | Returns |
| --- | --- |
| `terms` | the terms it holds |
| `frame` | the dimensions it is indexed over |
| `coords` | the coordinates of that frame |
| `materialise()` | its coefficients as a `nimblend` array |

```python
import numpy as np
from nimopt import Model, Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))

m = Model("transport")
x = m.var("x", (P, W))
y = m.var("y", (P, W))

combined = cost[P, W] * x[P, W] - y[P, W]
print(combined.frame)
print(len(combined.terms))
print(combined.terms[0].free_dims)
```

Output:

```text
('P', 'W')
2
('P', 'W')
```

## `Sum`

```
Sum(I, J, ..., expression, where=None)
```

The expression reduced over the named sets. Each set named leaves the
frame. `where=` takes a domain and restricts each term's entries before the
reduction, so the sum runs over the coordinates given rather than every
coordinate of the product.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

print(x[P, W].frame)
print(Sum(W, x[P, W]).frame)
print(Sum(P, W, x[P, W]).frame)
```

Output:

```text
('P', 'W')
('P',)
()
```

A sum runs over the members of a set, so it takes the set and not a lag of
it. The lag belongs on the variable reference.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum

T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (T,))

Sum(T - 1, x[T])
```

Raises ValueError:

```text
ValueError: a sum is over the members of ['T'], so it takes the set and not a lag of it; state the lag at the variable's reference
```

## `Relation`

An expression, a sense and a right-hand side, produced by comparing an
expression with `<=`, `>=` or `==`. `Model.eq` turns one into a
constraint.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

bounded = Sum(W, x[P, W]) <= 30.0
print(type(bounded).__name__, bounded.sense)
```

Output:

```text
Relation <=
```

A relation has no truth value. Python evaluates `0 <= expr <= 10` as two
comparisons joined by `and`, which keeps only the second, so the chained
form raises rather than letting the first bound be dropped.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

0.0 <= Sum(W, x[P, W]) <= 10.0
```

Raises TypeError:

```text
TypeError: a relation has no truth value; a chained comparison such as 0 <= expr <= 10 reads as two comparisons joined by `and` and keeps only the second, so state each bound separately
```

## Forms that are not linear

`nimopt` expresses linear terms. Each form below raises where it is
written, and the message gives the form to write instead.

A variable raised to a power is not linear; a coefficient takes the power
and a variable multiplies it.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
x[T] ** 2
```

Raises TypeError:

```text
TypeError: nimopt expresses a linear term, so a variable raised to a power is not one; a coefficient takes the power instead, and a variable multiplies it
```

A variable in a denominator is not linear either; the reciprocal is written
as a coefficient the variable multiplies.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
1.0 / x[T]
```

Raises TypeError:

```text
TypeError: nimopt expresses a linear term, so a variable in a denominator is not one; state the reciprocal as a coefficient the variable multiplies
```

The absolute value of an expression is not linear. A magnitude is written
with two rows bounding the expression, and a reduction with `Sum` over its
sets.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
abs(x[T])
```

Raises TypeError:

```text
TypeError: nimopt expresses a linear term, so the absolute value of one is not linear; reduce with `Sum` over its sets, or state the magnitude with two rows bounding the expression
```

An LP has no row for a strict inequality.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
x[T] < 5.0
```

Raises TypeError:

```text
TypeError: an LP has no row for a strict inequality; state `<=` or `>=`. `min` and `max` compare two expressions this way and are not linear either, so reduce with `Sum` over the sets instead
```

The built-in `sum` of expressions with no set to reduce over reaches
`Expression.sum`, which would otherwise return the expression unchanged
having reduced nothing.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
x[T].sum()
```

Raises TypeError:

```text
TypeError: an expression is reduced over the sets it is summed across; state them with `Sum(I, J, expression)`
```

A relation expresses one bound. Comparing it a second time raises rather
than dropping the first.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
(x[T] <= 5.0) >= 1.0
```

Raises TypeError:

```text
TypeError: a relation is already an equation and states one bound; compare the expression a second time in its own equation rather than comparing the relation
```

A dimension is reduced once; a second reduction has nothing left to
reduce.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum

T = Set("T", np.arange(3))
m = Model("m")
x = m.var("x", (T,))
Sum(T, Sum(T, x[T]))
```

Raises ValueError:

```text
ValueError: term 'x' already sums over ['T']; a dimension is reduced once, and a second reduction has nothing left to reduce
```

---

# /reference/files

# Files

## `load`, `loads`, `save`

| Verb | Does |
| --- | --- |
| `load(path, data=None)` | reads a file; a `Definition`, or a `Model` where the file carries data or `data=` gives it |
| `loads(text, data=None)` | the same over text; a sidecar name in text is refused, because text has no directory |
| `save(what, path, inline=False)` | writes a definition's file, or a model's with an `.npz` beside it, or one file with an inline block when `inline=True` |

`data=` is the mapping `build` takes or the path of an `.npz`. A file that
carries data and a `data=` together is refused.

`Definition.to_yaml()` and `Model.to_yaml(inline=False)` return the text
`save` writes, without a sidecar line: only `save` writes a sidecar and the
line that names it.

## The file

| Key | Carries |
| --- | --- |
| `version` | `2`; any other value is refused naming the one this reader understands |
| `name`, `sense` | the model's |
| `sets` | a list of names |
| `aliases` | each alias to the set it names; absent where the model declares none |
| `parameters` | each name to its dimensions |
| `variables` | each name to `sets`, and to `subset`, `lower`, `upper`, `integer` where they differ from no subset, `0`, infinity and `false` |
| `constraints` | each name to `relation`, and to `where` or `over` where given |
| `objective` | the expression's spelling; absent where the model states none |
| `data` | an inline mapping, or the name of an `.npz` beside the file |

`subset`, `where` and `over` take a parameter's name, meaning the coordinates
it carries, or a list of set names, meaning their full product. A symbol's
name is a Python identifier other than `Sum`. A dimension named in
`parameters`, in a variable's `sets`, or in a list of set names is a set or
an alias; an alias is declared after the set it names.

The expressions are spelled as they are typed in Python and read back through
the same operators, so the file is a fixed point: reading it and writing it
again gives the same text.

```python
from nimopt import Definition, Sum, loads

d = Definition("d")
S = d.set("S")
c = d.param("c", (S,))
x = d.var("x", (S,), integer=True)
d.eq("cap", 2 * c[S] * x[S] - 1 <= 5)
text = d.to_yaml()
print(text)
print(loads(text).to_yaml() == text)
```

Output:

```text
version: 2
name: d
sense: min
sets: [S]
parameters:
  c: [S]
variables:
  x:
    sets: [S]
    integer: true
constraints:
  cap:
    relation: (c[S] * 2) * x[S] - 1 <= 5

True
```

A key the format does not carry is refused, at the top level and inside an
entry.

```python raises=ValueError
from nimopt import loads

loads(
    "version: 2\nname: d\nsense: min\nsets: [S]\n"
    "variables:\n  x: {sets: [S], bound: 1}\n"
)
```

Raises ValueError:

```text
ValueError: variable 'x' carries ['bound'], which the format does not; it takes ('sets', 'subset', 'lower', 'upper', 'integer')
```

## Data

The data a file carries is the mapping `build` takes, in three shapes.

| Shape | Inline | In the `.npz` |
| --- | --- | --- |
| a set's members | a list | a one-dimensional label array |
| a dense parameter | nested lists in row-major order | its grid |
| a long parameter | `columns`, the dimensions then `value`, and `rows` | a structured array with one field per dimension and `value` |

A parameter is written dense where its array covers its full product and
long otherwise. The `.npz` is read with `allow_pickle=False`; an array of
object dtype is refused at save, naming the symbol, because the container
would pickle it silently and the reader would then refuse the file.

```python raises=ValueError
from nimopt import loads

loads(
    "version: 2\nname: d\nsense: min\nsets: [S]\nparameters:\n  c: [S]\n"
    "data:\n  S: [a]\n  c:\n    columns: [value, S]\n    rows:\n    - [1.0, a]\n"
)
```

Raises ValueError:

```text
ValueError: parameter 'c' is given columns ['value', 'S']; a table states the dimensions then value: ['S', 'value']
```

## What is refused before anything is written

| Written | Refused because |
| --- | --- |
| a model whose `subset`, `where` or `over` is a domain with no name | declare the members as a parameter |
| two parameter objects or two set objects sharing a name in one model | the file keys a symbol by name |
| a symbol whose name is not an identifier, or is `Sum` | the spelling cannot address it |
| an array of object dtype | the container would pickle it |

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, subset

P = Set("P", np.array(["a", "b"]))
m = Model("m")
x = m.var("x", (P,))
m.eq("cap", x[P] <= 1.0, where=subset((P,), {"P": np.array(["a"])}))
m.to_yaml()
```

Raises ValueError:

```text
ValueError: constraint 'cap' states where= with a domain that has no name; declare its members as a parameter and name that
```

---

# /reference/inspection

# Inspecting a built model

## `Row`

Returned by `Model.row(name, **coords)`. One row as the assembled matrix
holds it, not a second walk of the expression, so what is shown is what
the solver receives.

| Field | Holds |
| --- | --- |
| `constraint` | the equation this row belongs to |
| `coordinate` | the row's own coordinate, per free dimension |
| `index` | the solver's own row number |
| `terms` | one `RowTerm` per coefficient |
| `sense`, `lower`, `upper` | read from the row's bounds |

| `RowTerm` field | Holds |
| --- | --- |
| `column` | the solver's own column number |
| `variable` | the variable that column belongs to |
| `coordinate` | that column's coordinate, per dimension |
| `coefficient` | the value in the matrix |

A variable occupies a contiguous range of the column space from its
`start`, so a column resolves to its variable by that range and to a
coordinate through the variable's own numbering rule.

`sense` is read from the bounds: equal bounds are `==`, an infinite lower
bound is `<=`, an infinite upper bound is `>=`.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, cost[P, W] * x[P, W]) <= supply[P])

print(m.row("supply", P="porto"))
```

Output:

```text
supply[P='porto']  row 1
  3·x[porto,berlin] + 1·x[porto,paris] + 6·x[porto,rome] <= 25
```

A coordinate at which the constraint has no row raises `ValueError`; the
message points to the function that reports why it is missing.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set

P = Set("P", np.array(["p1", "p2", "p3"]))
m = Model("m")
x = m.var("x", (P,), upper=5.0)
one = Param.from_dense("one", (P,), np.ones(3))
rhs = Param.from_long("rhs", (P,), {"P": np.array(["p1", "p2"])}, np.ones(2))
m.eq("cap", one[P] * x[P] <= rhs[P])

m.row("cap", P="p3")
```

Raises ValueError:

```text
ValueError: constraint 'cap' states no row at {'P': 'p3'}; `absent('cap')` names the rule that dropped it
```

## `Absence`

Returned by `Model.absent(name)`. What a constraint set out to produce,
what it produced, and which coordinates were dropped.

| Field | Holds |
| --- | --- |
| `constraint` | the equation this is about |
| `stated_by` | `"terms"` where the rows are derived, `"over"` where given explicitly |
| `expected`, `standing` | rows set out, rows kept |
| `dropped_rows` | one `DroppedRow(coordinate, rule, detail)` per row lost |
| `dropped_terms` | one `DroppedTerm(coordinate, variable, rule, detail)` per term lost |

`expected - len(dropped_rows) == standing`.

| `dropped_rows` rule | Meaning |
| --- | --- |
| `term-does-not-reach` | a term has no value at that coordinate, so the row would express something unwritten |
| `where` | the condition excludes it |
| `absent-rhs` | the right-hand side has no value there |

| `dropped_terms` rule | Meaning |
| --- | --- |
| `absent-coefficient` | a coefficient absent inside a sum, so the row stands with one term fewer |

The row and term split is the distinction to hold onto: a coefficient
absent inside a sum removes a **term** and leaves the row standing; a term
absent along a **free** dimension removes the **row**.

Under `over=` the rows are given explicitly, so nothing is dropped and a
right-hand side that misses one raises instead. An empty `dropped_rows`
beside `stated_by="over"` is structural.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["p1", "p2"]))
W = Set("W", np.array(["w1", "w2", "w3"]))
m = Model("t")
flow = m.var("flow", (P, W))
cost = Param.from_long(
    "cost",
    (P, W),
    {"P": np.array(["p1", "p1", "p2"]), "W": np.array(["w1", "w2", "w1"])},
    np.array([1.0, 2.0, 3.0]),
)
supply = Param.from_dense("supply", (P,), np.array([3.0, 3.0]))
m.eq("supply", Sum(W, cost[P, W] * flow[P, W]) <= supply[P])

print(m.absent("supply"))
```

Output:

```text
supply  2 of 2 rows  stated by terms
  term absent P='p1', W='w3'  flow  absent-coefficient (cost)
  term absent P='p2', W='w2'  flow  absent-coefficient (cost)
  term absent P='p2', W='w3'  flow  absent-coefficient (cost)
```

---

# /reference/model

# Model

## `Model`

```
Model(name="model", sense="min")
```

A model holds one column space, the constraints declared against it, and an
objective. `name` labels it and is otherwise unused. `sense` is `"min"` or
`"max"`, set once here; any other value raises.

| Member | Returns |
| --- | --- |
| `var(name, sets, subset=None, lower=0.0, upper=inf, integer=False)` | a `Variable` occupying the next range of columns |
| `eq(name, relation, where=None, over=None)` | a `Constraint` occupying the next range of rows |
| `set_objective(expression)` | nothing; sets the objective |
| `sense` | `"min"` or `"max"`, as declared |
| `solve(solver="highs", options=None)` | a `Solution` |
| `assemble()` | an `Assembled`: the matrix, with no solver involved |
| `n_columns`, `n_rows`, `nnz` | the shape declared so far |
| `column_bounds()` | the lower and upper bound vectors, in column order |
| `integrality()` | one flag per column |
| `objective_coefficients()` | one coefficient per column |
| `explain()` | an `Explanation` of what the model built |
| `to_yaml(inline=False)` | the text of this model's file, with its data inline where asked |
| `objective` | the objective expression, or `None` |

Declaring costs shapes, not blocks: `n_rows` and `nnz` are known when a
constraint is added, and no matrix exists until `assemble` or `solve`.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= 30.0)
m.eq("total", Sum(P, W, x[P, W]) <= 100.0)
m.set_objective(Sum(P, W, x[P, W]))

print(m.n_columns, m.n_rows, m.nnz)
print(m.objective_coefficients())
```

Output:

```text
6 3 12
[1. 1. 1. 1. 1. 1.]
```

## `Assembled`

The model's matrix in CSR form, returned by `assemble`. `indices` and
`values` are views of the one buffer the model allocated; only `indptr` is
built.

| Member | Returns |
| --- | --- |
| `indptr`, `indices`, `values` | the matrix in CSR form |
| `n_rows`, `n_cols` | its shape |
| `row_lower`, `row_upper` | one bound per row |
| `col_lower`, `col_upper`, `col_cost`, `integrality` | one entry per column |
| `row_of(name)` | a constraint's rows, as a slice |
| `to_dense()` | the matrix as an ndarray |

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= 30.0)

assembled = m.assemble()
print(assembled.n_rows, assembled.n_cols)
print(assembled.indptr)
print(assembled.row_of("supply"))
print(assembled.to_dense())
```

Output:

```text
2 6
[0 3 6]
slice(0, 2, None)
[[1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1.]]
```

`to_dense` is for a small model. A model of any size is read through
`row_of` and the CSR arrays.

## What a model built

`explain()` reports every declaration with the count it built, and has
`built=True`. It returns the same record type a `Definition` returns with
every count absent, so one reader serves both.

A model holds variables and constraints; its sets and parameters are
collected from them, in order of first appearance. A dimension a
coefficient introduces belongs to no variable and is found through the
parameter that has it.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, cost[P, W] * x[P, W]) <= supply[P])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

print(m.explain())
```

Output:

```text
transport  min  6 columns · 2 rows · 6 nonzeros
  sets        P 2 · W 3
  parameters  cost (P,W) 6 · supply (P) 2
  variables   x (P×W) 6 cols [0.0, inf]
  constraint  supply (P)  Sum(W, cost[P, W] * x[P, W]) <= supply[P]  2 rows  6 nz
  objective   min  Sum(P, W, cost[P, W] * x[P, W])
```

---

# /reference/param

# Param

## `Param`

Coefficients over a set product. A parameter is data, not a model object:
it has no columns and produces no rows. It supplies a term's coefficient
and a constraint's right-hand side.

A parameter's array declares `absence="empty"`. A coordinate it does not
have is a coefficient that is not there, which is the additive identity a
sum needs.

### `Param.from_dense(name, sets, values)`

Every cell of `values` as a coefficient. `values.shape` must equal the
sizes of `sets`, in order; a mismatch raises `ValueError` with both shapes.

```python
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
print(cost.dims, cost.nnz)
```

Output:

```text
('P', 'W') 6
```

### `Param.from_long(name, sets, columns, values)`

Coefficients from one label column per set and one value column. `columns`
is a mapping keyed by set name; each column and `values` are read in
parallel, so all have the same length.

```python
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_long(
    "cost",
    (P, W),
    {"P": np.array(["lisbon", "porto"]), "W": np.array(["berlin", "paris"])},
    np.array([2.0, 1.0]),
)
print(cost.nnz)
```

Output:

```text
2
```

A label column of a different length raises `ValueError`; the message
gives the parameter, the column, its length and the value column's.

```python raises=ValueError
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

Param.from_long(
    "cost",
    (P, W),
    {"P": np.array(["lisbon"]), "W": np.array(["berlin", "paris"])},
    np.array([2.0, 1.0]),
)
```

Raises ValueError:

```text
ValueError: parameter 'cost': label column 'P' has length 1 and the value column has length 2; they name the same entries
```

### Members

| Member | Returns |
| --- | --- |
| `dims` | the names of the sets it is indexed over |
| `nnz` | the number of coefficients |
| `materialise()` | the coefficients as a `nimblend` array |
| `param[sets]` | a reference, with the sets given checked against `dims` |

A label in place of a set fixes that dimension at one member: the
coefficients at that member are read and the dimension leaves the
reference.

```python
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
print(cost[P, W].dims)
print(cost[P, "berlin"].dims)
```

Output:

```text
('P', 'W')
('P',)
```

## `Coefficient`

What a term reads as its coefficient: a `name` to report, the `dims` it is
indexed over, the array it `materialise()`s to, and a reading at its sets.
A parameter read at its sets is one, and so is an arithmetic combination of
coefficients, so a function that reports a coefficient handles either
through one interface.

`+`, `-`, `*`, `/` and a power by a number combine coefficients. The
combination is symbolic: it holds references, derives its dimensions from
its operands, and is evaluated once, when the term it multiplies is
materialised. It can therefore be written in a definition before any data
exists.

```python
import numpy as np
from nimopt import Param, Set

G = Set("G", np.array(["base", "peak"]))
T = Set("T", np.arange(3))
price = Param.from_dense("fuel_price", (G, T), np.full((2, 3), 30.0))
eta = Param.from_dense("efficiency", (G, T), np.array([[0.5] * 3, [0.4] * 3]))

unit_cost = price[G, T] / eta[G, T]
print(unit_cost.name, unit_cost.dims)
print(unit_cost[G, T].materialise().to_dense()[:, 0])
```

Output:

```text
(fuel_price / efficiency) ('G', 'T')
[60. 75.]
```

A parameter has no arithmetic of its own. It is read at its sets, and the
references combine.

```python raises=TypeError
import numpy as np
from nimopt import Param, Set

G = Set("G", np.array(["a", "b"]))
price = Param.from_dense("price", (G,), np.array([1.0, 2.0]))
eta = Param.from_dense("eta", (G,), np.array([0.5, 0.4]))
price / eta
```

Raises TypeError:

```text
TypeError: parameter 'price' carries ('G',) and states no coefficient until it is read; read it at its sets as price[G]
```

---

# /reference/sets

# Sets and domains

## `Set`

```
Set(name, labels)
```

A named dimension with labels. `labels` is an array; `name` is what every
reference to the dimension uses.

| Member | Returns |
| --- | --- |
| `name`, `labels` | the declared name and labels |
| `len(set)` | the number of members |
| `position_of(labels)` | the position of each label given |
| `coord` | the coordinate the labels resolve through |
| `cyclic` | the same set, with a lag that wraps rather than drops |
| `set - 1` | the set lagged, dropping the members a lag runs off |

```python
import numpy as np
from nimopt import Set

T = Set("T", np.array(["t0", "t1", "t2"]))

print(T.name, len(T))
print(T.labels)
print(T.position_of(np.array(["t2", "t0"])))
```

Output:

```text
T 3
['t0' 't1' 't2']
[2 0]
```

## `Alias`

```
Alias(name, set)
```

A second name for a set, sharing its labels and its coordinate. A
parameter over a set and its alias is an ordinary two-dimensional array,
so a model relating a set to itself does so without declaring a second
set. No labels are copied: the alias uses the coordinate the set already
built.

```python
import numpy as np
from nimopt import Alias, Param, Set

N = Set("N", np.array(["a", "b"]))
M = Alias("M", N)

flow = Param.from_dense("flow", (N, M), np.array([[0.0, 1.0], [1.0, 0.0]]))
print(flow.dims)
print(flow.materialise().to_dense())
```

Output:

```text
('N', 'M')
[[0. 1.]
 [1. 0.]]
```

## `product`

```
product(sets)
```

Every member of a set product, as a domain. Passed to `over=`, it states
the rows of a constraint explicitly, for a constraint whose terms each
cover some of its rows.

## `subset`

```
subset(sets, columns)
```

The members of a set product a model uses, given by label. `columns` holds
one label column per set, keyed by the set's name, read in parallel: the
k-th entry of each column belongs to the same member. It is a list of
members, not a cross product.

## `subset_of`

```
subset_of(sets, index)
```

The same, given by position. Each column of `index` is one member. A
caller holding positions passes them directly rather than building labels
to resolve back.

```python
import numpy as np
from nimopt import Set, product, subset, subset_of

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

print(product((P, W)).size)
print(
    subset(
        (P, W), {"P": np.array(["lisbon", "porto"]), "W": np.array(["berlin", "paris"])}
    ).size
)
print(subset_of((P, W), np.array([[0, 1], [0, 1]])).size)
```

Output:

```text
6
2
2
```

The product has six members. Both subsets have two, `lisbon` with `berlin`
and `porto` with `paris`, because the columns are read in parallel.

---

# /reference/solution

# Solution

## `Solution`

Returned by `Model.solve`. Primal and dual values, returned over the sets
they were declared over.

| Member | Returns |
| --- | --- |
| `status` | the outcome the solver reported |
| `objective` | the optimal objective value |
| `primal(name)` | the named variable's values over its own sets |
| `dual(name)` | the named constraint's duals over its free sets |

`status` is readable whatever the solver reported. `objective`, `primal`
and `dual` are not: a model the solver did not bring to an optimum has no
answer, and a vector it left behind would be indistinguishable from one.
Read `status` first.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

print(solution.status)
print(solution.objective)
print(solution.primal("x").to_dense())
print(solution.dual("demand").to_dense())
```

Output:

```text
optimal
135.0
[[20.  0. 10.]
 [ 0. 15.  5.]]
[3. 1. 6.]
```

Reading a value from a model that did not reach an optimum raises
`ValueError`; the message gives the status.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("infeasible")
x = m.var("x", (P, W))
m.eq("floor", Sum(W, x[P, W]) >= 10.0)
m.eq("ceiling", Sum(W, x[P, W]) <= 1.0)
m.set_objective(Sum(P, W, x[P, W]))

m.solve().objective
```

Raises ValueError:

```text
ValueError: the model's status is 'infeasible', so it carries no objective; read `status` before reading values
```

## The array type of a value

A variable over a full product has a value at every cell of its frame, and
the solver returns them in column order, so they reshape into a
`DenseArray` with no index built at all. A variable over a subset has
values at its members alone, and a dense frame would be the grid it was
declared to avoid, so those stay a `SparseArray`. A dual follows its
constraint's rows by the same rule.

Every array declares `absence="unknown"`: a coordinate the model did not
have has no value, and combining two models' results must not invent a zero
for it.

---

# /reference/solvers

# Solvers

## `available` and `capabilities`

`available()` lists every adapter whose backend can be imported in the
current environment, with what each declares. `capabilities(name)` reports
for an adapter whether or not its backend is installed, because a
descriptor describes what the adapter does as shipped, and reading one is
how a caller decides what to install.

A descriptor describes the adapter, not the library behind it: a solver
feature the adapter does not call is `absent`.

What `available()` lists depends on the machine; what `capabilities(name)`
reports does not.

```python
from nimopt import available, capabilities

print("highs" in available())
print(capabilities("highs"))
print(capabilities("gurobi"))
print(capabilities("mosek"))
```

Output:

```text
True
highs  integrality native · duals native · conflict native · ray native  refuses duals+integrality
gurobi  integrality native · duals native · conflict native · ray native  refuses duals+integrality
mosek  integrality native · duals native · conflict absent · ray native  refuses duals+integrality
```

## `Capabilities`

| Member | Returns |
| --- | --- |
| `solver` | the adapter's name |
| `support` | one of `"native"` or `"absent"` per capability |
| `refused` | the pairs this adapter refuses together |
| `supports(capability)` | whether the adapter handles it at all |
| `refuses(one, other)` | whether it refuses the two together |

The capabilities are `integrality`, `duals`, `conflict` and `ray`. A flat
set is insufficient: a solver can support two and refuse their
combination. Every adapter refuses `integrality` with `duals`, because a
mixed-integer model's duals are not the relaxation's, so a model with
integer columns has no duals at all and `Solution.dual` raises.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set, Sum

T = Set("T", np.arange(2))
one = Param.from_dense("one", (T,), np.ones(2))

m = Model("m")
x = m.var("x", (T,), upper=3.0, integer=True)
m.eq("cap", one[T] * x[T] <= 2.0)
m.set_objective(Sum(T, one[T] * x[T]))

m.solve().dual("cap")
```

Raises ValueError:

```text
ValueError: this model carries integer columns and 'highs' refuses duals for a model with integrality, so there is no dual for constraint 'cap' to read: a mixed-integer model's duals are not its relaxation's
```

## `Session`

Returned by `Model.session(solver="highs", options=None)`. One solver's
model, opened on one assembled model and kept: a solve hands the matrix
across, and a question asked afterwards is asked of the same solved
instance.

| Member | Returns |
| --- | --- |
| `assembled` | the matrix the session was opened on |
| `solver`, `capabilities` | which adapter, and what it can do |
| `status` | what the last solve reported, or `None` before one |
| `solve()` | a `Solution` |
| `close()` | releases the solver's model |

`Model.solve()` opens a session, solves and closes it, so a caller who
wants only a solution needs no session.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

SNAP = Set("snapshot", np.arange(3))
GEN = Set("generator", np.array(["wind", "gas"]))
p_max = Param.from_dense("p_max", (GEN,), np.array([10.0, 20.0]))
load = Param.from_dense("load", (SNAP,), np.array([25.0, 20.0, 5.0]))
cost = Param.from_dense("cost", (GEN,), np.array([1.0, 5.0]))

m = Model("dispatch", sense="min")
p = m.var("p", (SNAP, GEN), lower=0.0, upper=p_max)
m.eq("balance", Sum(GEN, p[SNAP, GEN]) == load[SNAP])
m.set_objective(Sum(SNAP, GEN, cost[GEN] * p[SNAP, GEN]))

with m.session() as session:
    solution = session.solve()
    print(solution.status, solution.objective)
    print(session.status)
```

Output:

```text
optimal 150.0
optimal
```

## `Diagnosis`

Returned by `Session.diagnose()`. Why a model did not solve, as the rows
and columns that explain it.

| Member | Returns |
| --- | --- |
| `status`, `solver` | what the solve reported, and which adapter |
| `method` | `"native"` where the solver computed the conflict |
| `conflict` | one `Row` per conflicting row, or `None` where the model is not infeasible |
| `columns` | one `ColumnBound` per column of the conflict, with the bounds the model declares |
| `ray` | one `RayTerm` per column the ray moves, or `None` where the model is not unbounded |

`conflict` holds the same `Row` that `Model.row` returns, so a conflicting
row reads in one format. Two solvers may return different irreducible
sets: what holds of each is that removing it makes the model feasible, not
that the two agree.

A conflict is a question asked of the backend that solved, so the session
is what keeps it available.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

SNAP = Set("snapshot", np.arange(3))
GEN = Set("generator", np.array(["wind", "gas"]))
p_max = Param.from_dense("p_max", (GEN,), np.array([10.0, 20.0]))
load = Param.from_dense("load", (SNAP,), np.array([25.0, 100.0, 5.0]))
cost = Param.from_dense("cost", (GEN,), np.array([1.0, 5.0]))

m = Model("dispatch", sense="min")
p = m.var("p", (SNAP, GEN), lower=0.0, upper=p_max)
m.eq("balance", Sum(GEN, p[SNAP, GEN]) == load[SNAP])
m.set_objective(Sum(SNAP, GEN, cost[GEN] * p[SNAP, GEN]))

with m.session() as session:
    print(session.solve().status)
    print(session.diagnose())
```

Output:

```text
infeasible
infeasible  highs  conflict native
balance[snapshot=1]  row 1
  1·p[1,wind] + 1·p[1,gas] == 100
  bound  p[snapshot=1, generator='wind']  [0, 10]
  bound  p[snapshot=1, generator='gas']  [0, 20]
```

HiGHS computes its conflict over the model's linear relaxation. A model
that is feasible as an LP and infeasible only through its integrality
therefore yields no conflict, and the adapter raises rather than naming
rows it did not prove. Gurobi's conflict covers the integrality. Mosek's
adapter computes no conflict, so a session on it refuses the question and
names `capabilities("mosek")` as what states so.

## `options` and `Option`

`options()` lists every option a caller can set, in `nimopt`'s own names. An
option outside the list raises rather than being passed to a solver that
would ignore it, so a misspelled name stops a solve instead of running a
different one.

`options(solver)` lists the same options with that solver's own name and
values, which is how a caller follows one into the solver's own
documentation.

```python
from nimopt import options

for option in options("highs"):
    print(f"{option.name:<16} {option.native}")
```

Output:

```text
time_limit       time_limit
iteration_limit  simplex_iteration_limit
node_limit       mip_max_nodes
mip_gap          mip_rel_gap
mip_abs_gap      mip_abs_gap
feasibility_tol  primal_feasibility_tolerance
optimality_tol   dual_feasibility_tolerance
threads          threads
seed             random_seed
log              output_flag
presolve         presolve
method           solver
newton_system    hipo_system
crossover        run_crossover
pdlp_tol         pdlp_optimality_tolerance
```

An `Option` has a `name`, the `kind` it takes, what it `does`, and its
`choices` where it takes one of a set. Read for a solver it also has
`native` and `native_choices`.

| Option | Takes | Does | `highs` | `gurobi` | `mosek` |
| --- | --- | --- | --- | --- | --- |
| `time_limit` | float | seconds the solver may run for | `time_limit` | `TimeLimit` | `optimizer_max_time` |
| `iteration_limit` | int | simplex iterations the solver may take | `simplex_iteration_limit` | `IterationLimit` | `sim_max_iterations` |
| `node_limit` | int | branch-and-bound nodes the solver may explore | `mip_max_nodes` | `NodeLimit` | `mio_max_num_branches` |
| `mip_gap` | float | relative gap at which a mixed-integer solve stops | `mip_rel_gap` | `MIPGap` | `mio_tol_rel_gap` |
| `mip_abs_gap` | float | absolute gap at which a mixed-integer solve stops | `mip_abs_gap` | `MIPGapAbs` | `mio_tol_abs_gap` |
| `feasibility_tol` | float | how far a primal solution may miss a row | `primal_feasibility_tolerance` | `FeasibilityTol` | `basis_tol_x` |
| `optimality_tol` | float | how far a dual solution may miss a bound | `dual_feasibility_tolerance` | `OptimalityTol` | `basis_tol_s` |
| `threads` | int | threads the solver may use; 0 leaves it the choice | `threads` | `Threads` | `num_threads` |
| `seed` | int | the seed the solver randomises from | `random_seed` | `Seed` | `mio_seed` |
| `log` | bool | whether the solver writes its own iteration log | `output_flag` | `OutputFlag` | `log` |
| `presolve` | `off` / `choose` / `on` | how hard the solver presolves | `presolve` | `Presolve` | `presolve_use` |
| `method` | `choose` / `simplex` / `barrier` / `hipo` / `pdlp` | the algorithm the solver runs | `solver` | `Method` | `optimizer` |
| `newton_system` | `choose` / `augmented` / `normaleq` | the Newton system an interior point method factorises | `hipo_system` | not carried | not carried |
| `crossover` | `choose` / `off` / `on` | whether an interior point is moved to a vertex after the solve | `run_crossover` | `Crossover` | `intpnt_basis` |
| `pdlp_tol` | float | relative tolerance at which the first-order method stops | `pdlp_optimality_tolerance` | not carried | not carried |

A choice each solver spells differently is written once and translated, so
the value a caller writes means one thing whichever solver reads it. Not
every solver carries every option or every choice: `newton_system` and
`pdlp_tol` are HiGHS's, as are `hipo` and `pdlp` under `method`, and asking
Gurobi or Mosek for one of them is refused by name rather than answered by
a different algorithm. Mosek runs only its mixed-integer optimizer on a
model with integer columns, so `method` stays at `choose` there and any
other choice is refused naming that cause. What each method holds in
memory, and how to install a HiGHS that carries HiPO and a GPU, is in the
guide on [interior point and first-order methods](/guides/highs-methods).

## Progress reporting

`build()`, `assemble()`, `session()` and `solve()` take `progress=`.
`progress=True` draws a report in a terminal and nothing where output is
redirected. A reporter of your own is used as given, so a notebook or an
interface writes there: it implements `start(total, what)`, `step(done,
what)` and `done()`, and that is the whole contract.

The report covers building. Its resolution is the model's own structure:
the measuring pass counts constraints and the writing pass counts nonzeros,
so a model with one constraint of one term reports one step and no
fraction.

A solver's own account of a solve is the solver's to give, and `log=True`
asks for it. Building finishes before a solver starts, so the report and
the log never interleave.

---

# /reference/variable

# Variable

## `Variable`

Returned by `Model.var`. A variable over a set product, or over a subset
of one.

```
Model.var(name, sets, subset=None, lower=0.0, upper=inf, integer=False)
```

| Argument | Meaning |
| --- | --- |
| `name` | the name `Solution.primal` reads it back by |
| `sets` | the dimensions it is declared over |
| `subset` | the members it has; the full product when omitted |
| `lower`, `upper` | the bound every one of its columns takes |
| `integer` | whether its columns are integral |

The variable's columns are a virtual coordinate: a member's column is
computed from its multi-index by stride arithmetic for a full product, or
is its rank among a subset's codes. Nothing stores a column index, which is
why a variable over millions of columns costs only its members.

| Member | Returns |
| --- | --- |
| `dims` | the names of the sets it is over |
| `n_columns` | the number of columns it occupies |
| `domain()` | the members it has |
| `terms()` | its coefficients over `(*dims, COLUMN)` |
| `variable[sets]` | a one-term expression referencing it |

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))
open_plant = m.var("open_plant", (P,), lower=0.0, upper=1.0, integer=True)

print(x.dims, x.n_columns)
print(open_plant.n_columns)
print(m.n_columns)
print(m.integrality())
```

Output:

```text
('P', 'W') 6
2
8
[0 0 0 0 0 0 1 1]
```

Each variable occupies the next range of the model's one column space, so
`m.n_columns` counts every column declared so far.

## A variable over no dimension

A variable's bracket lists the dimensions it carries, so a variable over
none carries no bracket and enters a row on its own. It is one column: a
value-at-risk level, a budget slack, a bound every row of a family shares.
`theta[()]` is the same term written out.

A variable that does carry dimensions states no term until it is read, and
using one bare raises `TypeError` naming the reading it wants. The same rule
holds for a parameter, which is read `cost[G, T]` and, over no dimension,
`k`.

Comparing a variable states a row, so `==` between two variables states one
too rather than answering true or false. A list of variables therefore cannot
be searched with `in` or `.index`, which compare their items: those raise the
reading refusal, naming whichever variable they reached first. Keep variables
in a dict or a set, which match on identity, or search them by `name`.

```python
import numpy as np
from nimopt import Model, Set, Sum

S = Set("S", np.array(["s1", "s2"]))

m = Model("cvar", sense="min")
theta = m.var("theta", (), lower=-np.inf)
p = m.var("p", (S,))

m.eq("tail", theta - Sum(S, p[S]) >= 0.0)
m.set_objective(theta)
print(m.n_columns, m.n_rows)
print(m.constraints["tail"].relation)
```

Output:

```text
3 1
theta - Sum(S, p[S]) >= 0
```

## `COLUMN` and `ROW`

The dimension names `nimopt` reserves. `COLUMN` is `"__column__"` and `ROW`
is `"__row__"`; both are spelled so that no ordinary set name collides
with them.

A variable's terms are an array over `(*dims, COLUMN)`, and a constraint's
block is one over `(ROW, COLUMN)`. That is the whole of the correspondence
between a model and its matrix: the column space is a dimension, so the
array is the matrix.

```python
import numpy as np
from nimopt import COLUMN, ROW, Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

print(COLUMN, ROW)
print(x.terms().dims)
print(x.domain().dims)
```

Output:

```text
__column__ __row__
('P', 'W', '__column__')
('P', 'W')
```

A caller writes neither name. They exist to be recognised when a `nimblend`
array from inside a model is inspected.

---

# /nimblend/arrays

# nimblend arrays

`nimblend` is the layer below `nimopt`. It knows dimensions, labels, entries
and alignment, and nothing about optimization. A `nimopt` caller meets these
names when reading a solution or inspecting what a model built.

**Import from `nimblend` itself, never from a submodule, and never read an
array's `.index` or `.data` or a domain's `.codes`.** Those are the raw
index matrix, value buffer and ravelled members of the layer below the
array layer; reading them bypasses the contract. A test in this repository
fails on any of them.

Nor assemble one. `SparseArray.from_canonical` below takes an index matrix
the caller built, which is the array layer's own work: `nimopt`'s modules
call it nowhere, and a test holds them to that. A `Domain` returns the
array over its own members instead.

## `Array`

A labeled N-dimensional array. `Array` is the contract both
implementations satisfy, and what a caller writes against.

`absence` declares what a coordinate the array does not have means:
`"empty"` that it contributes nothing, `"unknown"` that it was not
modelled. Operators and reductions follow from that declaration, so it is
part of the array's meaning and not a hint.

| Member | Returns |
| --- | --- |
| `dims`, `shape`, `nnz` | the dimensions, their sizes, and the number of entries |
| `coords` | each dimension's coordinate, which resolves its labels |
| `absence` | the meaning of a coordinate the array does not have |
| `as_empty()`, `as_unknown()` | the array under the other absence declaration |
| `values()`, `coordinates()` | the entries and their multi-indices |
| `domain()` | the coordinates the array has |
| `sum`, `min`, `max`, `mean` | reductions over named dimensions |
| `sel`, `restrict` | a selection by label, and by domain |
| `rename`, `transpose`, `expand`, `conform` | reshaping the dimensions |
| `shift`, `roll` | a lag that drops, and one that wraps |
| `group` | entries combined into a destination |
| `to_dense(fill=None)` | the entries as an ndarray |
| `+`, `-`, `*`, `/`, unary `-` | arithmetic over one frame, and with a scalar |

The arithmetic is part of the contract, not an implementation's own: a
caller writing `coefficient * columns` is writing against `Array`. Both
implementations behave alike, including over frames that differ. One frame
nested inside the other broadcasts over the wider; frames sharing some
dimensions align on those and multiply out the rest. Frames sharing no
dimension raise. A product mixing the two implementations returns a
`SparseArray`, because a product intersects presence and so has at most
what the sparse operand has.

Absence and zero stay distinct. A stored `0.0` is a coordinate that is
present with value zero, which is not the same as one the array does not
have.

```python
import numpy as np
import nimblend as nb

labels = {"A": np.array(["a0", "a1"]), "B": np.array(["b0", "b1"])}
array = nb.SparseArray.from_dense(np.array([[1.0, 0.0], [0.0, 2.0]]), labels)

print(isinstance(array, nb.Array))
print(array.dims, array.shape, array.absence)
print(array.nnz)
print(array.sum("B").to_dense())
```

Output:

```text
True
('A', 'B') (2, 2) empty
4
[1. 2.]
```

`from_dense` stores every cell it was given, so this array has four entries
and not two: the zeros are stored, and stored means present.

## `SparseArray`

Entries in canonical order, under a coordinate per dimension. An entry that
is not stored is absent.

It is what an operation returns when the result is sparse, and what
`Solution.primal` returns for a variable over a subset: a dense frame there
would be the grid the variable was declared to avoid.

| Member | Returns |
| --- | --- |
| `SparseArray.from_dense(values, labels)` | every cell of an ndarray |
| `SparseArray.from_canonical(index, values, coords, dims)` | entries already in order |
| `as_empty()`, `as_unknown()` | the same entries under the other declaration |
| `to_csr()` | the entries as compressed rows |

```python
import numpy as np
import nimblend as nb

labels = {"A": np.array(["a0", "a1"]), "B": np.array(["b0", "b1"])}
array = nb.SparseArray.from_dense(np.array([[1.0, 0.0], [0.0, 2.0]]), labels)

print(array.absence)
print(array.as_unknown().absence)
print(array.to_dense())
```

Output:

```text
empty
unknown
[[1. 0.]
 [0. 2.]]
```

## `DenseArray`

An ndarray over labeled dimensions, distinguishing absence from zero.

How presence is stored follows the absence declaration, because the two
declarations need opposite things from an operator. An `"unknown"` array
tags absence with NaN, which propagates through arithmetic at no cost and
needs no storage beside the values. An `"empty"` array carries a boolean
mask, because absence is the additive identity there and substituting it is
cheaper than tagging.

`Solution.primal` returns one of these for a variable over a full product:
the solver returns a value at every cell of the frame in column order, so
they reshape with no index built at all.

```python
import numpy as np
import nimblend as nb

coords = {
    "A": nb.StoredCoord(np.array(["a0", "a1"])),
    "B": nb.StoredCoord(np.array(["b0", "b1"])),
}
array = nb.DenseArray(np.array([[1.0, 0.0], [0.0, 2.0]]), coords, ("A", "B"))

print(array.absence, array.nnz)
print(array.present)
print(array.as_unknown().absence)
```

Output:

```text
empty 4
[[ True  True]
 [ True  True]]
unknown
```

## Building an array from columns

`nimblend` itself provides the two constructors a caller uses when the data is
not already an ndarray.

| Constructor | Returns |
| --- | --- |
| `from_long(dims, coords, labels, values)` | one label column per dimension and one value column |
| `from_dense(values, labels)` | every cell of an ndarray |
| `is_canonical(index, shape)` | whether buffers are in the order `from_canonical` takes |

`from_long` resolves each label through the coordinate that dimension
already has, so a caller holding coordinates, which is any caller with sets
of its own, gives its entries as labels rather than resolving them to
positions first. The columns are read in parallel, so they must have equal
length.

```python
import numpy as np
import nimblend as nb

coords = {
    "t": nb.StoredCoord(np.array([2030, 2040])),
    "r": nb.StoredCoord(np.array(["DE", "FR"])),
}
arr = nb.from_long(
    ("t", "r"),
    coords,
    {"t": np.array([2030, 2040, 2040]), "r": np.array(["DE", "DE", "FR"])},
    np.array([5.0, 6.0, 7.0]),
)
print(arr.nnz, arr.to_dense()[1, 1])
```

Output:

```text
3 7.0
```

A column of a different length raises `ValueError`; the message gives the
column and both lengths.

```python raises=ValueError
import numpy as np
import nimblend

nimblend.from_long(
    ("t",),
    {"t": nimblend.StoredCoord(np.array([2030, 2040]))},
    {"t": np.array([2030, 2040])},
    np.array([1.0]),
)
```

Raises ValueError:

```text
ValueError: label column 't' has length 2 and the value column has length 1; they name the same entries
```

`is_canonical` answers the question `SparseArray.from_canonical` asks a
caller to answer about buffers the caller built: entries sorted by ravel
key with no repeat. Verifying it inside `from_canonical` would cost the
ravel that path exists to avoid.

```python
import numpy as np
import nimblend as nb

ordered = np.array([[0, 0, 1], [0, 1, 0]], dtype=np.int32)
print(nb.is_canonical(ordered, (2, 2)))
print(nb.is_canonical(ordered[:, ::-1].copy(), (2, 2)))
```

Output:

```text
True
False
```

## The frame of a binary result

`combined_dims(left, right)` returns the dimensions a binary operator's
result has, from the two operands' dimensions alone. A caller reads it
before materialising either operand, which is what lets a combination
report its frame while its data is still unbound.

| Operands | Result |
| --- | --- |
| equal frames | that frame, in its order |
| one frame nested in the other | the wider |
| frames that overlap | the left, then the dimensions only the right has |
| frames sharing no dimension | raises |

```python
import nimblend as nb

print(nb.combined_dims(("P", "Q"), ("Q", "R")))
print(nb.combined_dims(("P",), ("P", "Q")))
```

Output:

```text
('P', 'Q', 'R')
('P', 'Q')
```

Frames sharing no dimension have nothing to align on, so `combined_dims`
raises: their combination would be an outer product no caller asked for.

```python raises=ValueError
import nimblend

nimblend.combined_dims(("P",), ("Q",))
```

Raises ValueError:

```text
ValueError: frames ('P',) and ('Q',) share no dimension; there is nothing to align them on
```

## Densifying an unknown array

An array declaring `"unknown"` that does not have every coordinate of its
frame raises on `to_dense()` without a fill. There is no value it can place
at the rest, and choosing one silently would invent a value.

```python raises=ValueError
import numpy as np
import nimblend

labels = {"A": np.array(["a0", "a1"]), "B": np.array(["b0", "b1"])}
partial = nimblend.SparseArray.from_canonical(
    np.array([[0], [0]], dtype=np.int32),
    np.array([1.0]),
    {
        "A": nimblend.StoredCoord(labels["A"]),
        "B": nimblend.StoredCoord(labels["B"]),
    },
    ("A", "B"),
    absence="unknown",
)

partial.to_dense()
```

Raises ValueError:

```text
ValueError: this array declares absence 'unknown' and does not carry every coordinate of its frame, so densifying must state fill=<value> to place at the rest
```

With a fill value, the grid is returned.

```python
import numpy as np
import nimblend as nb

labels = {"A": np.array(["a0", "a1"]), "B": np.array(["b0", "b1"])}
partial = nb.SparseArray.from_canonical(
    np.array([[0], [0]], dtype=np.int32),
    np.array([1.0]),
    {
        "A": nb.StoredCoord(labels["A"]),
        "B": nb.StoredCoord(labels["B"]),
    },
    ("A", "B"),
    absence="unknown",
)

print(partial.to_dense(fill=np.nan))
```

Output:

```text
[[ 1. nan]
 [nan nan]]
```

---

# /nimblend/domains

# nimblend domains

## `Domain`

A sorted, unique set of multi-indices over named dimensions.

A coordinate resolves labels for one dimension; a domain does so for a
tuple of them: which multi-indices it has, what position each occupies, and
which multi-index sits at a position. It records which coordinates are
present and nothing about where they are numbered from.

| Constructor | Returns |
| --- | --- |
| `Domain.full(dims, coords)` | every coordinate of the product `dims` spans |
| `Domain.from_labels(dims, coords, labels)` | a domain from one label column per dimension |
| `Domain.from_coordinates(dims, coords, index)` | a domain from an index matrix of one row per dimension |

| Member | Returns |
| --- | --- |
| `size`, `dims`, `shape`, `coords` | the number of members, the dimensions, their extents and their coordinates |
| `is_full` | whether every coordinate of the product is present |
| `coordinates()` | the multi-index of each member, as an int32 index matrix |
| `labels()` | each member's label, per dimension |
| `intersect(other)` | the members both have |
| `union(other)` | the members either has |
| `difference(other)` | the members this one has and `other` does not |
| `positions_of(array)` | each entry of `array` as its position here, `-1` where absent |
| `positions_of_coordinates(index)` | each column of an index matrix as its position here, `-1` where absent |
| `expand(dims, coords)` | every member crossed with the full extent of the named dimensions |
| `transpose(*dims)` | the same members, over the dimensions in the order given |
| `as_coord(start=0)` | the domain read as a coordinate, its members numbered from `start` |
| `array(values, absence="empty")` | the members with one value each, as a `SparseArray` |
| `identity(into, coord, start=0)` | each member paired with its own position along `into`, valued 1.0 |

That table is the whole surface. **A domain's `codes` are the raw ravelled
members of the layer below it, as an array's `.index` and `.data` are its
raw buffers; never read them.** `coordinates()` and `labels()` report which
members are present, `positions_of_coordinates` reports where one sits,
`as_coord` numbers them, and `array` and `identity` return an array over
them, so nothing above needs to build an index matrix either. A test in
this repository fails on a read of any of the three.

The label columns of `from_labels` are read in parallel: the k-th entry of
each column belongs to the same member. A domain is a list of members, not
a cross product.

```python
import numpy as np
import nimblend as nb

coords = {
    "P": nb.StoredCoord(np.array(["lisbon", "porto"])),
    "W": nb.StoredCoord(np.array(["berlin", "paris", "rome"])),
}

full = nb.Domain.full(("P", "W"), coords)
pairs = nb.Domain.from_labels(
    ("P", "W"),
    coords,
    {"P": np.array(["lisbon", "porto"]), "W": np.array(["berlin", "paris"])},
)

print(full.size, pairs.size)
print(pairs.coordinates())
print(pairs.labels())
print(full.intersect(pairs).size, full.difference(pairs).size)
```

Output:

```text
6 2
[[0 1]
 [0 1]]
{'P': array(['lisbon', 'porto'], dtype='= 0`.

`as_coord` reads the domain as a coordinate: a member's position is its
rank among the members present, numbered from `start`. This is how a
dimension spanning a subset of a product is numbered.

```python
import numpy as np
import nimblend as nb

coords = {
    "P": nb.StoredCoord(np.array(["lisbon", "porto"])),
    "W": nb.StoredCoord(np.array(["berlin", "paris", "rome"])),
}
pairs = nb.Domain.from_labels(
    ("P", "W"),
    coords,
    {"P": np.array(["lisbon", "porto"]), "W": np.array(["berlin", "paris"])},
)

asked = np.array([[0, 1, 1], [0, 1, 2]], dtype=np.int32)
print(pairs.positions_of_coordinates(asked))
print(pairs.positions_of_coordinates(asked) >= 0)

numbered = pairs.as_coord(100)
print(numbered.to_position(asked[:, :2]))
```

Output:

```text
[ 0  1 -1]
[ True  True False]
[100 101]
```

`("porto", "rome")` is not a member, so it returns `-1`. The other two are
the domain's first and second members, and `as_coord(100)` numbers them
from 100.

## Crossing a domain with further dimensions

`expand` replicates every member across the full extent of the named
dimensions, which is how the coordinates a term *could* have are enumerated
before asking which of them it does. The new dimensions are appended;
`transpose` reads the result in another order. A member's code is its own
scaled by the appended extent, plus each position within it, so the cross
product is arithmetic on the members and no index matrix is built to hold
it.

```python
import numpy as np
import nimblend as nb

coords = {
    "P": nb.StoredCoord(np.array(["lisbon", "porto"])),
    "W": nb.StoredCoord(np.array(["berlin", "paris", "rome"])),
    "H": nb.StoredCoord(np.array([0, 1])),
}
pairs = nb.Domain.from_labels(
    ("P", "W"),
    coords,
    {"P": np.array(["lisbon", "porto"]), "W": np.array(["berlin", "paris"])},
)

hourly = pairs.expand(("H",), coords)
print(hourly.dims, hourly.size)
print(hourly.labels())
print(hourly.transpose("H", "P", "W").dims)
```

Output:

```text
('P', 'W', 'H') 4
{'P': array(['lisbon', 'lisbon', 'porto', 'porto'], dtype='<U6'), 'W': array(['berlin', 'berlin', 'paris', 'paris'], dtype='<U6'), 'H': array([0, 1, 0, 1])}
('H', 'P', 'W')
```

Two members crossed with two hours are four, and `transpose` presents them
over the dimensions in another order without changing which members are
present.

## A domain returns an array

A caller holding one value per member, or wanting each member paired with
its own position along a new dimension, asks the domain rather than
building an index matrix. Both are readers standing above the raw members,
as `coordinates()` and `as_coord()` are.

`array(values)` assigns one value to each member, in the order they are
held. The members ascend, so the entries are canonical as written and no
sort runs.

```python
import numpy as np
import nimblend as nb

coords = {"t": nb.StoredCoord(np.array([2030, 2040, 2050]))}
members = nb.Domain.full(("t",), coords)
print(members.array(np.array([1.0, 2.0, 3.0])).values())
```

Output:

```text
[1. 2. 3.]
```

A full domain's members ascend with the ravel key, which is the order
`values.ravel()` reads a grid in, so a whole array is built in one call.

```python
import numpy as np
import nimblend as nb

coords = {
    "x": nb.StoredCoord(np.array(["a", "b"])),
    "y": nb.StoredCoord(np.array([10, 20, 30])),
}
values = np.arange(6, dtype=np.float64).reshape(2, 3)
grid = nb.Domain.full(("x", "y"), coords).array(values.ravel())
print(grid.to_dense())
```

Output:

```text
[[0. 1. 2.]
 [3. 4. 5.]]
```

One value per member is the whole rule. A column of another length raises
`ValueError`.

```python raises=ValueError
import numpy as np
import nimblend

coords = {"t": nimblend.StoredCoord(np.array([2030, 2040, 2050]))}
nimblend.Domain.full(("t",), coords).array(np.array([1.0, 2.0]))
```

Raises ValueError:

```text
ValueError: a domain of 3 member(s) takes one value each, as a column of that length; got shape (2,)
```

`identity(into, coord, start)` pairs each member with its own position
along a new dimension, valued 1.0. A member's position is its rank plus
`start`, which is the numbering `as_coord(start)` uses, so an array built
one way and a coordinate built the other place a member alike. `coord` is
the coordinate of the new dimension and spans the whole extent the
positions are numbered into, wider than these members where several
domains share one numbering.

```python
import numpy as np
import nimblend as nb

coords = {"t": nb.StoredCoord(np.array([2030, 2040, 2050]))}
members = nb.Domain.full(("t",), coords)
paired = members.identity("k", nb.ProductCoord((20,)), start=10)
print(paired.dims)
print(paired.coordinates())
```

Output:

```text
('t', 'k')
[[ 0  1  2]
 [10 11 12]]
```

A destination too short for the members it is asked to number raises
`ValueError` rather than writing a position outside it.

```python raises=ValueError
import numpy as np
import nimblend

coords = {"t": nimblend.StoredCoord(np.array([2030, 2040, 2050]))}
members = nimblend.Domain.full(("t",), coords)
members.identity("k", nimblend.ProductCoord((6,)), start=4)
```

Raises ValueError:

```text
ValueError: 3 member(s) numbered from 4 reach position 6, and dimension 'k' spans 6
```

## The three coordinates

A coordinate resolves where a label sits along one dimension. Which of the
three is used follows from what the dimension is.

| Coordinate | Is | Used for |
| --- | --- | --- |
| `StoredCoord(labels)` | labels held as an array | a dimension whose members are named |
| `ProductCoord(sizes, start=0)` | positions of a full product, numbered from `start` | a dimension whose positions are computed, such as a variable's columns |
| `SubsetCoord(codes, sizes, start=0)` | positions of a subset of a product, numbered from `start` in code order | a variable over a subset, where a position is a rank among the codes |

`SubsetCoord` gives an entry's position as its rank among the codes, so a
block already in canonical order needs no lookup at all.

```python
import numpy as np
import nimblend as nb

stored = nb.StoredCoord(np.array(["a", "b", "c"]))
print(stored.to_position(np.array(["c", "a"])))

product = nb.ProductCoord((2, 3))
print(product.to_position(np.array([[0, 1], [2, 0]])))

subset = nb.SubsetCoord(np.array([0, 4]), (2, 3))
print(subset.to_position(np.array([[0, 1], [0, 1]])))
```

Output:

```text
[2 0]
[2 3]
[0 1]
```

`StoredCoord` looks a label up among the ones it holds, so `"c"` resolves
to position 2. `ProductCoord` ravels a multi-index against the sizes, so
`(0, 2)` is position 2 and `(1, 0)` is position 3. `SubsetCoord` holds the
codes `0` and `4`, which are those same two members, and returns their
ranks.

`to_position` takes an index matrix of one row per dimension, and each
column is one entry.

## `EntryBuffer`

A fixed index and value buffer handing out successive slices.

A block computed into a reserved slice never exists as a separate object,
so assembling several of them holds one copy of the result rather than one
copy per block plus the result. It is the destination a model's assembly
writes into: each constraint writes its rows into its own slice of one
buffer.

| Member | Returns |
| --- | --- |
| `EntryBuffer(ndim, capacity)` | a buffer for `capacity` entries of `ndim` dimensions |
| `reserve(n)` | the next `n` index and value slices, to write into |
| `written()` | the index and values written so far |
| `array(coords, dims, absence="empty")` | what was written, as an array |

```python
import numpy as np
import nimblend as nb

buffer = nb.EntryBuffer(2, 4)
index, values = buffer.reserve(2)
index[:] = np.array([[0, 1], [0, 1]])
values[:] = np.array([5.0, 6.0])

coords = {
    "A": nb.StoredCoord(np.array(["a0", "a1"])),
    "B": nb.StoredCoord(np.array(["b0", "b1"])),
}
array = buffer.array(coords, ("A", "B"))

print(array.nnz)
print(array.to_dense())
```

Output:

```text
2
[[5. 0.]
 [0. 6.]]
```

The slices `reserve` hands out are views of the one allocation, so writing
into them is writing into the array that comes out.

---

# /nimblend

# nimblend

`nimblend` is a labeled sparse N-dimensional array library. Its vocabulary is
dimensions, labels, entries and alignment, and it knows nothing about
optimization. It depends on NumPy and nothing else.

`nimopt` imports it, never the reverse. A model reaches `nimblend` in two
places: a solution is returned as a `nimblend` array, and a constraint's rows
are a `nimblend` domain. A reader who wants labeled sparse data and no model
at all can use it on its own.

## Design

**Absence is distinct from zero.** An entry is either stored or absent, and
every array declares the meaning of absence: `"empty"` for a coordinate
that contributes nothing, `"unknown"` for one that was never modelled.
Division by an absent value raises instead of producing infinity.

**One contract, two implementations.** `Array` defines what an array does.
`SparseArray` stores only the entries it has; `DenseArray` stores a grid and
the presence its declaration implies. Both are tested against the same
conformance suite.

**A coordinate is computed, not stored.** A dimension spanning millions of
positions costs nothing to hold: `ProductCoord` computes a position by
stride arithmetic and `SubsetCoord` by rank among the members of a domain.

**A domain is a set of coordinates.** It reports which members it has and
where each sits, and through `array` and `identity` it returns an array
over them, so a caller never assembles an index matrix.

## Pages

- [Arrays](/nimblend/arrays): the contract, the two implementations, the
  constructors, and what an absence declaration means.
- [Domains](/nimblend/domains): the coordinates an array has, the three ways a
  position is computed, and the buffer assembly writes into.

---

# /explanation/a-variable-is-a-dimension

# A variable is a dimension

A model has one column space. `m.var` does not create an object with its
own numbering; it takes the next range of that space, and each subsequent
variable continues from where the previous one ended.

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("transport")
x = m.var("x", (P, W))
y = m.var("y", (P,))

print(x.n_columns, y.n_columns)
print(m.n_columns)
```

Output:

```text
6 2
8
```

The shared column space is what makes the column a dimension. A variable's
coefficients form an array over `(*dims, COLUMN)`, where `COLUMN` is the
model's column space. A variable does not own columns; it occupies a block
of one dimension that all variables share.

## A column is computed, not stored

A member's column is a virtual coordinate, obtained by arithmetic rather
than by lookup.

For a full product, the column is the member's multi-index ravelled against
the set sizes, offset by the start of the variable's block. `ProductCoord`
performs that computation.

```python
import numpy as np
import nimblend as nb

columns = nb.ProductCoord((2, 3))
print(columns.to_position(np.array([[0, 1], [2, 0]])))
```

Output:

```text
[2 3]
```

Member `(0, 2)` is column 2 and `(1, 0)` is column 3: stride arithmetic and
nothing else.

For a variable over a subset, the column is the member's rank among the
subset's codes. `SubsetCoord` holds the codes in order, so a block already
in canonical order needs no lookup at all.

```python
import numpy as np
import nimblend as nb

columns = nb.SubsetCoord(np.array([0, 4]), (2, 3))
print(columns.to_position(np.array([[0, 1], [0, 1]])))
```

Output:

```text
[0 1]
```

Codes `0` and `4` are members `(0, 0)` and `(1, 1)`, with ranks `0` and `1`.

## Consequences

Nothing stores a column index. A variable over a million members holds its
set sizes, the start of its block and, for a subset, the codes of its
members. It does not hold a million integers recording which column each
member is, because that number is recoverable from the member itself.

The cost of declaring a variable is therefore the cost of its members, not
of its columns. A variable over a full product costs nothing per column:
two set sizes and a start.

A subset variable is not a special case. Both kinds answer the same
question, which position a member occupies, and differ only in whether the
answer is arithmetic or a rank.

---

# /explanation/expressions-are-symbolic

# Expressions are symbolic

A `Term` is a recipe, not a block of numbers. It holds a reference to a
variable, an optional coefficient, the dimensions summed over, a scale
factor, and any lags, conditions or fixed members. An `Expression` is a list
of such terms and the frame they share.

Nothing in that list is an array. `cost[P, W] * x[P, W]` records which
parameter and which variable, and reads neither.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array([f"p{i}" for i in range(200)]))
W = Set("W", np.array([f"w{i}" for i in range(100)]))
cost = Param.from_dense("cost", (P, W), np.ones((200, 100)))

m = Model("transport")
x = m.var("x", (P, W))

expression = Sum(W, cost[P, W] * x[P, W])
print(expression.frame)
print(len(expression.terms))
```

Output:

```text
('P',)
1
```

One term over twenty thousand columns. The same line over twenty million
columns is still one term and costs the same to write.

## Materialisation

An expression becomes matrix entries when it is materialised, which is the
only point at which values are read. Materialisation walks the terms and
issues `nimblend` operations in order: the coefficient is aligned with the
variable's block, the summed dimensions are reduced, the scale is applied,
and the terms are combined over the shared frame.

The result is a `nimblend` array over the frame crossed with the column space:
the block of coefficients the constraint contributes to the matrix.

## A hundred constraints cost a hundred shapes

A constraint holds the recipe rather than the block, so adding one costs
computing its shape. `n_rows` and `nnz` are known at declaration, and no
coefficients exist yet.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array([f"p{i}" for i in range(200)]))
W = Set("W", np.array([f"w{i}" for i in range(100)]))

m = Model("transport")
x = m.var("x", (P, W))
for i in range(20):
    m.eq(f"cap{i}", Sum(W, x[P, W]) <= 1.0)

print(m.n_rows, m.nnz)
```

Output:

```text
4000 400000
```

Four thousand rows and four hundred thousand coefficients are declared, and
the model holds twenty term lists.

## The cost

An expression is materialised twice: once to compute its shape when the
constraint is added, and once to write its entries when the matrix is
assembled. Building twice costs build time. In exchange, one expression is
live at a time rather than all of them, so peak memory is set by the largest
constraint rather than by their sum.

A model that is cheap to declare and more expensive to assemble suits a
builder, because declaration is what a caller iterates on.

---

# /explanation/the-array-is-the-matrix

# The array is the matrix

A constraint is a `nimblend` array indexed over its free sets and the column
space, or over `(ROW, COLUMN)` once its frame has been grouped into rows.
The values of that array are the coefficients. No step converts a model
into a matrix, because the array already is one.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= 1.0)

print(x.terms().dims)
print(m.assemble().to_dense())
```

Output:

```text
('P', 'W', '__column__')
[[1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1.]]
```

The variable's coefficients are indexed over `('P', 'W', '__column__')`.
Grouping the frame into rows puts the same entries over `('__row__',
'__column__')`, which is a matrix in every sense that matters: a row index,
a column index and a value.

## One buffer

A model allocates one `nimblend.EntryBuffer` for its whole matrix. Each
constraint reserves the slice its coefficients need and groups its block
directly into that slice, so the block never exists as a second object.

The frame precedes the column dimension in canonical order, so the grouping
reads a leading prefix and the result is canonical as written: no sort
afterwards, no copy into place.

Handing the matrix to a solver is then a matter of returning views.
`indices` and `values` are the buffer; only `indptr` is built. A model of
four million nonzeros hands over its matrix without copying it.

## Why the shape is computed first

Reserving a slice requires its size, so a constraint computes its shape when
it is added and writes its entries when the model is assembled. The two must
agree, and the model checks that they do.

If a parameter's data changes between the two, the constraint computed one
number of coefficients and built another. The model raises rather than
writing a matrix that does not match the shape it reported.

## What the design gives up

Materialising the expression twice costs build time. That is the price of
holding one expression at a time rather than every constraint's block at
once, and it is a deliberate trade: peak memory is set by the largest
constraint, not by the sum of all of them.

It also means a model's declared shape is exact before anything is built.
`n_rows`, `n_columns` and `nnz` are facts, not estimates.

---

# /explanation/the-package-boundary

# The package boundary

Two packages, with the dependency in one direction. `nimopt` imports
`nimblend`; `nimblend` never imports `nimopt`.

**`nimblend`** is a labeled sparse N-dimensional array. Its vocabulary is
dimensions, labels, entries and alignment. It knows nothing about
optimization, and a function in it referring to a row, a column or a
constraint would be a boundary violation.

**`nimopt`** is an LP/MILP builder in which a variable is a dimension. Its
types are `nimblend` arrays with names attached, and the buffer its matrix
lives in is a `nimblend.EntryBuffer`.

The boundary is not a convention. It is enforced by tests that fail when it
moves.

## What the tests enforce

**`nimblend` never refers to `nimopt`.** A scan of `nimblend`'s source fails on any
mention. A second test imports `nimblend` alone and fails if `nimopt` is
imported with it. `nimblend` carries the other half of the rule in its own
suite: no class, function or parameter it declares is named for a
constraint, an objective, a solver or a variable, and no source file of it
mentions one even in prose. A package developed on its own needs its own
suite to fail, rather than waiting for a consumer's.

**Every `nimblend` import is a public name of the top-level module.** A model
imports `SparseArray` from `nimblend`, never from `nimblend.sparse`. Importing a
public name by its submodule path is how a dependency on an internal
starts. The modules `nimopt` ships are held tighter still: the set of `nimblend`
names they import is pinned, so widening it is a deliberate act rather than
drift.

**No array's `.index` or `.data` is read.** Those are the raw index matrix
and value buffer of the layer below the array layer. The rule is checked by
walking the syntax tree rather than by grepping, so it catches a read that
is not a subscript, and it allows `dims.index(name)`, which is a tuple
being asked for a position.

**And none is assembled.** Reading a buffer is one half of the bypass and
building one is the other. An index matrix assembled in `nimopt` is array
work done a layer too high, and it is exactly the seam that has to move
when the kernel below `nimblend` is replaced. A domain returns the array over
its own members instead: `array(values)` assigns a value to each member,
and `identity(into, coord, start)` pairs each with its position along a new
dimension. No module of `nimopt` calls `SparseArray(index, ...)` or
`from_canonical`.

## Where the bytes land

The sharpest boundary test concerns allocation, and it does not ask which
package holds more.

A build's cost per nonzero lands in `nimblend`: the matrix is an int32 column
and a float64 value per entry, and it lives in a `nimblend` buffer. Increase
the nonzeros ninefold and `nimblend`'s share grows by at least twelve bytes
for each one added.

`nimopt` stays flat across the same change. What it allocates is what a
solver takes beside the matrix: a lower bound, an upper bound and an
integrality flag per column, and a pair of bounds per row. Those are fixed
by the sets, not by the density.

That is why the test is written as a gradient rather than a comparison.
Asserting that `nimblend` simply holds more bytes would measure the workload
instead of the boundary: `nimopt` holds a vector per column whatever the
density, so a sparse enough model puts `nimopt` above `nimblend` with nothing
having drifted.

## Why the site lives in `nimopt`

`nimopt` imports `nimblend`, so a site inside `nimopt` documenting both runs
with the dependency. A site inside `nimblend` documenting `nimopt` would invert
it, and the inversion would be real: `nimblend`'s own tests would need a model
to describe.

---

# /explanation/what-the-numbers-measure

# What the numbers measure

Every figure here is given against a stated baseline. A number without its
denominator says nothing, and a ratio can be made to sound like anything by the choice of denominator.

## Who owns the bytes

Every surviving allocation of a 200x100 transport build, attributed to the package that allocated it. Both rows are the same model over the same sets, so
rows and columns are fixed at 200 and 20 000, and the coefficient's density
alone decides the nonzeros.

| nonzeros | nimblend | nimopt |
|---|---|---|
| 20 000 | 0.490 MB | 0.406 MB |
| 2 209 | 0.204 MB | 0.406 MB |

`nimblend` absorbs 16.0 bytes per added nonzero — an int32 column and a float64
value — because the buffer a model's matrix lives in is a `nimblend.EntryBuffer`
and the CSR handoff returns views of it. `nimopt` moves by 48 bytes across a
ninefold change in the matrix: what it allocates is the per-column vectors,
bounds and integrality, and the row bounds, which the sets fix.

**What it does not claim.** Not that `nimblend` holds more bytes than `nimopt`.
That is a property of a model's nonzeros per column, and a model carrying one
entry per column puts `nimopt` above `nimblend` with nothing having drifted. The
measure is the gradient, not the totals.

## Two models

`benchmarks/bench_transport.py` ships from plants to warehouses over a
network in which each plant serves a band of nearby warehouses, so the flow
variable spans the arcs rather than the full product.

| plants | warehouses | arcs/plant | rows | columns | nonzeros | matrix | peak | ratio | build |
|---|---|---|---|---|---|---|---|---|---|
| 200 | 100 | 10 | 300 | 2 000 | 4 000 | 0.05 MB | 0.39 MB | 8.09x | 8 ms |
| 2 000 | 500 | 20 | 2 500 | 40 000 | 80 000 | 0.96 MB | 6.95 MB | 7.24x | 45 ms |
| 10 000 | 2 000 | 40 | 12 000 | 400 000 | 800 000 | 9.60 MB | 69.71 MB | 7.26x | 314 ms |

`benchmarks/bench_storage.py` dispatches a thermal fleet, a solar fleet and a
set of batteries against an hourly demand. The batteries couple neighbouring
hours through a cyclic state of charge, and the generators through an upward
ramp limit.

| generators | batteries | hours | rows | columns | nonzeros | matrix | peak | ratio | build |
|---|---|---|---|---|---|---|---|---|---|
| 10 | 2 | 168 | 4 862 | 2 688 | 9 724 | 0.12 MB | 0.95 MB | 8.10x | 26 ms |
| 40 | 8 | 720 | 81 320 | 46 080 | 166 960 | 2.00 MB | 14.76 MB | 7.37x | 69 ms |
| 80 | 20 | 8 760 | 2 111 080 | 1 226 400 | 4 379 840 | 52.56 MB | 370.45 MB | 7.05x | 1 552 ms |

Both solve through HiGHS. The two lag rules are visible in the row counts: over 168 hours the ramp block has 10 x 167 rows, because the first hour has no predecessor and its row is not produced, while the cyclic state of charge wraps to the last hour and keeps all 2 x 168 of its rows.

## Peak against the matrix

The peak is reached while the model is built, not while it is assembled. On
the 400 000-column transport model, building peaks at 69.72 MB and the
assembly that follows peaks at 68.23 MB, because every constraint's
expression is built once to measure its shape and the model that produces the
matrix stays live while the matrix is written.

The denominator decides how large the ratio sounds. The CSR matrix is
9.60 MB; the data a solver takes — matrix, row pointer, row and column
bounds, cost, integrality — is 21.04 MB, and the model that produced it is
21.74 MB more. Peak is 7.26x the matrix and 3.31x the full LP data, and
45.99 MB of the 69.72 MB is still live when the build returns.

The ratio is close to flat across the rungs, 8.09x at 4 000 nonzeros and
7.26x at 800 000, and the smallest rung's figure is stable rather than a
first-call artefact: three consecutive measurements of it in one process give
0.39, 0.38 and 0.38 MB. What the small rung has above the large is the
model's own fixed structures, which do not shrink with the matrix.

Where the transient bytes go, attributed by the frame that allocated them at
the moment the 400 000-column build peaks:

| bytes | allocated by |
|---|---|
| 16.26 MB | the caller's own arc labels and costs |
| 14.40 MB | the broadcast product of a parameter against a variable |
| 12.80 MB | the sorted copies the two unordered blocks need |
| 8.00 MB | the variable's own index, values and column coordinate |
| 6.40 MB | the subset's codes and its index in code order |
| 3.20 MB | the positions a probe resolves to |
| 3.20 MB | a ravelled key set |

An `int64` ravel key costs 8 bytes per entry and an argsort permutation
another 8, against the 12 bytes a matrix entry finally occupies, so an
operation that aligns costs more than the result it produces. That is the
floor the ratio rests on, and it is per-operation and transient rather than
retained.

A second constraint costs its own rows and little else. Over a 500 000-cell
model, one constraint peaks at 65.71 MB against an 8.00 MB buffer and two
peak at 75.85 MB against a 16.00 MB buffer: 2.14 MB beyond the rows the
second adds, because a constraint holds its term list rather than a block, so
an expression exists while its shape is measured and again while it is
written, and never between.

**What it does not claim.** Not that peak memory is 7x the matrix in any
absolute sense. It is 7.26x *the CSR matrix* and 3.31x *the full LP data* on
this model, and which of those a reader cares about depends on what they were
going to compare it with.

## Against linopy

`benchmarks/bench_vs_linopy.py` builds the same three models both ways.
Inputs — the arc list, the hourly profiles, the costs — are prepared outside
the measured region and handed to both. The measured build runs from an empty
model to the matrix a solver would be given: `assemble()` for nimopt,
`m.matrices` for linopy. Each side runs in a process of its own, because
resident memory carries the import cost of whichever library is loaded, and a
row is printed only once the two agree on rows, columns, nonzeros and the
solved objective.

Resident memory is sampled rather than traced. `tracemalloc` sees only what
passes through Python's allocator, and two libraries that allocate memory by different routes would be compared on the route rather than on the memory.

**A variable over a sparse subset.** The flow spans the arcs in nimopt and the
full plant-by-warehouse product under a mask in linopy.

| arcs | matrix | nimopt RSS | linopy RSS | nimopt build | linopy build |
|---|---|---|---|---|---|
| 2 000 | 0.05 MB | 3.4 MB | 27.3 MB | 3.8 ms | 177.7 ms |
| 40 000 | 0.96 MB | 10.7 MB | 107.5 MB | 27.9 ms | 206.5 ms |
| 400 000 | 9.60 MB | 72.9 MB | 1 682.9 MB | 292.9 ms | 844.5 ms |

At 400 000 arcs the mask spans a 20 000 000-cell product: nimopt builds the
same matrix in a twenty-third of the memory and a third of the time. This is
the axis the design is for.

**Temporal coupling on a dense model.** Ramp limits and a cyclic state of
charge over a full generator-by-hour grid.

| rows | matrix | nimopt RSS | linopy RSS | nimopt build | linopy build |
|---|---|---|---|---|---|
| 4 862 | 0.12 MB | 4.3 MB | 28.8 MB | 8.5 ms | 267.5 ms |
| 81 320 | 2.00 MB | 18.9 MB | 41.4 MB | 58.6 ms | 280.4 ms |
| 2 111 080 | 52.56 MB | 383.4 MB | 409.3 MB | 1 511.0 ms | 507.9 ms |

Here the advantage narrows and then reverses on time: at 2 111 080 rows
linopy builds the same matrix **3.0x faster** for 7% more resident memory.
Nothing is sparse in this model, so what remains is what the design costs where it gains nothing — every operation aligns by key, while xarray
broadcasts over dense grids and aligns by position — and the model is built
twice, once to measure each constraint's shape and once to write it. Resident
memory stays close because both end up holding the same dense coefficient
grids.

**Integrality.** Making the flow an integer column moves neither side's
build: nimopt 8.6 MB and 27.3 ms against its own 10.7 MB and 27.9 ms as an LP,
linopy 108.4 MB and 204.4 ms against 107.5 MB and 206.5 ms. Integrality is a
column vector, not a matrix.

**Reading the answer back.** Primals onto their sets and duals onto their
rows: nimopt 0.7–6.9 ms across every rung, linopy 1.4–64.5 ms, the 64.5 ms
being the 400 000-arc transport model, where the answer is unpacked onto the same masked product the build used.

The matrices differ in width: linopy hands over a scipy matrix with `int64`
column indices, so the same 800 000 nonzeros occupy 12.80 MB against nimopt's
9.60 MB.

**What it does not claim.** Not that nimopt is faster than linopy. On a dense
temporally-coupled model at two million rows it is three times slower, and
that row is in the table for the same reason the others are. The claim the
numbers support is narrower: where a model is sparse in its variables, not
materialising the grid is worth a great deal, and where it is not, the alignment work is a cost with no corresponding saving.

---

# /tutorial/sets-and-parameters

# Sets and parameters

The tutorial builds one model over six pages, the transport problem from
[Get started](/get-started), one concept per page. This page declares the
index sets and the data.

## The problem

Two plants, Lisbon and Porto, ship to three warehouses, Berlin, Paris and
Rome. Plant `p` has supply `s[p]`, warehouse `w` has demand `d[w]`, and one
unit shipped on route `(p, w)` costs `c[p, w]`. The decision is the quantity
`x[p, w]` shipped on each of the six routes, and the objective is total
cost.

| | Berlin | Paris | Rome | Supply |
| --- | --- | --- | --- | --- |
| Lisbon | 2 | 4 | 5 | 30 |
| Porto | 3 | 1 | 6 | 25 |
| Demand | 20 | 15 | 15 | |

## Sets

A `Set` is a named index dimension with labels. Parameters, variables and
constraints are indexed over sets, and solution values are returned over the
same sets.

```python
import numpy as np
from nimopt import Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

print(P.labels)
print(len(W))
print(P.position_of(np.array(["porto"])))
```

Output:

```text
['lisbon' 'porto']
3
[1]
```

`P` has two members and `W` three. `position_of` maps labels to their
integer positions, which are the indices used internally.

## Parameters

A `Param` is data indexed over a set product: one value per combination of
members. `Param.from_dense` takes an array whose shape equals the sizes of
the sets, in order. Cost is indexed over `(P, W)`, supply over `P`, and
demand over `W`.

```python
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

print(cost.dims, cost.nnz)
print(cost.materialise().to_dense())
print(supply.dims, demand.dims)
```

Output:

```text
('P', 'W') 6
[[2. 4. 5.]
 [3. 1. 6.]]
('P',) ('W',)
```

`cost` has six entries. `materialise()` returns the parameter as a `nimblend`
array, and `to_dense()` renders it as a NumPy array with axes in the
declared set order.

A shape mismatch raises `ValueError`; the message gives the expected and the
actual shape.

```python raises=ValueError
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

Param.from_dense("cost", (P, W), np.array([[2.0, 4.0], [3.0, 1.0]]))
```

Raises ValueError:

```text
ValueError: parameter 'cost' is over sets of shape (2, 3); got values of shape (2, 2)
```

## Sparse data

In a sparse network, a plant serves a subset of the warehouses, and the cost
parameter has entries only on existing routes. `Param.from_long` takes the
entries in long form: one label column per set and one value column, read in
parallel, so the k-th entry of each column belongs to the same route.

```python
import numpy as np
from nimopt import Param, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

cost = Param.from_long(
    "cost",
    (P, W),
    {
        "P": np.array(["lisbon", "lisbon", "porto"]),
        "W": np.array(["berlin", "rome", "paris"]),
    },
    np.array([2.0, 5.0, 1.0]),
)

print(cost.nnz)
print(cost.materialise().to_dense())
```

Output:

```text
3
[[2. 0. 5.]
 [0. 1. 0.]]
```

Three routes, three entries. `to_dense()` prints zeros at the three missing
routes, but the parameter stores nothing there: an unlisted route is absent,
not zero. The distinction matters on the last page of the tutorial, where a
variable declared over exactly these routes has no column for the others.

Next: [Variables](/tutorial/variables).

---

# /tutorial/variables

# Variables

The decision is the quantity shipped on each route: one variable indexed
over plants and warehouses.

## Declaring a variable

A `Model` holds variables, constraints and the objective. `m.var(name,
sets)` declares a variable indexed over a tuple of sets and returns a handle
for use in expressions.

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

print(x.dims)
print(x.n_columns)
```

Output:

```text
('P', 'W')
6
```

Two plants by three warehouses gives six members, so `x` occupies six
columns of the coefficient matrix. A column index is computed from a
member's positions in each set; nothing stores a column per member, so a
variable over a million members costs the same to declare as one over six.

## Bounds and integrality

A variable has a lower bound of 0 and no upper bound unless declared
otherwise. `lower=` and `upper=` take a number that applies to every column.
`integer=True` restricts the columns to integer values, which makes the
model a MILP.

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W), upper=20.0)
trucks = m.var("trucks", (P,), integer=True)

lower, upper = m.column_bounds()
print(lower)
print(upper)
print(m.integrality())
print(m.n_columns)
```

Output:

```text
[0. 0. 0. 0. 0. 0. 0. 0.]
[20. 20. 20. 20. 20. 20. inf inf]
[0 0 0 0 0 0 1 1]
8
```

A model has one column space shared by all its variables: `x` occupies
columns 0 to 5 and `trucks` columns 6 and 7. `column_bounds()` returns the
lower and upper bound vectors in column order, and `integrality()` returns
one flag per column.

A parameter in place of a number gives each column its own bound; see
[Bounds from a parameter](/guides/bounds-from-parameters).

Next: [Expressions](/tutorial/expressions).

---

# /tutorial/expressions

# Expressions

Constraints and the objective are stated over sums of variables: the total
shipped from a plant, the total received by a warehouse, the total cost. An
expression is such a sum. It is symbolic: writing one records the variables,
coefficients and sets involved, and computes nothing.

## Referencing a variable

`x[P, W]` references the variable over its sets and returns an expression
with one term. The **frame** of an expression is the tuple of dimensions it
is still indexed over. `x[P, W]` has frame `(P, W)`: one value per route.

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

shipped = x[P, W]
print(type(shipped).__name__)
print(shipped.frame)
```

Output:

```text
Expression
('P', 'W')
```

## Sum

`Sum(S, expression)` sums over the members of `S` and removes `S` from the
frame. Summing over `W` gives the total shipped from each plant, indexed
over `P`. Summing over both sets gives a scalar.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

print(x[P, W].frame)
print(Sum(W, x[P, W]).frame)
print(Sum(P, W, x[P, W]).frame)
```

Output:

```text
('P', 'W')
('P',)
()
```

The frame determines the shape of a constraint built on the expression: an
expression with frame `(P,)` produces one row per plant. An expression with
an empty frame is a scalar, which is the form an objective takes.

## Coefficients

The cost of a plan is `Σ_{p,w} c[p, w] · x[p, w]`. Multiplying a reference
by a parameter over the same sets gives the term a coefficient. The frame is
unchanged until the sum is taken.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))

m = Model("transport")
x = m.var("x", (P, W))

per_route = cost[P, W] * x[P, W]
total_cost = Sum(P, W, per_route)
print(per_route.frame)
print(total_cost.frame)
print(len(total_cost.terms))
```

Output:

```text
('P', 'W')
()
1
```

`total_cost` is a single term. The same expression over a million routes is
still one term, because it holds references to `cost` and `x` rather than
their values. Values are read when the matrix is assembled.

## Addition and subtraction

Expressions add and subtract, producing one expression over the frame both
share. A balance, inflow minus outflow, is written this way. With a second
variable for returned goods, the net shipment on a route is the outbound
quantity minus the returned quantity.

```python
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))
returned = m.var("returned", (P, W))

net = x[P, W] - returned[P, W]
print(net.frame)
print(len(net.terms))
```

Output:

```text
('P', 'W')
2
```

Two terms, one per variable, over the same frame.

Next: [Constraints](/tutorial/constraints).

---

# /tutorial/constraints

# Constraints

The model has two constraint families: a supply limit per plant and a
demand requirement per warehouse.

```text
Σ_w x[p,w] ≤ s[p]        for each plant p
Σ_p x[p,w] ≥ d[w]        for each warehouse w
```

Each family is one line of code and produces one row per member of its
frame.

## Relations

Comparing an expression with `<=`, `>=` or `==` produces a `Relation`: the
expression, the sense, and the right-hand side. A relation is not yet part
of the model.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))

m = Model("transport")
x = m.var("x", (P, W))

rule = Sum(W, x[P, W]) <= supply[P]
print(type(rule).__name__, rule.sense)
```

Output:

```text
Relation <=
```

## Adding a constraint

`m.eq(name, relation)` adds the relation to the model under a name and
returns the `Constraint`. The name identifies the constraint's rows in the
matrix and its dual values in the solution. A constraint produces one row
per member of its expression's frame.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))

supply_rows = m.eq("supply", Sum(W, x[P, W]) <= supply[P])
demand_rows = m.eq("demand", Sum(P, x[P, W]) >= demand[W])

print(supply_rows.n_rows, supply_rows.nnz)
print(demand_rows.n_rows, demand_rows.nnz)
print(m.n_rows, m.nnz)
```

Output:

```text
2 6
3 6
5 12
```

The supply expression has frame `(P,)` and produces two rows; the demand
expression has frame `(W,)` and produces three. Each supply row has three
nonzeros, one per route out of its plant, and each demand row two, one per
route into its warehouse: twelve nonzeros in total.

## The right-hand side

The right-hand side is a scalar, applied to every row, or a parameter read
at exactly the constraint's frame, giving each row its own value. The
parameter is read at its sets here as it is anywhere else: `supply[P]`, not
`supply`. A bare name raises `TypeError` and names the reading it wants,
because a row bounded by a bare name reads as one bounded by a scalar.

Supply is indexed over `P`, as are the supply rows. A parameter read over any
other index set raises `ValueError`; the message gives both index sets.

```python raises=TypeError
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))

m = Model("transport")
x = m.var("x", (P, W))

m.eq("supply", Sum(W, x[P, W]) <= supply)
```

Raises TypeError:

```text
TypeError: parameter 'supply' carries ('P',) and states no coefficient until it is read; read it at its sets as supply[P]
```

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))

m.eq("supply", Sum(W, x[P, W]) <= demand[W])
```

Raises ValueError:

```text
ValueError: constraint 'supply' has free dimensions ('P',); its right-hand side 'demand' is over ('W',)
```

## One bound per constraint

Python evaluates the chained comparison `0 <= expr <= 10` as
`(0 <= expr) and (expr <= 10)`, which discards the first relation. `nimopt`
raises `TypeError` on the chained form rather than dropping a bound. Each
bound is written as its own constraint.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))

m = Model("transport")
x = m.var("x", (P, W))

0.0 <= Sum(W, x[P, W]) <= 10.0
```

Raises TypeError:

```text
TypeError: a relation has no truth value; a chained comparison such as 0 <= expr <= 10 reads as two comparisons joined by `and` and keeps only the second, so state each bound separately
```

Next: [Solving](/tutorial/solving).

---

# /tutorial/solving

# Solving

The remaining pieces are the objective function and the solver call.

## The objective

`m.set_objective(expression)` takes an expression with an empty frame. The
model's `sense`, `"min"` by default, sets the direction. `m.solve()`
assembles the matrix, calls HiGHS and returns a `Solution`.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

solution = m.solve()
print(solution.status)
print(solution.objective)
```

Output:

```text
optimal
135.0
```

The status is `optimal` and the objective value is 135.

## Status

`status` reports the outcome of the solve and is always readable.
`objective`, `primal` and `dual` are defined only for an optimal solve;
reading one after any other status raises `ValueError`.

Raising Berlin's demand to 40 makes total demand 70 against total supply
55, so the model is infeasible.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([40.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

solution = m.solve()
print(solution.status)
solution.objective
```

Raises ValueError:

```text
infeasible
ValueError: the model's status is 'infeasible', so it carries no objective; read `status` before reading values
```

For an infeasible model whose cause is not evident, `m.session()` keeps the
solver instance open and `diagnose()` returns the conflicting rows. See
[Solvers](/reference/solvers).

## The matrix

`m.assemble()` builds the coefficient matrix without a solver call and
returns it in CSR form. `to_dense()` renders it for a model of this size,
and `row_of(name)` gives the row range of a named constraint.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])

assembled = m.assemble()
print(m.n_rows, m.n_columns, m.nnz)
print(assembled.row_of("demand"))
print(assembled.to_dense())
```

Output:

```text
5 6 12
slice(2, 5, None)
[[1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1.]
 [1. 0. 0. 1. 0. 0.]
 [0. 1. 0. 0. 1. 0.]
 [0. 0. 1. 0. 0. 1.]]
```

The six columns are the routes, Lisbon's three followed by Porto's. Rows 0
and 1 are the supply rows, each with a 1 under its plant's three routes.
Rows 2 to 4 are the demand rows, each with a 1 under the two routes into its
warehouse. `row_of("demand")` returns that range.

Next: [Reading the answer](/tutorial/reading-the-answer).

---

# /tutorial/reading-the-answer

# Reading the answer

A solver returns primal and dual values as flat vectors. `nimopt` returns
them as arrays over the sets each variable and constraint was declared on.

## Primals and duals

`solution.primal(name)` returns a variable's values over its sets.
`solution.dual(name)` returns a constraint's dual values over its frame: for
the demand constraint, one value per warehouse.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

shipped = solution.primal("x")
print(shipped.dims)
print(shipped.to_dense())
print(solution.dual("demand").dims)
print(solution.dual("demand").to_dense())
```

Output:

```text
('P', 'W')
[[20.  0. 10.]
 [ 0. 15.  5.]]
('W',)
[3. 1. 6.]
```

Rows are plants and columns are warehouses: Lisbon ships 20 to Berlin and
10 to Rome, Porto ships 15 to Paris and 5 to Rome. The dual of a demand row
is the increase in total cost per additional unit of demand at that
warehouse: 3 in Berlin, 1 in Paris, 6 in Rome. Lisbon's supply constraint
binds, so each marginal unit is served from Porto at the cost of Porto's
route.

## A variable over a subset

If Porto cannot ship to Rome, the variable is declared over the five
existing routes with `subset=`. The sixth route has no column, and the
solution has values at the five members only.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum, subset

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))

routes = subset(
    (P, W),
    {
        "P": np.array(["lisbon", "lisbon", "lisbon", "porto", "porto"]),
        "W": np.array(["berlin", "paris", "rome", "berlin", "paris"]),
    },
)

m = Model("transport")
x = m.var("x", (P, W), subset=routes)
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

print(x.n_columns)
print(solution.objective)
print(type(solution.primal("x")).__name__)
```

Output:

```text
5
135.0
SparseArray
```

Five columns instead of six. The objective is unchanged at 135, because
Rome is served from Lisbon in both solutions. The array type differs: a
variable over a full product returns a `DenseArray`, a variable over a
subset a `SparseArray` with an entry per member and nothing elsewhere.

## Absence is not zero

The model contains no decision for the route Porto to Rome. Every array a
solution returns declares `absence="unknown"`, and `to_dense()` raises
`ValueError` rather than filling the missing coordinate with a value the
model never produced.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set, Sum, subset

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))
routes = subset(
    (P, W),
    {
        "P": np.array(["lisbon", "lisbon", "lisbon", "porto", "porto"]),
        "W": np.array(["berlin", "paris", "rome", "berlin", "paris"]),
    },
)

m = Model("transport")
x = m.var("x", (P, W), subset=routes)
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

solution.primal("x").to_dense()
```

Raises ValueError:

```text
ValueError: this array declares absence 'unknown' and does not carry every coordinate of its frame, so densifying must state fill=<value> to place at the rest
```

`to_dense(fill=...)` supplies the value for missing coordinates. `nan`
distinguishes a missing route from a route with zero shipment.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum, subset

P = Set("P", np.array(["lisbon", "porto"]))
W = Set("W", np.array(["berlin", "paris", "rome"]))
cost = Param.from_dense("cost", (P, W), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))
supply = Param.from_dense("supply", (P,), np.array([30.0, 25.0]))
demand = Param.from_dense("demand", (W,), np.array([20.0, 15.0, 15.0]))
routes = subset(
    (P, W),
    {
        "P": np.array(["lisbon", "lisbon", "lisbon", "porto", "porto"]),
        "W": np.array(["berlin", "paris", "rome", "berlin", "paris"]),
    },
)

m = Model("transport")
x = m.var("x", (P, W), subset=routes)
m.eq("supply", Sum(W, x[P, W]) <= supply[P])
m.eq("demand", Sum(P, x[P, W]) >= demand[W])
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))
solution = m.solve()

values = solution.primal("x")
print(values.absence, values.nnz)
print(values.to_dense(fill=np.nan))
```

Output:

```text
unknown 5
[[15.  0. 15.]
 [ 5. 15. nan]]
```

Lisbon serves Rome alone, and the Porto to Rome cell reads `nan`. A stored
`0.0` would denote a route that exists and carries nothing.

## Summary

The tutorial declared index sets and parameters, a decision variable,
expressions, two constraint families and an objective, solved the model,
and read the solution back over its sets. The [guides](/guides/subsets)
cover variables over sparse networks, conditions on rows, lags, bounds from
data, and models with millions of rows. The [worked models](/models) show
eight complete formulations.

---

# /guides/at-scale

# Stating a model at scale

Models with millions of columns are declared the same way as small ones.
Declaring costs shapes, not blocks: a constraint computes its row and
nonzero counts when it is added, and no matrix exists until `assemble()` or
`solve()` is called.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array([f"p{i}" for i in range(200)]))
W = Set("W", np.array([f"w{i}" for i in range(100)]))
cost = Param.from_dense("cost", (P, W), np.ones((200, 100)))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= 1.0)
m.eq("demand", Sum(P, x[P, W]) >= 1.0)
m.set_objective(Sum(P, W, cost[P, W] * x[P, W]))

print(m.n_columns)
print(m.n_rows)
print(m.nnz)
```

Output:

```text
20000
300
40000
```

Twenty thousand columns, three hundred rows and forty thousand nonzeros are
declared, and the model holds no matrix. The parameter is the caller's own
data; the model has added two constraints and an objective, each a symbolic
expression recording which variable, which coefficient and which sets are
summed.

## Building the matrix once

`assemble()` builds the matrix into one buffer and returns it in CSR form.
Each constraint writes its rows into its own slice, so the matrix exists
once and one expression at a time is materialised.

```python
import numpy as np
from nimopt import Model, Set, Sum

P = Set("P", np.array([f"p{i}" for i in range(200)]))
W = Set("W", np.array([f"w{i}" for i in range(100)]))

m = Model("transport")
x = m.var("x", (P, W))
m.eq("supply", Sum(W, x[P, W]) <= 1.0)
m.eq("demand", Sum(P, x[P, W]) >= 1.0)

assembled = m.assemble()
print(assembled.n_rows, assembled.n_cols)
print(assembled.values.shape)
print(assembled.row_of("demand"))
```

Output:

```text
300 20000
(40000,)
slice(200, 300, None)
```

`indices` and `values` are views of that buffer; only `indptr` is built.
`row_of(name)` gives a constraint's rows as a slice, which is how a named
row block is located in a matrix too large to print.

## Reading a large solution

`to_dense()` is for a model small enough to print. At scale, read the CSR
arrays, or read the solution through `primal()` and `dual()`, which return
labeled arrays over the sets rather than offsets into a vector.

A variable over a subset keeps its values sparse, so reading a solution over
a sparse network does not build the grid the model avoided.

## What a variable costs

A member's column is computed from its multi-index rather than stored, so a
variable over millions of columns costs its members, not its columns.
`n_columns` above is twenty thousand while the variable holds a few
numbers: the start of its block and the sizes of its sets.

## Progress

A model of size takes seconds to build. `progress=True` reports the build
in a terminal.

```python skip="the report draws to a terminal, and a page is not one"
from nimopt.models import storage

model = storage.definition().build(storage.data(60), progress=True)
solution = model.solve(options={"time_limit": 300.0, "log": True}, progress=True)
```

---

# /guides/bounds-from-parameters

# Bounds from a parameter

A capacity per generator, an energy limit per battery, a flow limit per
line: bounds usually come from data and differ by member. `lower=` and
`upper=` accept a number, applied to every column, or a `Param`, giving
each member its own value.

```python
import numpy as np
from nimopt import Model, Param, Set

G = Set("G", np.array(["g0", "g1"]))
cap = Param.from_dense("cap", (G,), np.array([5.0, 7.0]))

m = Model("schedule")
m.var("x", (G,), upper=cap)

lower, upper = m.column_bounds()
print(lower)
print(upper)
```

Output:

```text
[0. 0.]
[5. 7.]
```

## Broadcasting a narrower parameter

A parameter indexed over fewer dimensions than the variable is broadcast
over the rest. A capacity per unit bounds every period of that unit.

```python
import numpy as np
from nimopt import Model, Param, Set

G = Set("G", np.array(["g0", "g1"]))
T = Set("T", np.array(["t0", "t1", "t2"]))
cap = Param.from_dense("cap", (G,), np.array([5.0, 7.0]))

m = Model("schedule")
m.var("x", (G, T), upper=cap)

print(m.column_bounds()[1])
```

Output:

```text
[5. 5. 5. 7. 7. 7.]
```

Six columns; each unit's three periods take that unit's capacity.

## A bound covers every column

A bound with no value for some member of the variable raises `ValueError`;
the message gives the member. A dense parameter covers its product by
construction; a long-form parameter can omit a member.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set

S = Set("S", np.array(["a", "b", "c"]))
cap = Param.from_long("cap", (S,), {"S": np.array(["a", "c"])}, np.array([1.0, 3.0]))

m = Model("bounds")
m.var("x", (S,), upper=cap)

m.column_bounds()
```

Raises ValueError:

```text
ValueError: the upper bound 'cap' carries no value for member ('b',) of variable 'x'; a bound covers every column of the variable it bounds
```

The check runs when the bound vectors are built. `column_bounds()`,
`assemble()` and `solve()` all build them.

## A bound over a dimension the variable lacks

A parameter indexed over a dimension the variable is not declared on raises
`ValueError`. A bound is per column, and a dimension the variable does not
have selects no column.

```python raises=ValueError
import numpy as np
from nimopt import Model, Param, Set

G = Set("G", np.array(["g0", "g1"]))
W = Set("W", np.array(["w0", "w1"]))
cap = Param.from_dense("cap", (W,), np.array([5.0, 7.0]))

m = Model("schedule")
m.var("x", (G,), upper=cap)
```

Raises ValueError:

```text
ValueError: variable 'x' is declared over ('G',) and does not carry ['W']; its upper bound 'cap' is declared over ('W',)
```

---

# /guides/coefficient-arithmetic

# Coefficient arithmetic

A coefficient is often derived from several parameters: fuel price divided
by efficiency, a cost scaled by a factor. In `nimopt` a coefficient is a
parameter read at its sets or an arithmetic combination of such readings.
`+`, `-`, `*`, `/` and a power by a number combine them. The combination is
symbolic: it holds references, derives its dimensions from its operands, and
is evaluated once, when the term it multiplies is materialised. A derived
coefficient can therefore appear in a definition before any data exists.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

G = Set("G", np.array(["base", "peak"]))
T = Set("T", np.arange(3))
price = Param.from_dense("fuel_price", (G, T), np.full((2, 3), 30.0))
eta = Param.from_dense("efficiency", (G, T), np.array([[0.5] * 3, [0.4] * 3]))
load = Param.from_dense("load", (T,), np.full(3, 100.0))
cap = Param.from_dense("capacity", (G, T), np.full((2, 3), 80.0))

unit_cost = price[G, T] / eta[G, T]

m = Model("dispatch", sense="min")
gen = m.var("gen", (G, T), lower=0.0, upper=cap)
m.eq("balance", Sum(G, gen[G, T]) == load[T])
m.set_objective(Sum(G, T, unit_cost[G, T] * gen[G, T]))

print(unit_cost.name, unit_cost.dims)
print(m.solve().objective)
```

Output:

```text
(fuel_price / efficiency) ('G', 'T')
18900.0
```

## Reading a derived coefficient

`unit_cost[G, T]` reads a combination the same way `price[G, T]` reads a
parameter, and the sets given are checked against the combination's
dimensions. A transposed or incomplete index raises `ValueError`.

```python raises=ValueError
import numpy as np
from nimopt import Param, Set

G = Set("G", np.array(["base", "peak"]))
T = Set("T", np.arange(3))
price = Param.from_dense("fuel_price", (G, T), np.full((2, 3), 30.0))
eta = Param.from_dense("efficiency", (G, T), np.array([[0.5] * 3, [0.4] * 3]))

(price[G, T] / eta[G, T])[T, G]
```

Raises ValueError:

```text
ValueError: coefficient (fuel_price / efficiency) is over ('G', 'T'); got ('T', 'G')
```

A bare parameter has no arithmetic: `price * 2.0` raises `TypeError`. Read
the parameter at its sets first and combine the references.

## Alignment

Two operands with the same dimensions align entry by entry. Operands whose
dimensions nest or overlap align on the shared dimensions and broadcast over
the rest, with the left operand's order first. Operands sharing no dimension
raise `ValueError`: their product would be an outer product, which a linear
model does not require. The same rule applies when a coefficient multiplies
a variable. A number has no dimensions and scales.

```python raises=ValueError
import numpy as np
from nimopt import Param, Set

G = Set("G", np.array(["base", "peak"]))
T = Set("T", np.arange(3))
over_g = Param.from_dense("over_g", (G,), np.ones(2))
over_t = Param.from_dense("over_t", (T,), np.ones(3))

over_g[G] * over_t[T]
```

Raises ValueError:

```text
ValueError: frames ('G',) and ('T',) share no dimension; there is nothing to align them on
```

**A coefficient never adds a column the variable does not have.** The
variable's members define the columns; a coefficient can only reduce which
of them receive a nonzero. A coefficient with dimensions the variable lacks
defines rows over those dimensions: this is how a term maps rows to columns.

## Division by zero

A divisor that is zero raises `ZeroDivisionError`; the message gives the
coordinate. Substituting infinity would hand the solver a model nobody
wrote. The check covers a Python number, a NumPy scalar and a coefficient
with a zero at any coordinate, because a NumPy scalar divides to infinity
where a Python number raises.

```python raises=ZeroDivisionError
import numpy as np
from nimopt import Param, Set

G = Set("G", np.array(["base", "peak"]))
T = Set("T", np.arange(3))
price = Param.from_dense("fuel_price", (G, T), np.full((2, 3), 30.0))
holed = Param.from_dense("holed", (G, T), np.array([[0.5, 0.0, 0.5], [0.4] * 3]))

price[G, T] / holed[G, T]
```

Raises ZeroDivisionError:

```text
ZeroDivisionError: divisor holed carries a zero at 1 coordinate(s), the first at {'G': 'base', 'T': 1}; a quotient there states a coefficient no solver can read
```

## Constants

An expression consists of terms and a constant. `x + 1 <= 5` produces the
row `x <= 4`: the constant moves to the right-hand side when the constraint
is built, and into the reported objective value as a fixed cost. A constant
adds no column. A constant alone is neither an objective nor a constraint,
and raises.

```python
import numpy as np
from nimopt import Model, Param, Set, Sum

P = Set("P", np.array(["a"]))
one = Param.from_dense("one", (P,), np.ones(1))

m = Model("m", sense="max")
x = m.var("x", (P,), upper=100.0)
m.eq("cap", one[P] * x[P] + 1.0 <= 5.0)
m.set_objective(Sum(P, one[P] * x[P]) + 7.0)

print(m.assemble().n_cols, m.assemble().row_upper)
print(m.solve().objective)
```

Output:

```text
1 [4.]
11.0
```

## Refused forms

A line of modelling arithmetic produces a linear term, or raises with a
message that gives the form to write instead. There is no third outcome.

| Written | The message says to write |
| --- | --- |
| `x[P] * y[P]` | a linear term has one variable; a coefficient multiplies it |
| `x[P] ** 2` | a coefficient takes the power, and a variable multiplies it |
| `x[P] / y[P]` | a variable in a denominator is not linear |
| `2.0 / x[P]` | the same: write the reciprocal as a coefficient |
| `x[P] / 0.0` | a divisor of zero is handled before it reaches an expression |
| `x[P] < 1.0` | `<=` and `>=`; an LP has no row for a strict inequality |
| `x[P] > 1.0` | the same |
| `x[P] != 1.0` | one bound per equation |
| `0.0 <= x[P] <= 1.0` | each bound as its own equation |
| `x[P] + c[P]` | a coefficient has no row until a variable multiplies it |
| `abs(x[P])`, `min` and `max` | reduce with `Sum`, or bound the expression with two rows |
| `np.sum(x[P])` | `Sum` and its sets |
| `np.array([...]) * x[P]` | `Param.from_dense`, read at its sets |
| `Sum(P, Sum(P, x[P]))` | a dimension is reduced once |
| `Sum(P - 1, x[P])` | `Sum(P, x[P - 1])`: the lag belongs on the reference |
| `T - 1.7` | a lag is a whole number of members |

Each raises a `nimopt` message rather than a bare Python error, and the test
suite executes both the table and the rewrites the messages name.

```python raises=TypeError
import numpy as np
from nimopt import Model, Set

P = Set("P", np.array(["a", "b"]))
m = Model("m")
x = m.var("x", (P,))

np.array([1.0, 2.0]) * x[P]
```

Raises TypeError:

```text
TypeError: a coefficient is a parameter; build one with `Param.from_dense` or `Param.from_long` and read it at its sets. A product of two expressions is not linear.
```

A NumPy array multiplying a term would let NumPy handle the operator and
return an array of expressions. The expression types refuse the ufunc, and
the message says how a coefficient is built.

---

# /guides/conditions

# Conditions on a sum and on a constraint

Two situations call for a condition. A constraint may sum over part of a
variable's members, such as the arcs of a network when the variable is
indexed over the full product. And a constraint may apply to some members
of its frame only, such as a capacity limit on one plant. `where=` handles
both. `over=` states a constraint's rows explicitly.

## Restricting a sum

`Sum(..., where=domain)` restricts each term to the members of `domain`
before summing.

```python
import numpy as np
from nimopt import Model, Set, Sum, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

arcs = subset(
    (P, W),
    {"P": np.array(["p0", "p0", "p1"]), "W": np.array(["w0", "w1", "w2"])},
)

m = Model("network")
x = m.var("x", (P, W))
rows = m.eq("capacity", Sum(W, x[P, W], where=arcs) <= 10.0)

print(x.n_columns)
print(rows.n_rows, rows.nnz)
print(m.assemble().to_dense())
```

Output:

```text
6
2 3
[[1. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1.]]
```

The variable is over the full product and has six columns; the condition
leaves three of them with a coefficient. A subset variable would have three
columns in the first place. Use a condition when the variable is over the
product and one constraint reads part of it; declare a subset when the
model never uses the other members.

## Restricting the rows

`m.eq(..., where=domain)` takes a domain over the constraint's frame and
keeps the rows in it. A row outside the condition is not produced.

```python
import numpy as np
from nimopt import Model, Set, Sum, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("network")
x = m.var("x", (P, W))
only_p0 = subset((P,), {"P": np.array(["p0"])})
rows = m.eq("capacity", Sum(W, x[P, W]) <= 10.0, where=only_p0)

print(rows.n_rows)
print(m.assemble().to_dense())
```

Output:

```text
1
[[1. 1. 1. 0. 0. 0.]]
```

One row, for `p0`. `p1` has no capacity row.

A condition over dimensions other than the constraint's frame raises
`ValueError`; the message gives both index sets.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("network")
x = m.var("x", (P, W))
by_warehouse = subset((W,), {"W": np.array(["w0"])})

m.eq("capacity", Sum(W, x[P, W]) <= 10.0, where=by_warehouse)
```

Raises ValueError:

```text
ValueError: constraint 'capacity' has free dimensions ('P',); its condition is over ('W',)
```

## Stating the rows explicitly

By default the rows of a constraint are derived from its terms: a row
exists where every term has a value and the right-hand side has a value. A
term with no value along a frame dimension removes the row, because a row
missing one of its terms would express a constraint that was not written.

`over=domain` states the rows instead of deriving them. A term with values
at some of the rows contributes where it has them, and every row in the
domain is produced.

```python
import numpy as np
from nimopt import Model, Set, Sum, product

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("network")
x = m.var("x", (P, W))
rows = m.eq("capacity", Sum(W, x[P, W]) <= 10.0, over=product((P,)))

print(rows.n_rows)
```

Output:

```text
2
```

## `over` or `where`, not both

`over=` states the rows and `where=` restricts them, so passing both raises
`ValueError`.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum, product, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

m = Model("network")
x = m.var("x", (P, W))

m.eq(
    "capacity",
    Sum(W, x[P, W]) <= 10.0,
    where=subset((P,), {"P": np.array(["p0"])}),
    over=product((P,)),
)
```

Raises ValueError:

```text
ValueError: constraint 'capacity' states its rows with over= and narrows them with where=; state one
```

A condition on a sum and a condition on the constraint compose: the first
restricts what is summed, the second which rows exist.

---

# /guides/fixed-members

# A member fixed at a label

Initial conditions, terminal conditions and boundary rows reference one
member of a set: the state at the first period, the level at the last. A
label in place of a set in a reference fixes that dimension at one member
and removes it from the frame.

```python
import numpy as np
from nimopt import Model, Set

G = Set("G", np.array(["g0", "g1"]))
T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (G, T))

print(x[G, T].frame)
print(x[G, "t0"].frame)
```

Output:

```text
('G', 'T')
('G',)
```

`T` is fixed at `t0`, so the reference is indexed over `G` alone. An
initial condition is one row per unit, referencing that unit's column at the
first period.

```python
import numpy as np
from nimopt import Model, Set

G = Set("G", np.array(["g0", "g1"]))
T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (G, T))
rows = m.eq("start", x[G, "t0"] <= 1.0)

print(rows.n_rows, rows.nnz)
print(m.assemble().to_dense())
```

Output:

```text
2 2
[[1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0.]]
```

Two rows over six columns, each with one nonzero: the `t0` column of its own
unit.

## A coefficient at a member

A parameter takes a label the same way and yields the coefficients at that
member.

```python
import numpy as np
from nimopt import Model, Param, Set

G = Set("G", np.array(["g0", "g1"]))
T = Set("T", np.array(["t0", "t1", "t2"]))
rate = Param.from_dense("rate", (G, T), np.array([[2.0, 4.0, 5.0], [3.0, 1.0, 6.0]]))

m = Model("schedule")
x = m.var("x", (G, T))
rows = m.eq("start", rate[G, "t0"] * x[G, "t0"] <= 1.0)

print(rate[G, "t0"].dims)
print(m.assemble().to_dense())
```

Output:

```text
('G',)
[[2. 0. 0. 0. 0. 0.]
 [0. 0. 0. 3. 0. 0.]]
```

The coefficients are the `t0` column of `rate`: 2.0 and 3.0.

## The label must be a member

A label that is not a member of the set raises `ValueError`; the message
gives the label and the dimension.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set

G = Set("G", np.array(["g0", "g1"]))
T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (G, T))

x[G, "t9"]
```

Raises ValueError:

```text
ValueError: variable 'x' is read at member 't9' of dimension 'T', which that set does not carry
```

---

# /guides/highs-methods

# Interior point and first-order methods

`method=` names the algorithm HiGHS runs on the assembled matrix. Four
values name one each, and three further options shape the interior point
and first-order methods. Every setting on this page runs in the test
suite against the bundled models, and a solve carrying it reports the
optimum the model has; which method is fastest or smallest on a given
model is the model's to show.

| Option | Choices | Applies to |
| --- | --- | --- |
| `method` | `choose`, `simplex`, `barrier`, `hipo`, `pdlp` | every solve |
| `newton_system` | `choose`, `augmented`, `normaleq` | `hipo` |
| `crossover` | `choose`, `off`, `on` | `barrier` and `hipo` |
| `pdlp_tol` | a relative tolerance | `pdlp` |

`barrier` runs IPX, HiGHS's interior point method with a preconditioned
conjugate gradient at its core. `hipo` runs HiPO, an interior point method
built on a direct factorisation of the Newton system, parallel across the
elimination tree; `newton_system` chooses between the augmented system and
the normal equations, and `choose` leaves that to the solver. HiPO holds
that factorisation in memory, so it needs more of it than IPX on the same
model, and `threads` reaches its iterations but not the crossover that may
follow them. `pdlp` runs cuPDLP-C, a primal-dual hybrid gradient method that
touches the matrix only through matrix-vector products. `crossover` decides
whether an interior point method hands its solution to the simplex method
to reach a vertex.

```python skip="needs a HiGHS built with HiPO"
solution = model.solve(
    options={
        "method": "hipo",
        "newton_system": "augmented",
        "crossover": "off",
        "threads": 8,
    }
)
```

Gurobi and Mosek carry `crossover` and `method` up to `barrier`.
`newton_system`, `pdlp_tol`, `hipo` and `pdlp` are HiGHS's, and asking either
of the others for one of them is refused by name, as the
[solvers reference](/reference/solvers) states.

## Crossover

An interior point method stops at a point inside the feasible region,
within tolerance of the optimum on every constraint. Crossover moves that
point to a vertex with the simplex method, which is what a basis, an exact
active set and duals at a vertex require. It runs serially, and on a large
model it can cost more than the interior point iterations before it.

`crossover="off"` returns the interior point as the solution. Primals and
duals are read back the same way; what changes is that a constraint holding
with equality at the optimum may sit a tolerance away from it, and a
variable at a bound may sit a tolerance inside it. A model whose answer is
read as quantities and prices, rather than as a basis, is served by the
interior point.

## PDLP on a GPU

A first-order method holds no factorisation. Its memory is the matrix in
two orientations, one for each product, and a set of working vectors of the
row and column dimensions, so it grows linearly with the problem and a
problem too large for a factorisation still fits.

**Where the memory goes.** HiGHS presolves on the CPU, in host memory,
before any method runs: the original LP and its reduced form are both held
there while presolve works, which is the host's peak. The reduced problem is
what moves to the GPU, and it is smaller than the problem stated, because
presolve removes the rows that only bound a single column and the columns
presolve can fix. The GPU therefore needs memory for the reduced matrix
twice, the working vectors, and the sparse kernels' buffers, and nothing
else. The host keeps the original and the reduced LP throughout, and
postsolve maps the answer back through them.

A HiGHS built without CUDA runs the same method on the CPU. The solver's
log, asked for with `log=True`, names the device it runs on.

**The tolerance.** `pdlp_tol` is the relative tolerance at which PDLP
stops: the duality gap and the primal and dual residuals, each relative to
the scale of the problem, must all fall below it. A looser tolerance stops
sooner and returns a point farther from the optimum and farther from
feasibility; a tighter one costs more iterations, and each iteration is a
pass over the matrix. The method's iterations are cheap and numerous, so
the tolerance is the lever on both the time and the quality of the answer.

HiGHS checks the point PDLP returns against its own `feasibility_tol` and
`optimality_tol` after postsolve, and that check is stricter than PDLP's
own criterion: PDLP measures its residuals relative to the scaled problem
it iterates on, and postsolve maps the point back onto the original rows,
where a residual that passed can exceed the tolerance. A point that meets
`pdlp_tol` and misses HiGHS's check reports the status HiGHS calls unknown,
and `nimopt` refuses that status rather than reading values from it: the
solve raises, naming the status. On a model of any size, the default
`pdlp_tol` is therefore usually not enough for HiGHS to accept the point,
and a tolerance one to two orders tighter is where the two agree.
Tightening `pdlp_tol` until the point passes, or loosening the feasibility
and optimality tolerances to what the model needs, are the two ways
through.

PDLP reports two iterates in its log, the running average marked `[A]` and
the last marked `[L]`, and stops on whichever meets the tolerance first.
Progress is not monotone: the gap can close and open again as the method
restarts, so the last row of the log, not the best one, is what it returns.

```python skip="needs a HiGHS built with CUDA to run on a GPU"
solution = model.solve(options={"method": "pdlp", "pdlp_tol": 1e-8, "log": True})
```

## Installing a HiGHS that carries HiPO and a GPU

HiGHS 1.15.1 keeps HiPO's orderings and its BLAS in a library of its own,
`libhighs_extras`, which `libhighs` loads at run time by name. The `highspy`
wheel on PyPI and the conda-forge package ship without that library and
without CUDA. Asked for `hipo`, such a HiGHS logs an error and runs simplex;
`nimopt` refuses the request instead, naming the missing library. Asked for
`pdlp`, it runs the method on the CPU.

Both come from a source build of the same version, with three pieces
installed separately:

1. **The extras library**, from the `extern` directory of the HiGHS
   repository, which is a CMake project of its own. It builds METIS, AMD
   and RCM from the tree and links a BLAS; OpenBLAS from a conda
   environment serves.

   ```bash
   cmake -S extern -B build-extras -G Ninja -DCMAKE_BUILD_TYPE=Release \
     -DHIPO=ON -DBLA_VENDOR=OpenBLAS
   cmake --build build-extras
   ```

2. **The wheel**, from the repository root, where HiGHS's own build
   configuration turns HiPO on. The conda environment's prefix is given
   through the `CMAKE_PREFIX_PATH` environment variable, not through
   `CMAKE_ARGS`, because the latter replaces the path the Python build adds
   for pybind11.

   ```bash
   CMAKE_PREFIX_PATH=$CONDA_PREFIX CMAKE_ARGS="-DBLA_VENDOR=OpenBLAS" \
     uv build --wheel --python <venv>/bin/python -o dist .
   ```

   For the GPU, the same command with a CUDA toolkit in the environment and
   the card's compute capability:

   ```bash
   CMAKE_PREFIX_PATH=$CONDA_PREFIX CUDACXX=$CONDA_PREFIX/bin/nvcc CUDAToolkit_ROOT=$CONDA_PREFIX \
   CMAKE_ARGS="-DCUPDLP_GPU=ON -DCMAKE_CUDA_ARCHITECTURES=75 -DBLA_VENDOR=OpenBLAS" \
     uv build --wheel --python <venv>/bin/python -o dist .
   ```

3. **The placement.** The wheel's `libhighs` searches its own directory
   for the extras library, so the built `libhighs_extras.so` goes into the
   `highspy` package directory of the environment the wheel is installed
   in. The extras library links the BLAS it was built against, and the GPU
   wheel links the CUDA runtime, cuBLAS and cuSPARSE, so the conda
   environment that provided them stays. The GPU wheel's `libhighs` and
   `libcudalin` must name that environment's `lib` in their own run path,
   which `patchelf --set-rpath` sets on the unpacked wheel before it is
   packed again; the extension module's run path does not reach the CUDA
   libraries through them.

   ```bash
   uv pip install --python <venv>/bin/python dist/highspy-1.15.1-*.whl
   cp build-extras/libhighs_extras.so <venv>/lib/python3.13/site-packages/highspy/
   ```

The check that HiPO is in place is a solve asking for it: a HiGHS without
the extras library is refused before anything is solved, and one with it
returns the optimum.

---

# /guides/lags

# Lags

Time-coupled constraints reference the previous period. A storage balance
relates the state of charge at `t` to that at `t-1`; a ramp limit bounds
the change in output between consecutive periods. `T - 1` is the set `T`
lagged by one member, and `x[T - 1]` references the variable at the
previous member.

## A lag that drops the boundary row

```python
import numpy as np
from nimopt import Model, Set

T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (T,))
rows = m.eq("carry", x[T] - x[T - 1] <= 0.0)

print(rows.n_rows, rows.nnz)
print(m.assemble().to_dense())
```

Output:

```text
2 4
[[-1.  1.  0.]
 [ 0. -1.  1.]]
```

Each row references its own column and the previous one. The first member
has no predecessor, so its row is not produced: three members give two rows.

`T + 1` references the following member by the same rule.

## A lag that wraps

`T.cyclic` lags with wrap-around: the member before the first is the last.
No row is dropped.

```python
import numpy as np
from nimopt import Model, Set

T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (T,))
rows = m.eq("carry", x[T] - x[T.cyclic - 1] <= 0.0)

print(rows.n_rows, rows.nnz)
print(m.assemble().to_dense())
```

Output:

```text
3 6
[[ 1.  0. -1.]
 [-1.  1.  0.]
 [ 0. -1.  1.]]
```

Three rows, and the first references the last column: the `-1` in row 0 is
in the final position. A storage balance over a repeating horizon is written
this way, so that the level at the end of the horizon carries into the
beginning.

## A lag applies to a reference, not to a sum

`Sum` runs over the members of a set and takes the set itself. Passing a
lagged set raises `ValueError`.

```python raises=ValueError
import numpy as np
from nimopt import Model, Set, Sum

T = Set("T", np.array(["t0", "t1", "t2"]))

m = Model("schedule")
x = m.var("x", (T,))

Sum(T - 1, x[T])
```

Raises ValueError:

```text
ValueError: a sum is over the members of ['T'], so it takes the set and not a lag of it; state the lag at the variable's reference
```

The lag belongs on the variable reference: `Sum(T, x[T - 1])`.

## A lag applies to a variable, not to a parameter

Reading a parameter at a lag raises `ValueError`. A coefficient is indexed
by the row it appears in, and a lag selects which column a row references.
`rate[T] * x[T - 1]` applies the rate at `t` to the variable at `t-1`.

```python raises=ValueError
import numpy as np
from nimopt import Param, Set

T = Set("T", np.array(["t0", "t1", "t2"]))
rate = Param.from_dense("rate", (T,), np.array([1.0, 2.0, 3.0]))

rate[T - 1]
```

Raises ValueError:

```text
ValueError: parameter 'rate' is read at a lag ['T']; state the lag at the variable's reference, where a coefficient multiplies the row it lands on
```

A lag is an integer number of members. A fractional lag raises `ValueError`
rather than being truncated to a different lag.

```python raises=ValueError
import numpy as np
from nimopt import Set

T = Set("T", np.arange(3))
T - 1.7
```

Raises ValueError:

```text
ValueError: a lag is a whole number of members; got 1.7
```

---

# /guides/saving-and-loading

# Saving and loading a model

A definition writes itself to a YAML file whose expressions are spelled as
they are typed in Python. The file is the model's final form: every derived
coefficient spelled out, every term with its own sum and sign, the constant
last.

```python
from nimopt import Definition, Sum

d = Definition("dispatch", sense="min")
G, T = d.set("G"), d.set("T")
price, eta = d.param("price", (G, T)), d.param("eta", (G, T))
cap, load = d.param("cap", (G, T)), d.param("load", (T,))
gen = d.var("gen", (G, T), upper=cap)
d.eq("balance", Sum(G, gen[G, T]) == load[T])
d.set_objective(Sum(G, T, 2 * (price[G, T] / eta[G, T]) * gen[G, T]))

print(d.to_yaml())
```

Output:

```text
version: 2
name: dispatch
sense: min
sets: [G, T]
parameters:
  price: [G, T]
  eta: [G, T]
  cap: [G, T]
  load: [T]
variables:
  gen:
    sets: [G, T]
    upper: cap
constraints:
  balance:
    relation: Sum(G, gen[G, T]) == load[T]
objective: Sum(G, T, ((price[G, T] / eta[G, T]) * 2) * gen[G, T])
```

The structure section is the data's schema: every set, and every parameter
with its dimensions. `build` refuses a mapping that misses any of them by name.

## Reading a file back

`loads` reads text and `load` reads a path. A file without data returns a
`Definition`, whose next step is `build(data)`.

```python
import numpy as np
from nimopt import loads

text = """
version: 2
name: dispatch
sense: min
sets: [G, T]
parameters:
  cost: [G]
  load: [T]
variables:
  gen:
    sets: [G, T]
    upper: 10.0
constraints:
  balance:
    relation: Sum(G, gen[G, T]) == load[T]
objective: Sum(G, T, cost[G] * gen[G, T])
"""

d = loads(text)
m = d.build(
    {
        "G": np.array(["a", "b"]),
        "T": np.arange(2),
        "cost": np.array([1.0, 3.0]),
        "load": np.array([12.0, 15.0]),
    }
)
print(m.solve().objective)
```

Output:

```text
41.0
```

## Editing by hand

The text is read through the same operators a Python model is built from,
so an edit is accepted where Python accepts it and normalised the same way.
Adding a scalar, reordering terms, or writing a comparison the other way
round all read, and the file written back is the canonical form.

```python
from nimopt import loads

edited = """
version: 2
name: dispatch
sense: min
sets: [G, T]
parameters:
  cost: [G]
  load: [T]
variables:
  gen:
    sets: [G, T]
constraints:
  balance:
    relation: load[T] == Sum(G, gen[G, T]) * 2 + 1 - 1
objective: Sum(G, T, cost[G] * gen[G, T])
"""
print(loads(edited).to_yaml().splitlines()[-2])
```

Output:

```text
relation: 2 * Sum(G, gen[G, T]) == load[T]
```

An edit Python refuses is refused here with the same message, and a
construct outside the spelling is refused naming it.

```python raises=ValueError
from nimopt import loads

loads(
    """
version: 2
name: dispatch
sense: min
sets: [G, T]
variables:
  gen:
    sets: [G, T]
constraints:
  peak:
    relation: max(gen[G, T]) <= 10
"""
)
```

Raises ValueError:

```text
ValueError: 'max(gen[G, T]) <= 10': Sum is the one call the spelling carries
```

## Data inline, for a model small enough to read

A built model writes its data into the file with `inline=True`. A set is a
list, a dense parameter is nested lists, and a parameter carrying some
coordinates of its product is a table of the dimensions then `value`. A
file carrying data loads to a built `Model`.

```python
from nimopt import loads
from nimopt.models import transport

m = transport.definition().build(transport.data())
text = m.to_yaml(inline=True)
print(text[text.index("data:") :])
print(loads(text).solve().objective)
```

Output:

```text
data:
  P: [p0, p1, p2, p3]
  W: [w0, w1, w2, w3, w4, w5]
  cost:
    columns: [P, W, value]
    rows:
    - [p0, w0, 1.1322210842282328]
    - [p0, w1, 7.506161913602179]
    - [p0, w4, 1.3277881914895575]
    - [p1, w0, 6.835972487871987]
    - [p1, w3, 8.302044618221775]
    - [p1, w4, 5.853086206137439]
    - [p2, w2, 5.348999931723383]
    - [p2, w3, 8.480579390302147]
    - [p2, w4, 7.526828432972257]
    - [p3, w1, 1.0219080013611848]
    - [p3, w2, 7.859234212700555]
    - [p3, w3, 1.2686846024437148]
  supply: [60.0, 60.0, 60.0, 60.0]
  demand: [10.0, 10.0, 10.0, 10.0, 10.0, 10.0]

100.99601811246072
```

## Data beside the file, for a model of size

`save` writes a model's file and its data as an `.npz` beside it, under the
file's stem, and the file names that sidecar. `load` reads both. The pair
moves together; the file names no path and no machine.

```python
import tempfile
from pathlib import Path
from nimopt import load, save
from nimopt.models import storage

m = storage.definition().build(storage.data())
with tempfile.TemporaryDirectory() as held:
    path = Path(held) / "storage.yaml"
    save(m, path)
    print(sorted(p.name for p in Path(held).iterdir()))
    print(path.read_text().splitlines()[-1])
    print(load(path).solve().status)
```

Output:

```text
['storage.npz', 'storage.yaml']
data: storage.npz
optimal
```

A definition's file takes data from the caller instead: `load(path, data=...)`
with the mapping `build` takes or the path of an `.npz`. A file that carries
data and a `data=` together is refused, because two sources for one model is
a choice the library does not make.

---

# /guides/subsets

# A variable over a subset

In a network model, most pairs of a set product are not connected: a plant
serves some warehouses, a line joins two of many buses. A variable declared
over the full product has a column for every pair, including the ones that
do not exist. `subset=` restricts the variable to the members that do, and
the others have no column at all.

## Declaring the subset

`subset(sets, columns)` lists members of a set product by label, one column
per set, read in parallel. `m.var(..., subset=arcs)` declares the variable
over those members.

```python
import numpy as np
from nimopt import Model, Set, Sum, product, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

arcs = subset(
    (P, W),
    {"P": np.array(["p0", "p0", "p1"]), "W": np.array(["w0", "w1", "w2"])},
)

m = Model("network")
x = m.var("x", (P, W), subset=arcs)

print(product((P, W)).size)
print(x.n_columns)
print(m.assemble().to_dense())
```

Output:

```text
6
3
[]
```

The product has six members and the subset three, so `x` has three columns.
The assembled matrix is empty because no constraint has been added. A model
over a sparse network pays for its arcs, not for the grid that contains
them.

## By label or by position

`subset` takes labels and `subset_of` takes integer positions. Both read
their columns in parallel: the k-th entry of each column belongs to the same
member. A subset is a list of members, not a cross product of its columns.

```python
import numpy as np
from nimopt import Set, subset, subset_of

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

by_label = subset(
    (P, W),
    {"P": np.array(["p0", "p1"]), "W": np.array(["w0", "w2"])},
)
by_index = subset_of((P, W), np.array([[0, 1], [0, 2]]))

print(by_label.size, by_index.size)
print(by_label.labels())
```

Output:

```text
2 2
{'P': array(['p0', 'p1'], dtype='<U2'), 'W': array(['w0', 'w2'], dtype='<U2')}
```

Both list the same two members, `(p0, w0)` and `(p1, w2)`. `subset_of`
avoids resolving labels when positions are already at hand.

## Constraints over a subset variable

A sum over a subset variable runs over the members the variable has, so a
row contains the arcs at that member and nothing else.

```python
import numpy as np
from nimopt import Model, Set, Sum, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

arcs = subset(
    (P, W),
    {"P": np.array(["p0", "p0", "p1"]), "W": np.array(["w0", "w1", "w2"])},
)

m = Model("network")
x = m.var("x", (P, W), subset=arcs)
rows = m.eq("capacity", Sum(W, x[P, W]) <= 10.0)

print(rows.n_rows, rows.nnz)
print(m.assemble().to_dense())
```

Output:

```text
2 3
[[1. 1. 0.]
 [0. 0. 1.]]
```

Two rows over three columns: `p0` has two arcs and `p1` one.

## Bounds over a subset

A bound applies to every column of the variable. A parameter indexed over
fewer dimensions than the variable is broadcast over the rest, so a bound
per plant applies to each of that plant's arcs.

```python
import numpy as np
from nimopt import Model, Param, Set, subset

P = Set("P", np.array(["p0", "p1"]))
W = Set("W", np.array(["w0", "w1", "w2"]))

arcs = subset(
    (P, W),
    {"P": np.array(["p0", "p0", "p1"]), "W": np.array(["w0", "w1", "w2"])},
)
cap = Param.from_dense("cap", (P,), np.array([4.0, 9.0]))

m = Model("network")
m.var("x", (P, W), subset=arcs, upper=cap)

print(m.column_bounds()[1])
```

Output:

```text
[4. 4. 9.]
```

The two `p0` arcs take 4.0 and the `p1` arc 9.0.

---

# /models/commitment

# Commitment

`nimopt.models.commitment` is unit commitment. A committed unit runs between
its minimum and its maximum output and pays a no-load cost for being on; an
uncommitted unit produces nothing. The `capacity` and `minimum` rows are
written against the binary column. This is the corpus's only MILP.

```text
minimise    Σ_{t,g} cost[g] · gen[t,g] + Σ_{t,g} no_load[g] · on[t,g]
subject to  gen[t,g] − p_max[g] · on[t,g] ≤ 0
            gen[t,g] − p_min[g] · on[t,g] ≥ 0
            Σ_g gen[t,g] == load[t]           for each snapshot t
            on[t,g] ∈ {0, 1},  gen[t,g] ≥ 0
```

```python
from nimopt.models import commitment

print(commitment.definition().explain())
```

Output:

```text
commitment  min  not built
  sets        T · G
  parameters  p_max (G) · p_min (G) · cost (G) · no_load (G) · load (T)
  variables   on (T×G) [0.0, 1.0] integer · gen (T×G) [0.0, inf]
  constraint  capacity (T,G)  gen[T, G] - p_max[G] * on[T, G] <= 0
  constraint  minimum (T,G)  gen[T, G] - p_min[G] * on[T, G] >= 0
  constraint  balance (T)  Sum(G, gen[T, G]) == load[T]
  objective   min  Sum(T, G, cost[G] * gen[T, G]) + Sum(T, G, no_load[G] * on[T, G])
```

Snapshots are uncoupled, so `reference` enumerates every on-off subset per
snapshot and takes the cheapest feasible one: an optimum computed without a
solver.

```python
from nimopt.models import commitment

inputs = commitment.data()
model = commitment.definition().build(inputs)
solution = model.solve()
print(int(model.integrality().sum()), "binary columns of", model.n_columns)
print(solution.objective, commitment.reference(inputs))
```

Output:

```text
12 binary columns of 24
13800.0 13800.0
```

Both unit rows are produced for every generator and snapshot.

```python
from nimopt.models import commitment

model = commitment.definition().build(commitment.data())
print(model.absent("capacity"))
```

Output:

```text
capacity  12 of 12 rows  stated by terms
```

---

# /models/dispatch

# Dispatch

`nimopt.models.dispatch` is least-cost dispatch of a generator fleet against
a load. One variable `p` is indexed over snapshots and generators, there is
one balance row per snapshot, and each generator has a cost. Every other
model in the corpus adds one axis to this one.

```text
minimise    Σ_{t,g} cost[g] · p[t,g]
subject to  Σ_g p[t,g] == load[t]        for each snapshot t
            0 ≤ p[t,g] ≤ p_max[g]
```

The balance row has no coefficient: a sum over a dimension needs none, and
the corpus writes no coefficient a model does not need.

```python
from nimopt.models import dispatch

print(dispatch.definition().explain())
```

Output:

```text
dispatch  min  not built
  sets        snapshot · generator
  parameters  p_max (generator) · load (snapshot) · cost (generator)
  variables   p (snapshot×generator) [0.0, p_max]
  constraint  balance (snapshot)  Sum(generator, p[snapshot, generator]) == load[snapshot]
  objective   min  Sum(snapshot, generator, cost[generator] * p[snapshot, generator])
```

Snapshots are independent, so the optimum is the merit order per snapshot,
and `reference` computes it without a solver.

```python
from nimopt.models import dispatch

inputs = dispatch.data()
solution = dispatch.definition().build(inputs).solve()
print(solution.status)
print(solution.objective, dispatch.reference(inputs))
```

Output:

```text
optimal
1920.0 1920.0
```

The balance row is produced for every snapshot in the load, so `absent`
reports no dropped rows.

```python
from nimopt.models import dispatch

model = dispatch.definition().build(dispatch.data())
print(model.absent("balance"))
```

Output:

```text
balance  6 of 6 rows  stated by terms
```

---

# /models/expansion

# Expansion

`nimopt.models.expansion` is a two-stage stochastic program. The first stage
builds capacity; the second stage dispatches it against a demand and a fuel
price that the scenario reveals. What the first stage cannot see is stated by
the shape of its variable: `cap` is indexed by technology alone, so one
capacity serves every scenario, while `p` and `shed` carry the scenario and
may differ across it.

```text
minimise    Σ_g capital[g] · cap[g]
            + Σ_{s,g,t} weight[s] · cost[s,g] · p[s,g,t]
            + Σ_{s,t}   weight[s] · voll[s]  · shed[s,t]
subject to  p[s,g,t] ≤ cap[g]                              for each s, g, t
            Σ_g p[s,g,t] + shed[s,t] == demand[s,t]        for each s, t
            cap, p, shed ≥ 0
```

`weight` is the probability of a scenario, so the second and third sums are
an expectation and the objective is the capital committed plus the expected
cost of the recourse. There is no row tying one scenario's capacity to
another's: non-anticipativity is the missing dimension.

```python
from nimopt.models import expansion

print(expansion.definition().explain())
```

Output:

```text
expansion  min  not built
  sets        S · G · T
  parameters  capital (G) · cost (S,G) · demand (S,T) · weight (S) · voll (S)
  variables   cap (G) [0.0, inf] · p (S×G×T) [0.0, inf] · shed (S×T) [0.0, inf]
  constraint  capacity (S,G,T)  p[S, G, T] - cap[G] <= 0
  constraint  balance (S,T)  Sum(G, p[S, G, T]) + shed[S, T] == demand[S, T]
  objective   min  Sum(G, capital[G] * cap[G]) + Sum(S, G, T, (weight[S] * cost[S, G]) * p[S, G, T]) + Sum(S, T, (weight[S] * voll[S]) * shed[S, T])
```

`cap` has one column per technology and `p` has one per scenario, technology
and hour. The capacity row is stated over all three dimensions even though
the variable it bounds carries one, so a single column is read by every
scenario's rows — which is what makes the capacity a shared decision.

```python
from nimopt.models import expansion

model = expansion.definition().build(expansion.data())
print(model.explain())
```

Output:

```text
expansion  min  75 columns · 72 rows · 180 nonzeros
  sets        G 3 · S 3 · T 6
  parameters  demand (S,T) 18 · capital (G) 3 · weight (S) 3 · cost (S,G) 9 · voll (S) 3
  variables   cap (G) 3 cols [0.0, inf] · p (S×G×T) 54 cols [0.0, inf] · shed (S×T) 18 cols [0.0, inf]
  constraint  capacity (S,G,T)  p[S, G, T] - cap[G] <= 0  54 rows  108 nz
  constraint  balance (S,T)  Sum(G, p[S, G, T]) + shed[S, T] == demand[S, T]  18 rows  72 nz
  objective   min  Sum(G, capital[G] * cap[G]) + Sum(S, G, T, (weight[S] * cost[S, G]) * p[S, G, T]) + Sum(S, T, (weight[S] * voll[S]) * shed[S, T])
```

A technology is available in full wherever it is built and nothing couples
one hour to the next, so the recourse in each scenario-hour is the merit
order of the built capacity against that demand, with the remainder unserved
at `voll`. Write the capacity as bands — the cheapest technology's band, then
the next, and last the band between the dearest technology and lost load —
and the total separates into one term per band, each convex in that band's
level and turning only at a demand. `reference` minimises them one at a time
and adds them up, so the answer is arithmetic over the inputs rather than a
second solve.

```python
from nimopt.models import expansion

inputs = expansion.data()
solution = expansion.definition().build(inputs).solve()
print(solution.objective, expansion.reference(inputs))
```

Output:

```text
29703.899999999994 29703.899999999994
```

The data states a cost frontier: capital falls as marginal cost rises, and
lost load is dearer than the dearest technology. All three are built, and the
capacity that stacks to 180 leaves the cold scenario's peak hour of 225
short — shedding 45 is cheaper than a fourth band that earns its capital in
one hour of one scenario.

```python
import numpy as np

from nimopt.models import expansion

inputs = expansion.data()
solution = expansion.definition().build(inputs).solve()
print(np.cumsum(np.asarray(solution.primal("cap").values())))
print(np.asarray(solution.primal("shed").values()).reshape(3, 6))
```

Output:

```text
[119. 153. 180.]
[[ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0. 45.  0.  0.]]
```

The separation holds only for data whose technologies are a frontier and
whose bands stack. `reference` refuses anything else rather than returning a
number that is not the optimum.

```python raises=ValueError
from nimopt.models import expansion

inputs = expansion.data()
inputs["capital"] = inputs["capital"][::-1]
print(expansion.reference(inputs))
```

Raises ValueError:

```text
ValueError: capital does not fall from base to what follows it
```

---

# /models/fleet

# Fleet

`nimopt.models.fleet` is the same problem as `dispatch`, declared as one
variable per unit over the snapshots alone, with the units' terms added into
one balance row. The answer is the same merit order; what differs is the
cost of declaring it.

```text
minimise    Σ_t Σ_u cost_u[t] · u[t]
subject to  Σ_u u[t] == load[t]          for each snapshot t
            0 ≤ u[t] ≤ p_max_u[t]        for each unit u
```

`definition` takes a scale here, because the number of variables is a
property of the declaration rather than of the data.

```python
from nimopt.models import fleet

print(fleet.definition().explain())
```

Output:

```text
fleet  min  not built
  sets        T
  parameters  load (T) · p_max_g0_0 (T) · cost_g0_0 (T) · p_max_g1_0 (T) · cost_g1_0 (T) · p_max_g2_0 (T) · cost_g2_0 (T)
  variables   g0_0 (T) [0.0, p_max_g0_0] · g1_0 (T) [0.0, p_max_g1_0] · g2_0 (T) [0.0, p_max_g2_0]
  constraint  balance (T)  g0_0[T] + g1_0[T] + g2_0[T] == load[T]
  objective   min  Sum(T, cost_g0_0[T] * g0_0[T]) + Sum(T, cost_g1_0[T] * g1_0[T]) + Sum(T, cost_g2_0[T] * g2_0[T])
```

```python
from nimopt.models import fleet

inputs = fleet.data()
model = fleet.definition().build(inputs)
solution = model.solve()
print(len(model.variables), "variables,", model.n_columns, "columns")
print(solution.objective, fleet.reference(inputs))
```

Output:

```text
3 variables, 12 columns
8900.0 8900.0
```

One balance row per hour, and every unit appears in every row.

```python
from nimopt.models import fleet

model = fleet.definition().build(fleet.data())
print(model.absent("balance"))
```

Output:

```text
balance  4 of 4 rows  stated by terms
```

---

# /models

# Worked models

`nimopt.models` contains ten models. Each is a module with three functions.

| Name | Returns |
| --- | --- |
| `definition()` | a `Definition`: the formulation with no data bound |
| `data(scale=1)` | the inputs, at a given size |
| `reference(data)` | the optimal objective, computed by direct arithmetic |

The three serve different readers. This documentation calls `explain()`, a
benchmark calls `build(data(100))`, and a test compares a solve against
`reference(data(1))`. A reference is arithmetic over the inputs that uses
nothing from `nimopt`, so a formulation error is not checked against a copy
of itself.

| Model | Exercises |
| --- | --- |
| [`dispatch`](/models/dispatch) | the baseline formulation |
| [`transport`](/models/transport) | a sparse network over a subset of a product |
| [`storage`](/models/storage) | temporal coupling and a cyclic state |
| [`nodal`](/models/nodal) | grouping through a lookup parameter |
| [`commitment`](/models/commitment) | binary columns |
| [`fleet`](/models/fleet) | many small declarations rather than one large one |
| [`profiled`](/models/profiled) | a bound that varies by hour |
| [`sector`](/models/sector) | mixed density: dense in one axis, sparse in another |
| [`expansion`](/models/expansion) | a two-stage stochastic program: capacity before the scenario, dispatch after |
| [`recourse`](/models/recourse) | a binary first stage taken before the scenario is known |

```python
from nimopt.models import dispatch

inputs = dispatch.data()
solution = dispatch.definition().build(inputs).solve()
print(solution.status, solution.objective, dispatch.reference(inputs))
```

Output:

```text
optimal 1920.0 1920.0
```

---

# /models/nodal

# Nodal

`nimopt.models.nodal` groups generators into buses through a lookup
parameter. `at[G, B]` has an entry where generator `g` sits at bus `b`, and
multiplying the generation by it maps a row over generators to a row over
buses. The coefficient introduces `B`, a dimension no variable has, so the
balance is indexed over the dimensions the lookup defines.

```text
minimise    Σ_{t,g} cost[g] · gen[t,g]
subject to  Σ_g at[g,b] · gen[t,g] == demand[b,t]    for each bus b and hour t
            0 ≤ gen[t,g] ≤ p_max[g]
```

```python
from nimopt.models import nodal

print(nodal.definition().explain())
```

Output:

```text
nodal  min  not built
  sets        T · G · B
  parameters  at (G,B) · p_max (G) · cost (G) · demand (B,T)
  variables   gen (T×G) [0.0, p_max]
  constraint  balance (B,T)  Sum(G, at[G, B] * gen[T, G]) == demand[B, T]
  objective   min  Sum(T, G, cost[G] * gen[T, G])
```

Each bus meets its own demand from the generators sited at it, so the
optimum is a merit order per bus and hour.

```python
from nimopt.models import nodal

inputs = nodal.data()
solution = nodal.definition().build(inputs).solve()
print(solution.objective, nodal.reference(inputs))
```

Output:

```text
11850.0 11850.0
```

Every bus-hour has a row. A generator sited elsewhere is a term the lookup
removes from that row, not a row that is dropped.

```python
from nimopt.models import nodal

model = nodal.definition().build(nodal.data())
print(model.absent("balance"))
```

Output:

```text
balance  6 of 6 rows  stated by terms
  term absent G='g0_0', B='b1_0', T=0  gen  absent-coefficient (at)
  term absent G='g0_0', B='b1_0', T=1  gen  absent-coefficient (at)
  term absent G='g0_0', B='b1_0', T=2  gen  absent-coefficient (at)
  term absent G='g1_0', B='b1_0', T=0  gen  absent-coefficient (at)
  term absent G='g1_0', B='b1_0', T=1  gen  absent-coefficient (at)
  term absent G='g1_0', B='b1_0', T=2  gen  absent-coefficient (at)
  term absent G='g2_0', B='b0_0', T=0  gen  absent-coefficient (at)
  term absent G='g2_0', B='b0_0', T=1  gen  absent-coefficient (at)
  term absent G='g2_0', B='b0_0', T=2  gen  absent-coefficient (at)
  term absent G='g3_0', B='b0_0', T=0  gen  absent-coefficient (at)
  term absent G='g3_0', B='b0_0', T=1  gen  absent-coefficient (at)
  term absent G='g3_0', B='b0_0', T=2  gen  absent-coefficient (at)
```

---

# /models/profiled

# Profiled

`nimopt.models.profiled` is a dispatch whose capacity varies by hour.
`dispatch` bounds a generator by a single number; here `p_max` is a
parameter over generators and snapshots, so a solar unit is bounded by its
hourly availability and a thermal unit by its rating.

```text
minimise    Σ_{t,g} cost[g] · gen[t,g]
subject to  Σ_g gen[t,g] == load[t]          for each snapshot t
            0 ≤ gen[t,g] ≤ profile[g,t]
```

The profile is indexed `(G, T)` and the variable `(T, G)`. A bound is read
in the dimension order of the variable it bounds, so both orderings select
the same columns.

```python
from nimopt.models import profiled

print(profiled.definition().explain())
```

Output:

```text
profiled  min  not built
  sets        T · G
  parameters  profile (G,T) · cost (G) · load (T)
  variables   gen (T×G) [0.0, profile]
  constraint  balance (T)  Sum(G, gen[T, G]) == load[T]
  objective   min  Sum(T, G, cost[G] * gen[T, G])
```

Each snapshot is independent, so the optimum is the merit order against that
hour's capacities.

```python
from nimopt.models import profiled

inputs = profiled.data()
solution = profiled.definition().build(inputs).solve()
print(solution.objective, profiled.reference(inputs))
```

Output:

```text
28119.536003699297 28119.536003699293
```

Every hour has a balance row. A generator whose profile is zero is a column
bounded to zero, not a dropped row.

```python
from nimopt.models import profiled

model = profiled.definition().build(profiled.data())
print(model.absent("balance"))
```

Output:

```text
balance  8 of 8 rows  stated by terms
```

---

# /models/recourse

# Recourse

`nimopt.models.recourse` is `commitment` under uncertainty. The demand and the
fuel price are a scenario, and the on-off decision is taken before either is
revealed: `on` is indexed by hour and unit alone, while `p` and `shed` carry
the scenario. One commitment has to serve every scenario, which is what makes
the binary a hedge rather than a schedule.

```text
minimise    Σ_{t,g}   no_load[g] · on[t,g]
            + Σ_{s,g,t} weight[s] · cost[s,g] · p[s,g,t]
            + Σ_{s,t}   weight[s] · voll[s]  · shed[s,t]
subject to  p[s,g,t] ≤ p_max[g] · on[t,g]                  for each s, g, t
            p[s,g,t] ≥ p_min[g] · on[t,g]                  for each s, g, t
            Σ_g p[s,g,t] + shed[s,t] == demand[s,t]        for each s, t
            on ∈ {0,1},  p, shed ≥ 0
```

The first stage costs what it costs whatever happens; the second and third
sums are weighted by the probability of a scenario, so the objective is a
commitment charge plus the expected cost of the recourse.

```python
from nimopt.models import recourse

print(recourse.definition().explain())
```

Output:

```text
recourse  min  not built
  sets        S · G · T
  parameters  p_max (G) · p_min (G) · no_load (G) · cost (S,G) · demand (S,T) · weight (S) · voll (S)
  variables   on (T×G) [0.0, 1.0] integer · p (S×G×T) [0.0, inf] · shed (S×T) [0.0, inf]
  constraint  capacity (S,G,T)  p[S, G, T] - p_max[G] * on[T, G] <= 0
  constraint  minimum (S,G,T)  p[S, G, T] - p_min[G] * on[T, G] >= 0
  constraint  balance (S,T)  Sum(G, p[S, G, T]) + shed[S, T] == demand[S, T]
  objective   min  Sum(T, G, no_load[G] * on[T, G]) + Sum(S, G, T, (weight[S] * cost[S, G]) * p[S, G, T]) + Sum(S, T, (weight[S] * voll[S]) * shed[S, T])
```

Nothing couples one hour to the next, so the commitment is chosen hour by
hour. `reference` enumerates every on-off subset of the fleet and scores each
by its expected recourse across the scenarios, which is exact and cheap:
three units make eight subsets.

```python
from nimopt.models import recourse

inputs = recourse.data()
solution = recourse.definition().build(inputs).solve()
print(solution.objective, recourse.reference(inputs))
```

Output:

```text
16744.8125 16744.8125
```

A committed unit runs at least its minimum and there is nowhere to put
unwanted energy, so a unit whose minimum exceeds the mildest scenario's
demand cannot be committed at all — the row it would break is that scenario's
balance. The first hour asks 38.25 in the mildest scenario, and `base`, whose
minimum is 40.0 and whose fuel is the cheapest of the three, is left off it.
The third hour goes the other way: the fleet carries 260.0 against a coldest
demand of 268.75, so everything is committed and the remainder is shed.

```python
import numpy as np

from nimopt.models import recourse

inputs = recourse.data()
solution = recourse.definition().build(inputs).solve()
print(inputs["demand"].round(2))
print(np.asarray(solution.primal("on").values()).reshape(4, 3))
print(np.asarray(solution.primal("shed").values()).reshape(3, 4))
```

Output:

```text
[[ 38.25 119.   182.75  80.75]
 [ 45.   140.   215.    95.  ]
 [ 56.25 175.   268.75 118.75]]
[[0. 1. 0.]
 [1. 1. 0.]
 [1. 1. 1.]
 [1. 0. 0.]]
[[0.   0.   0.   0.  ]
 [0.   0.   0.   0.  ]
 [0.   0.   8.75 0.  ]]
```

The mean demand in the first hour is 43.875, which `base` can serve and is
cheapest at. A model given that one number commits it, and that commitment
breaks the balance of the mildest scenario, whose probability is 0.5. The
scenario dimension on `demand` is what rules the commitment out; the missing
scenario dimension on `on` is what makes the ruling out bind.

---

# /models/sector

# Sector

`nimopt.models.sector` has mixed density. The region-technology map is
sparse, since a technology exists in some regions and not others, while
every sited pair runs in every hour. The generation variable takes its
members from a parameter over the sited pairs crossed with the whole
horizon, so it is sparse in one axis and dense in the other.

```text
minimise    Σ_{(r,k) sited, t} cost[r,k] · gen[r,k,t]
subject to  Σ_k gen[r,k,t] == demand[r,t]     for each region r and hour t
            0 ≤ gen[r,k,t] ≤ capacity[r,k]    for each sited (r,k) and hour t
```

```python
from nimopt.models import sector

print(sector.definition().explain())
```

Output:

```text
sector  min  not built
  sets        R · K · T
  parameters  sited (R,K,T) · capacity (R,K) · cost (R,K) · demand (R,T)
  variables   gen (R×K×T) over sited [0.0, capacity]
  constraint  balance (R,T)  Sum(K, gen[R, K, T]) == demand[R, T]
  objective   min  Sum(R, K, T, cost[R, K] * gen[R, K, T])
```

Each region meets its own demand from the technologies sited in it, so the
optimum is a merit order per region and hour.

```python
from nimopt.models import sector

inputs = sector.data()
model = sector.definition().build(inputs)
solution = model.solve()
print(
    model.n_columns,
    "columns of a possible",
    len(inputs["R"]) * len(inputs["K"]) * len(inputs["T"]),
)
print(solution.objective, sector.reference(inputs))
```

Output:

```text
16 columns of a possible 24
18350.0 18350.0
```

Every region-hour has a balance row, and the capacity bound applies only to
the sited pairs, so no balance row is dropped.

```python
from nimopt.models import sector

model = sector.definition().build(sector.data())
print(model.absent("balance"))
```

Output:

```text
balance  8 of 8 rows  stated by terms
```

---

# /models/storage

# Storage

`nimopt.models.storage` dispatches a generator fleet and a set of batteries
against an hourly load. The state-of-charge row references the previous
hour through `T.cyclic - 1`, so the row at the first hour references the
last hour and every hour has a row. The ramp row references `T - 1`; the
first hour has no predecessor, so that row is not produced. The model
exercises both lag rules.

```text
minimise    Σ_{g,t} cost[g,t] · gen[g,t]
subject to  Σ_g gen[g,t] + Σ_s discharge[s,t] − Σ_s charge[s,t] == load[t]
            soc[s,t] − soc[s,t−1] − charge_eta[s,t] · charge[s,t]
                + discharge_eta[s,t] · discharge[s,t] == 0     (t−1 wraps)
            gen[g,t] − gen[g,t−1] ≤ ramp_limit[g,t]              (t=0 dropped)
            gen[g,t] ≤ capacity[g,t]
            charge[s,t] ≤ power[s,t]
            discharge[s,t] ≤ power[s,t]
            soc[s,t] ≤ energy[s,t]
```

Every limit is a constraint rather than a bound, so that each has a dual
value.

```python
from nimopt.models import storage

print(storage.definition().explain())
```

Output:

```text
storage  min  not built
  sets        T · G · S
  parameters  cost (G,T) · capacity (G,T) · ramp_limit (G,T) · load (T) · power (S,T) · energy (S,T) · charge_eta (S,T) · discharge_eta (S,T)
  variables   gen (G×T) [0.0, inf] · charge (S×T) [0.0, inf] · discharge (S×T) [0.0, inf] · soc (S×T) [0.0, inf]
  constraint  balance (T)  Sum(G, gen[G, T]) + Sum(S, discharge[S, T]) - Sum(S, charge[S, T]) == load[T]
  constraint  state_of_charge (S,T)  soc[S, T] - soc[S, T.cyclic - 1] - charge_eta[S, T] * charge[S, T] + discharge_eta[S, T] * discharge[S, T] == 0
  constraint  generation_limit (G,T)  gen[G, T] <= capacity[G, T]
  constraint  ramp (G,T)  gen[G, T] - gen[G, T - 1] <= ramp_limit[G, T]
  constraint  charge_limit (S,T)  charge[S, T] <= power[S, T]
  constraint  discharge_limit (S,T)  discharge[S, T] <= power[S, T]
  constraint  energy_limit (S,T)  soc[S, T] <= energy[S, T]
  objective   min  Sum(G, T, cost[G, T] * gen[G, T])
```

The batteries are lossy, with a round-trip efficiency of `0.95 · 0.93`, and
the fleet's costs span 50.0 to 55.0. Shifting energy through the store never
pays, so the store stays idle and the optimum is the hourly merit order. A
store that cycles requires data with a wider cost spread; the benchmarks
supply it.

```python
from nimopt.models import storage

inputs = storage.data()
solution = storage.definition().build(inputs).solve()
print(solution.objective, storage.reference(inputs))
print("store moved:", abs(solution.primal("charge").to_dense()).max())
```

Output:

```text
145449.3739227024 145449.37392270242
store moved: 0.0
```

Every hour has a state-of-charge row, because that lag wraps. The ramp row
at the first hour is not produced, because that lag does not.

```python
from nimopt.models import storage

model = storage.definition().build(storage.data())
print(model.absent("state_of_charge"))
print()
print(model.absent("ramp"))
```

Output:

```text
state_of_charge  24 of 24 rows  stated by terms

ramp  69 of 72 rows  stated by terms
  row absent  G='base0', T=0  term-does-not-reach (gen)
  row absent  G='mid0', T=0  term-does-not-reach (gen)
  row absent  G='peak0', T=0  term-does-not-reach (gen)
```

---

# /models/transport

# Transport

`nimopt.models.transport` ships from plants to warehouses over an incomplete
network: a plant serves a band of nearby warehouses rather than all of
them. The cost parameter has one entry per arc, and the flow variable takes
its members from that parameter, so the model has one column per arc rather
than one per cell of the plant-warehouse product.

```text
minimise    Σ_{(p,w) ∈ arcs} cost[p,w] · flow[p,w]
subject to  Σ_w flow[p,w] ≤ supply[p]     for each plant p
            Σ_p flow[p,w] ≥ demand[w]     for each warehouse w
            flow[p,w] ≥ 0                 for each arc (p,w)
```

```python
from nimopt.models import transport

print(transport.definition().explain())
```

Output:

```text
transport  min  not built
  sets        P · W
  parameters  cost (P,W) · supply (P) · demand (W)
  variables   flow (P×W) over cost [0.0, inf]
  constraint  supply (P)  Sum(W, flow[P, W]) <= supply[P]
  constraint  demand (W)  Sum(P, flow[P, W]) >= demand[W]
  objective   min  Sum(P, W, cost[P, W] * flow[P, W])
```

Supply is twice the total demand of a plant's band, so no supply row binds
and each warehouse buys from the cheapest plant that reaches it. `reference`
computes that sum.

```python
from nimopt.models import transport

inputs = transport.data()
model = transport.definition().build(inputs)
solution = model.solve()
print(model.n_columns, "columns for", len(inputs["cost"][1]), "arcs")
print(solution.objective, transport.reference(inputs))
```

Output:

```text
12 columns for 12 arcs
100.99601811246072 100.99601811246073
```

Arcs are drawn from every warehouse but the last, so the last warehouse is
reached by no plant and has no demand row. `absent` reports the row and the
rule that dropped it.

```python
from nimopt.models import transport

model = transport.definition().build(transport.data())
print(model.absent("demand"))
```

Output:

```text
demand  5 of 6 rows  stated by terms
  row absent  W='w5'  term-does-not-reach (flow)
```
