# /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 and is not stored. A variable over millions of columns costs
its members, not its columns.

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

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

**A definition is a model without its data.** `Definition` provides the
vocabulary a model is written in, `set`, `param`, `var`, `constraint` 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`. One definition builds a model for each dataset it is given.

**A built model is inspected through three methods.** `explain()` reports
what it built, and `row(name, **coords)` reads one row out of the assembled
matrix. `absent(name)` reports which coordinates were dropped from a
constraint and by which rule. All three read the assembled model, and none
walks the expression a second time.

**A session keeps the solver open.** `model.session()` assembles once and
keeps the solver's model. `diagnose()` queries the solved instance for the
conflicting rows of an infeasible model, or for the ray of an unbounded one.
`available()` lists the installed adapters. `capabilities(name)` reports what
one adapter supports and which capabilities it rejects together. A model with
integer columns has no duals: a mixed-integer model's duals are not its
relaxation's.

**`nimblend` is the layer below.** Its vocabulary is dimensions, labels,
entries and alignment, and it contains no optimization term. 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 of the
three has a reader above it:
`coordinates()`, `values()`, `positions_of_coordinates()` and `as_coord()`.
Do not build an index matrix either: 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`, `Piecewise`, `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 contains a
constant: `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` and keeps only the second.
A relation has no truth value, and the chained form raises instead of
dropping the first bound. 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)`
returns the correct expression 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. Building a model
that way runs 24 times slower at 25 members and 275 times slower at 400, and
the factor grows with the member count. 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` and
`primal` raise where `feasible` is False. They raise at status `unbounded`
and `unbounded_or_infeasible` whatever `feasible` reports. `bound` and `gap`
are `None` at those two statuses. `dual` raises where `status` is not
`"optimal"`. A solve stopped at a limit reports `feasible` True where the
solver found a point, with `bound` and `gap` beside it. 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 keeps the row; a term
absent along a **free** dimension removes the **row**.

**Reading a MILP's duals.** A model with integer columns has no duals.
`dual()` raises; it does not return the relaxation's duals. Read `primal`.

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

**Two operands that share no dimension.** Every binary operator combines two
dimensioned operands only where they share a dimension. The rule applies to a
coefficient multiplied by a variable and to two coefficients alike. Frames
sharing no dimension raise `ValueError`; their combination would be an outer
product. A number has no dimension and scales every entry.

**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 is passed
to 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 and reports
`('G', 'T')`.

**Bypassing the `nimblend` interface.** A test fails on an import from a
`nimblend` submodule. It also fails on a read of an array's `.index` or
`.data` or a domain's `.codes`, and on a module of the package that assembles
an index matrix of its own.

## Every error, and where it is shown

The prose above covers the common mistakes. The table lists every error the
documentation demonstrates. Each row is executed to produce the message
beside it.

| Raises | Message | Shown at |
| --- | --- | --- |
| `ValueError` | the upper bound 'cap' has no value at member ('b',) of variable 'x'; give the bound a value at every member of the variable | [/guides/bounds-from-parameters](/guides/bounds-from-parameters) |
| `ValueError` | variable 'x' is declared over ('G',) and is not over ['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 is zero at 1 coordinate(s), first at {'G': 'base', 'T': 1}; remove the zeros or divide by another parameter | [/guides/coefficient-arithmetic](/guides/coefficient-arithmetic) |
| `ValueError` | frames ('G',) and ('T',) share no dimension; pass operands that share a dimension | [/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' is given over= and where= together; pass one of them | [/guides/conditions](/guides/conditions) |
| `ValueError` | variable 'x' is read at member 't9' of dimension 'T'; read it at a member that set contains | [/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'] and takes the set, not a lag of it; write the lag at the variable's reference | [/guides/lags](/guides/lags) [/reference/expression](/reference/expression) |
| `ValueError` | parameter 'rate' is read at a lag ['T']; write the lag at the variable's reference | [/guides/lags](/guides/lags) |
| `ValueError` | piecewise 'fuel' has points that are not convex, required by sign '>=' at {'G': 'a'}; use method 'incremental' | [/guides/piecewise](/guides/piecewise) |
| `ValueError` | 'max(gen[G, T]) <= 10': the syntax supports one call; write Sum | [/guides/saving-and-loading](/guides/saving-and-loading) |
| `ValueError` | member '2030-01-01T00:30' does not convert exactly to datetime64[h] at dimension 'T' of variable 'gen'; write a member in the unit of that dimension | [/guides/saving-and-loading](/guides/saving-and-loading) |
| `ValueError` | capital does not fall from base to what follows it; pass a capital cost that falls across the merit order | [/models/expansion](/models/expansion) |
| `ValueError` | absence is 'unknown' and the array has no value at 3 of 4 coordinates; pass fill=<value> to to_dense() | [/nimblend/arrays](/nimblend/arrays) |
| `ValueError` | frames ('P',) and ('Q',) share no dimension; pass operands that share a dimension | [/nimblend/arrays](/nimblend/arrays) |
| `ValueError` | label column 't' has length 2 and the value column has length 1; pass columns of equal length | [/nimblend/arrays](/nimblend/arrays) |
| `ValueError` | 3 member(s) numbered from 4 end at position 6, and dimension 'k' has extent 6; pass a smaller start or a larger coord | [/nimblend/domains](/nimblend/domains) |
| `ValueError` | a domain of 3 member(s) requires values of shape (3,); 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']; add an entry for each | [/reference/definition](/reference/definition) |
| `ValueError` | parameter 'S' is already declared as a set; declare another name | [/reference/definition](/reference/definition) |
| `TypeError` | a relation already has one bound; compare the expression again in its own constraint | [/reference/expression](/reference/expression) |
| `TypeError` | a relation has no truth value; write each bound in its own constraint | [/reference/expression](/reference/expression) [/tutorial/constraints](/tutorial/constraints) |
| `TypeError` | an LP has no row for a strict inequality; write `<=` or `>=`, and reduce with `Sum` in place of `min` or `max` | [/reference/expression](/reference/expression) |
| `TypeError` | an expression has no absolute value: expressions are linear; bound the expression with two rows, or reduce it with `Sum` over its sets | [/reference/expression](/reference/expression) |
| `TypeError` | an expression is reduced over the sets it is summed across; specify them with `Sum(I, J, expression)` | [/reference/expression](/reference/expression) |
| `TypeError` | cannot divide by an expression: expressions are linear; declare the reciprocal as a coefficient the variable multiplies | [/reference/expression](/reference/expression) |
| `TypeError` | cannot raise an expression to a power: expressions are linear; raise a coefficient to the power and multiply it by a variable | [/reference/expression](/reference/expression) |
| `ValueError` | term 'x' already sums over ['T']; sum over each dimension once | [/reference/expression](/reference/expression) |
| `ValueError` | constraint 'cap' gives where= a domain with no name; declare its members as a parameter and refer to that parameter | [/reference/files](/reference/files) |
| `ValueError` | parameter 'c' is given columns ['value', 'S']; a table lists the dimensions then value: ['S', 'value'] | [/reference/files](/reference/files) |
| `ValueError` | variable 'x' contains the unknown key 'bound'; write only 'sets', 'subset', 'lower', 'upper', 'integer' | [/reference/files](/reference/files) |
| `ValueError` | constraint 'cap' has no row at {'P': 'p3'}; read `absent('cap')` for 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; pass columns of equal length | [/reference/param](/reference/param) |
| `TypeError` | parameter 'price' is over ('G',) and expresses no coefficient until it is read; read it at its sets as price[G] | [/reference/param](/reference/param) |
| `ValueError` | status is 'infeasible' and the solver reports no feasible point; read `status` before reading values | [/reference/solution](/reference/solution) [/tutorial/solving](/tutorial/solving) |
| `ValueError` | model 'm' has integer columns and 'highs' reports no duals for it; read primal values only | [/reference/solvers](/reference/solvers) |
| `TypeError` | parameter 'supply' is over ('P',) and expresses no coefficient until it is read; read it at its sets as supply[P] | [/tutorial/constraints](/tutorial/constraints) |
| `ValueError` | absence is 'unknown' and the array has no value at 1 of 6 coordinates; pass fill=<value> to to_dense() | [/tutorial/reading-the-answer](/tutorial/reading-the-answer) |
| `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, in the form of parameters, variables and constraints. The declaration is expanded into a coefficient matrix at assembly or at solve. Solutions are returned as arrays over the same index sets. A primal value is read by label, not by column position.

`nimopt` is built on `nimblend`, a labeled sparse N-dimensional array library. `nimblend` contains no optimization vocabulary and does not import `nimopt`. It 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: 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 modeled. A missing result is never counted as zero. Division by an absent value raises an error; it does not return an infinity.

**A subset determines the columns.** 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 after it.

**Expressions are symbolic.** An expression contains references to variables and parameters, not their values. `cost[P, W] * x[P, W]` is the same expression over a million routes and 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; it is not written incompletely. `absent()` lists every dropped row with the rule that dropped it. `row()` returns one row of the assembled matrix in the form passed to the solver.

## Install

```bash
pip install "nimopt[highs]"
```

The extra installs `nimopt` and its dependency `nimblend` from PyPI, with HiGHS as the solver backend.

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. `capabilities(name)` reports what one adapter supports, whether or not its backend is installed.

Development installs come from a checkout. `nimblend` is a dependency and installs first, from its clone. `nimopt` then installs from its own root:

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

## 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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` declares the same model without binding data. Its sets and parameters are declared by name and its constraints use the same expression syntax. `explain()` reports the whole declaration before any value is 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
d.constraint("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. One definition builds a model for each dataset it is given, and no build changes the definition. The built model is inspected the same way. `row()` reads one row out of the assembled matrix in the form passed to the solver.

```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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
d.constraint("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 instead of a full product, it has one 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()` reporting what each adapter supports and which capabilities it rejects together
- One option vocabulary translated into each solver's own option names, with a time limit written the same way for every solver
- 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 with its formulation, its inputs at any size, and an objective computed by arithmetic instead of by a solver

## Performance

Where a variable's columns are a subset of a set product, not materialising the product saves memory and time. On a transport model of 400 000 arcs over a 20 000 000-cell product, `nimopt` builds the matrix in 72.9 MB of resident memory against linopy's 1 682.9 MB. The build takes 292.9 ms against 844.5 ms.

Where nothing is sparse, the alignment work costs time and saves nothing. On a fully dense temporally coupled model at 2 111 080 rows, the comparison reverses: linopy builds the matrix three times faster, for seven percent more resident memory.

The benchmark suite measures both models. 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.

## License

MIT. See `LICENSE`.

## Citing

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

> Gaete-Morales, Carlos. *nimopt* (version 0.4.0). 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
minimize    Σ_{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 and the parameters contain the
data. `m.var` declares the decision variable, `m.constraint` 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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. The second block
therefore repeats the model. Each block runs in a Python session as it is, or
opens 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, and so 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.constraint`.

**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 is indexed over exactly
the frame of the constraint. Any other frame raises `ValueError`. Each row
then has its own value.

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

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

**Materialise.** Evaluate a parameter or an expression into an array of
values. Materialisation runs 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 modeled, used by solutions.

**Session.** A solver instance kept open on one assembled model. A caller
queries it after the solve, for the conflicting rows of an infeasible
model.

**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.constraint`. Rows over an expression's frame, bounded by a
right-hand side.

```
Model.constraint(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 contain |
| `row_of(name)` on the `Assembled` | the position of those rows 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 keeps the row. A term absent along a free dimension removes the row: a
row missing one of its terms would express a constraint that was not
written.

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

The expression is symbolic, and the constraint stores the term list and no
block. The expression is materialised once to compute its shape and once to
write it, and it stores nothing between the two.

```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.constraint("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 free dimensions of the constraint. A parameter gives each row its
own value. A parameter over other dimensions raises `ValueError`, and the
message gives both sets of dimensions.

```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.constraint("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 contains
no data. A set declared here identifies a dimension and has no members, and a
parameter identifies a shape and has no values.

An expression contains references, not arrays. The free dimensions of an
equation and its sense are read from the relation, and neither is 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` |
| `constraint(name, relation, where=None, over=None)` | nothing; registers the constraint |
| `piecewise(name, x, x_points, y, y_points, sign, method, active=None, relaxed=False, where=None)` | a `Piecewise`; `build` generates its declarations and checks its breakpoints |
| `build(data)` | a `Model` over the declarations, bound to `data` |
| `explain()` | an `Explanation` of what is declared |
| `to_yaml(instructions=False, version=4)` | the text of this definition's file, structure and no data; `instructions=True` adds the comment block that describes the format; `version=3` raises `ValueError` for a definition with a piecewise declaration |
| `set_objective(expression)` | nothing; sets the objective |
| `sense` | `"min"` or `"max"`, as declared |
| `sets`, `aliases`, `parameters`, `variables`, `constraints`, `piecewise_declarations` | 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.constraint("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. The data a definition is built from
is keyed by declared name, and one key identifies one symbol. Declaring a
parameter under the name of a set 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; declare another name
```

Equations are in no data mapping. A constraint may therefore 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.constraint("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 sets of the
definition. A model relates a set to itself through an alias. The alias has
no data of its own and reads the labels bound to its base set. `build` takes
members for the set and none for the alias, and an alias in `data` raises
`ValueError`.

```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.constraint("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 the coordinate of each set, and a
declared set has none. `where=` and `over=` on `constraint`, and `subset=` on
`var`, therefore take a tuple of the sets of the definition, meaning their
full product. They also take one of its parameters, whose coefficients are
the coordinates. Both forms resolve to the same domain. 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 the name of a declared set to its
members and the name of a declared parameter to its values. The definition is
unchanged, and it builds one model per dataset it is given.

The values of a parameter are given dense over its product, as an array of
one value per cell. They are also given long over its entries, as one mapping
of label columns and one value column. The long form gives a parameter with
coefficients at some coordinates and none at the rest. A variable declared
with `subset=` that parameter takes its members from those coordinates.

```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.constraint("supply", Sum(W, cost[P, W] * flow[P, W]) <= supply[P])
d.constraint("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 omits a declaration, or contains a key 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']; add an entry for each
```

---

# /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. A reader takes a field and parses no rendered text.

| Field | Contains |
| --- | --- |
| `name`, `sense` | the name of the declaration and the direction it optimizes |
| `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 |
| `piecewise` | one `PiecewiseShape` per piecewise declaration |
| `objective` | the objective expression, 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 value a caller acts on, and a declaration with no data reports no
count.

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

`VariableShape.members` identifies the parameter a sparse variable took its
members from, and is `None` for one over the full product. Columns are absent
until data binds. Without that field a sparse declaration and a dense one
read identically.

`ConstraintShape.free` and `.sense` are read off the relation. Neither is
declared beside it. An expression contains 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.constraint("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])
```

`PiecewiseShape.generated` lists the variables and constraints a model
generated for the declaration. A definition generates none, and the tuple is
empty.

---

# /reference/expression

# Expressions

## `Term`

One variable, an optional coefficient, the dimensions summed over, and a
scale factor. A term describes a block of coefficients and contains
references, not arrays. Writing it allocates 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 contains |
| `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 the entries of each term before
the reduction. The sum then runs over the coordinates given, not over 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 is over the members of a set and takes the set, not a lag of it. Write
the lag at 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'] and takes the set, not a lag of it; write 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.constraint` 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` and keeps only the second. The chained form
raises, and the first bound is not 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; write each bound in its own constraint
```

## 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: cannot raise an expression to a power: expressions are linear; raise a coefficient to the power and multiply it by a variable
```

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: cannot divide by an expression: expressions are linear; declare 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: an expression has no absolute value: expressions are linear; bound the expression with two rows, or reduce it with `Sum` over its sets
```

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; write `<=` or `>=`, and reduce with `Sum` in place of `min` or `max`
```

The built-in `sum` of expressions with no set to reduce over calls
`Expression.sum`. That call raises, and it returns no unreduced expression.

```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; specify them with `Sum(I, J, expression)`
```

A relation expresses one bound. Comparing it a second time raises, and the
first bound is not dropped.

```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 already has one bound; compare the expression again in its own constraint
```

A dimension is reduced once. A second reduction over the same dimension
raises.

```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']; sum over each dimension once
```

---

# /reference/files

# Files

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

| Function | Does |
| --- | --- |
| `load(path, data=None)` | reads a file; returns a `Definition`, or a `Model` where the file contains data or `data=` gives it |
| `loads(text, data=None)` | the same over text; a sidecar name in text raises `ValueError`. Text has no directory |
| `save(what, path, inline=False, instructions=False, version=4)` | writes a definition's file, or a model's with an `.npz` beside it, or one file with an inline block when `inline=True`; `instructions=True` writes the comment block that describes the format at the top of the file; `version` is `4` or `3` |

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

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

`version=4` writes each piecewise declaration under `piecewise`, and omits
the sets, parameters, variables and constraints it generated. `version=3`
writes a model's generated declarations as ordinary declarations, and omits
the piecewise declarations and any parameter only they read. A model loaded
from that file contains the same rows and no piecewise declaration. A
definition with a piecewise declaration raises `ValueError` at version 3. Any
other version raises `ValueError`.

With `instructions=True`, every writer puts a fixed comment block at the top
of the text. The block describes the format: the keys and their order, the
defaults, the rules that determine which rows a constraint has, and the
expression syntax. The text is the same in every file and describes the
format, not the model. A reader with one file interprets it without this
package. The block is a YAML comment: a file with it and a file without it
load to the same model.

## The file

| Key | Contains |
| --- | --- |
| `version` | `4`, or `3` where the caller asks for it; files of version `2` and `3` also load; any other value raises `ValueError` and reports the versions this reader accepts |
| `name`, `sense` | the model's |
| `sets` | a list of names |
| `aliases` | each alias to its base set; 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 |
| `piecewise` | each name to `x`, `x_points`, `y`, `y_points`, `sign`, `method`, and to `active`, `relaxed` and `where` where given; version 4 only; absent where the model declares none |
| `objective` | the objective expression; absent where the model declares 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 contains, or a list of set names, meaning their full product. A symbol's
name is a Python identifier other than `Sum`. A dimension listed 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 its base set.

The expressions are written as they are typed in Python and are read back
through the same operators. 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.constraint("cap", 2 * c[S] * x[S] - 1 <= 5)
text = d.to_yaml()
print(text)
print(loads(text).to_yaml() == text)
```

Output:

```text
version: 4
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 define raises `ValueError`, at the top level and
inside an entry.

```python raises=ValueError
from nimopt import loads

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

Raises ValueError:

```text
ValueError: variable 'x' contains the unknown key 'bound'; write only 'sets', 'subset', 'lower', 'upper', 'integer'
```

## Data

The data in a file 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 set of `datetime64` or `timedelta64` members | `dtype` and `members` | 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 raises `ValueError` at save and reports the symbol: the
container would pickle it, and the reader rejects a pickled array.

```python raises=ValueError
from nimopt import loads

loads(
    "version: 3\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 lists the dimensions then value: ['S', 'value']
```

## Datetime members

A set whose members are `datetime64` or `timedelta64` is written inline as a
mapping of `dtype` and `members`, instead of a list.

```yaml
data:
  T:
    dtype: datetime64[s]
    members: ['2030-01-01T00:00:00', '2030-01-01T01:00:00']
```

A `datetime64` member is written as its ISO 8601 string. A `timedelta64`
member is written as its integer count of the unit in the `dtype`. A label
column of a long table is written in the same text and takes no marker: the
column belongs to a set, and the reader converts it to that set's dtype. The
`.npz` stores the dtype of every array and uses no separate form.

A member fixed in a relation is written as text. A `datetime64` member is its
quoted ISO 8601 string, such as `x['2030-01-01T00:00:00']`. A `timedelta64`
member is a quoted count and numpy unit code, such as `x['3 h']`.

Every member is converted to the dtype of its set. A string is parsed as ISO
8601. A `datetime.datetime`, a `datetime.date` and a `datetime64` of another
unit are converted. An integer against a `timedelta64` set is a count of that
set's own unit. A conversion that is not exact raises `ValueError`:
`'2030-01-01T00:30'` against a set in hours raises instead of truncating to
the hour.

A member specifies no time zone. `datetime64` represents no offset, and a
conversion to UTC would move the member. A string with an offset or a trailing
`Z` raises `ValueError`, and a `datetime.datetime` with a `tzinfo` raises the
same error. Write the naive form, `'2030-01-01T00:00:00'`.

A member outside the range its dtype represents raises `ValueError` and
reports the first and the last member of that range. A `NaT` member raises
`ValueError`.

## What raises before anything is written

| Written | Reason |
| --- | --- |
| 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 expression syntax 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.constraint("cap", x[P] <= 1.0, where=subset((P,), {"P": np.array(["a"])}))
m.to_yaml()
```

Raises ValueError:

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

---

# /reference/inspection

# Inspecting a built model

## `Row`

Returned by `Model.row(name, **coords)`. One row as the assembled matrix
stores it. The row is read from the matrix, not from a second walk of the
expression, and it shows what is passed to the solver.

| Field | Contains |
| --- | --- |
| `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 | Contains |
| --- | --- |
| `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`.
A column resolves to its variable through that range, and to a coordinate
through the numbering rule of that variable.

`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.constraint("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. `row` and
`absent` raise `KeyError` for a name that is not a declared constraint, and
the message lists the declared constraints. `row` raises `KeyError` for a
label that is not a member of its dimension, and the message identifies the
dimension.

```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.constraint("cap", one[P] * x[P] <= rhs[P])

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

Raises ValueError:

```text
ValueError: constraint 'cap' has no row at {'P': 'p3'}; read `absent('cap')` for 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 | Contains |
| --- | --- |
| `constraint` | the equation this is about |
| `stated_by` | `"terms"` where the rows are derived, `"over"` where given explicitly |
| `expected`, `standing` | rows expected, 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; the row would express a constraint that was not written |
| `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; the row is kept with one term fewer |

The two rules differ in what they remove. A coefficient absent inside a sum
removes a **term** and keeps the row. A term absent along a **free**
dimension removes the **row**.

Under `over=` the rows are given explicitly. Nothing is dropped, and a
right-hand side that omits one raises instead. An empty `dropped_rows` beside
`stated_by="over"` follows from that rule.

```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.constraint("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 contains 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 `ValueError`.

| Member | Returns |
| --- | --- |
| `var(name, sets, subset=None, lower=0.0, upper=inf, integer=False)` | a `Variable` occupying the next range of columns |
| `constraint(name, relation, where=None, over=None)` | a `Constraint` occupying the next range of rows |
| `piecewise(name, x, x_points, y, y_points, sign, method, active=None, relaxed=False, where=None)` | a `Piecewise`; declares the variables and constraints of its method |
| `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, instructions=False, version=4)` | the text of this model's file, with its data inline where asked and the comment block that describes the format where asked; `version=3` writes the declarations a piecewise declaration generated in its place |
| `piecewise_declarations` | the piecewise declarations, keyed by name |
| `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.constraint("supply", Sum(W, x[P, W]) <= 30.0)
m.constraint("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.]
```

## `Piecewise`

```
Model.piecewise(name, x, x_points, y, y_points, sign, method, active=None, relaxed=False, where=None)
Definition.piecewise(name, x, x_points, y, y_points, sign, method, active=None, relaxed=False, where=None)
```

A piecewise-linear relation of the expression `y` to the expression `x`.
`x` is on the curve through `x_points` and `y_points`. `sign` compares `y`
with the curve: `"=="`, `"<="` or `">="`. The two points are parameters read
at their sets. Each is over some or all of the sets of `x` and over one
breakpoint set, the one set `x` is not over. An entity lists its first
breakpoints, and its last breakpoints may be absent. An entity with no
breakpoint has no generated rows and no generated columns.

`where` restricts the declaration to some entities: a parameter, a tuple of
sets or a domain over the sets of `x_points` other than the breakpoint set.
The breakpoint checks, the generated columns and the generated rows cover
the entities at its coordinates. `x`, `y` and `active` are compared at those
coordinates only.

| `method` | Generates | Requires |
| --- | --- | --- |
| `"incremental"` | per segment, one continuous and one integer column and their rows | breakpoints strictly increasing or strictly decreasing |
| `"tangent"` | one row per segment, and two rows that keep `x` between the first and the last breakpoint | points convex under `>=`, concave under `<=`; no `active`; no `==`; no constant in `x` |

`active` is a binary variable over the sets of `x`, or a sum of them. Where
it is 0, `x` is 0 and `y` is compared with 0. A term that is scaled or
bounded outside 0 and 1 raises ValueError, and so does a continuous term
under the default. `relaxed=True` accepts a continuous `active` between 0
and 1 and scales the curve by its value, which is the linear relaxation of
the switch. `relaxed=True` with no `active` raises ValueError. `Model.piecewise` generates the declarations at
once. `Definition.piecewise` stores the declaration, and `build` generates
them. A generated name is `name`, an underscore and a suffix:

| `method` | Sets | Parameters | Variables | Constraints |
| --- | --- | --- | --- | --- |
| `"incremental"` | `segment` | `members`, `x_step`, `y_step`, `x_first`, `y_first` | `fill`, `order` | `x`, `y`, `order_bound`, `fill_order`, `order_link`, `active` |
| `"tangent"` | `segment` | `slope`, `intercept`, `x_low`, `x_high` | none | `tangent`, `x_min`, `x_max` |

`{name}_active` exists only where `active` is given. The members of
`{name}_segment` are the breakpoint set's members without the first. A
segment is identified by its end breakpoint.

| Member | Contains |
| --- | --- |
| `name`, `x`, `x_points`, `y`, `y_points`, `sign`, `method`, `active`, `relaxed`, `where` | the arguments |
| `breakpoints` | the name of the breakpoint set |
| `names()` | the names the declaration generates, keyed by `"sets"`, `"parameters"`, `"variables"` and `"constraints"` |
| `generated` | the names a model generated, keyed the same way; empty on a definition |
| `generated_names()` | every generated name, as a frozenset |

An argument error raises when the declaration is made. `TypeError` is
raised for an `x`, `y` or `active` that is not an expression, and for points
that are not a parameter read at its sets. `ValueError` is raised for a name
that is not a Python identifier, an unknown `method` or `sign`, expressions
over different sets, points without exactly one breakpoint set, points over
different sets, `"tangent"` with `"=="`, with `active` or with a constant in
`x`, an `active` with a constant, a `where` of another type or over other sets
than the entity sets of `x_points`, and a generated name the model or
definition declares.

A breakpoint error raises `ValueError` when the data is bound, before any
declaration, and identifies the first entity at fault: points with no
breakpoint, points present at different breakpoints, an entity with one
breakpoint, an absent breakpoint before a present one, a value that is not
finite, breakpoints that are not strictly monotonic, and, for `"tangent"`,
points whose curvature does not match `sign`.

## `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.constraint("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 record type a `Definition` returns with every
count absent, and one reader covers both.

A model contains variables and constraints. Its sets and parameters are
collected from them, in order of first appearance. A dimension introduced by
a coefficient 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.constraint("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.

The array of a parameter declares `absence="empty"`. A coordinate it does
not have contributes no coefficient. Absence is the additive identity a sum
requires.

### `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 and 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; pass columns of equal length
```

### Members

| Member | Returns |
| --- | --- |
| `dims` | 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 a coefficient, and so is an arithmetic
combination of coefficients. One interface therefore covers both.

`+`, `-`, `*`, `/` and a power by a number combine coefficients. The
combination is symbolic: it contains 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' is over ('G',) and expresses 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 instead of dropping |
| `set - 1` | the set lagged, dropping the members outside the set |

```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. A model
relates a set to itself 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 declares
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`
contains one label column per set, keyed by the name of the set. The columns
are read in parallel: the k-th entry of each column belongs to the same
member. The result 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
with positions passes them directly and builds no 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`. 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 |
| `feasible` | whether the solver reports a primal-feasible point |
| `objective` | the objective value of that point |
| `bound` | the bound on the optimal objective the solver proved, or `None` |
| `gap` | the relative distance from the objective to the bound, or `None` |
| `primal(name)` | the named variable's values over its own sets |
| `dual(name, kind=None)` | a constraint's duals over its free sets, or a variable's reduced costs over its own sets |

`status` and `feasible` are readable whatever the solver reported.
`objective` and `primal` raise `ValueError` where `feasible` is False. They
raise at status `unbounded` and `unbounded_or_infeasible` whatever `feasible`
reports. An unbounded model has no optimal value, and `bound` and `gap` are
`None` there.
A solve stopped at a limit reports `feasible` True where the solver found a
point, and those reads then return it. `dual` raises `ValueError` where
`status` is not `optimal`. Read `status` first.

`primal` and `dual` raise `KeyError` for a name the model does not declare.
`primal` takes a variable, and its message reports a constraint name as one
`dual` reads. `dual` takes either, and its message lists the declared
constraints and variables. Both raise `KeyError` for a name declared after
the solve.

A model declares its constraints and its variables in two registries, so one
name identifies one of each. `dual` raises `ValueError` for such a name and
reads it under `kind="constraint"` or `kind="variable"`. Any other `kind`
raises `ValueError`.

`dual` returns a reduced cost for a variable: its objective coefficient less
the duals of the rows it appears in, weighted by its coefficients in them,
in the model's own objective under either sense. nimopt derives the value
from the row duals the solver reports, so the convention does not vary by
solver. The values follow the dual solution the solver returns. A degenerate
model has more than one such solution, and two solvers can report different
reduced costs for it. A reduced cost follows the
variable's members by the rule `primal` follows: a `DenseArray` over a full
product, a `SparseArray` over a subset.

`bound` is a lower bound on the optimal objective under sense `min` and an
upper bound under sense `max`. It is `None` where the solver reports none.
For a model without integer columns it is the objective at status `optimal`
and `None` at any other status. `bound` is readable at every status, and the
solvers report none at status `unbounded` and `unbounded_or_infeasible`.
`gap` is `abs(objective - bound) / abs(objective)`. It is `None` where
`feasible` is False, where `bound` is `None`, and at status `unbounded` and
`unbounded_or_infeasible`.

```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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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 where the solver reports no feasible point 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.constraint("floor", Sum(W, x[P, W]) >= 10.0)
m.constraint("ceiling", Sum(W, x[P, W]) <= 1.0)
m.set_objective(Sum(P, W, x[P, W]))

m.solve().objective
```

Raises ValueError:

```text
ValueError: status is 'infeasible' and the solver reports no feasible point; 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. The
solver returns those values in column order, and they reshape into a
`DenseArray` with no index built. A variable over a subset has values at its
members alone. A dense frame would be the grid the declaration avoids, and
those values remain a `SparseArray`. A dual follows the rows of its
constraint by the same rule.

Every array declares `absence="unknown"`. A coordinate the model does not
have has no value, and combining the results of two models adds no zero for
it.

---

# /reference/solvers

# Solvers

## `available` and `capabilities`

`available()` lists every adapter whose backend can be imported in the
current environment, with the capabilities each declares. `capabilities(name)`
reports for an adapter whether or not its backend is installed. A descriptor
describes the adapter as shipped. A caller reads one to choose 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  rejects duals+integrality
gurobi  integrality native · duals native · conflict native · ray native  rejects duals+integrality
mosek  integrality native · duals native · conflict absent · ray native  rejects duals+integrality
```

## `Capabilities`

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

The capabilities are `integrality`, `duals`, `conflict` and `ray`. A flat
set is insufficient: a solver can support two and reject their
combination. Every adapter rejects `integrality` with `duals`. A
mixed-integer model's duals are not the relaxation's. 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.constraint("cap", one[T] * x[T] <= 2.0)
m.set_objective(Sum(T, one[T] * x[T]))

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

Raises ValueError:

```text
ValueError: model 'm' has integer columns and 'highs' reports no duals for it; read primal values only
```

## `Session`

Returned by `Model.session(solver="highs", options=None)`. The model of one
solver, opened on one assembled model and kept open. A solve passes the
matrix to the solver, and a later query reads 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. A caller who needs
only a solution needs no session.

The session assembles the matrix when it opens. `solve()` and `diagnose()`
raise `ValueError` where the model declares a variable or a constraint after
that. Open a new session on the changed model.

```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.constraint("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` contains the same `Row` that `Model.row` returns, and a
conflicting row reads in one format. Two solvers may return different
irreducible sets. Removing either set makes the model feasible. The two
sets need not be equal.

A conflict is computed on the solved backend, and the session keeps that
backend 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.constraint("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. It reports no row it
did not prove. Gurobi's conflict covers the integrality. Mosek's adapter
computes no conflict. A session on it raises, and the message refers to
`capabilities("mosek")`.

## `options` and `Option`

`options()` lists every option a caller can set, under the names of
`nimopt`. An option outside the list raises `ValueError`. No option is passed
to a solver that ignores it, and an unrecognized name stops the solve.

`options(solver)` lists the same options with the name and the values of that
solver. A caller follows those names into the documentation of the solver.

```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 randomizes 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 factorizes | `hipo_system` | not supported | not supported |
| `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 supported | not supported |

A choice each solver writes differently is given once and translated. The
value a caller writes has one meaning for every solver. No solver supports
every option or every choice. `newton_system` and `pdlp_tol` are specific to
HiGHS, and so are `hipo` and `pdlp` under `method`. Passing one of them to
Gurobi or Mosek raises, and the message identifies it. No solver runs a
different algorithm in its place. Mosek runs only its mixed-integer
optimizer on a model with integer columns. `method` is `choose` there, and
any other choice raises. The guide on [interior point and first-order
methods](/guides/highs-methods) gives the memory each method requires. It
also gives the installation of a HiGHS with HiPO and a GPU.

## 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 given by the caller is used as passed, and a notebook
or another interface writes through it. A reporter implements
`start(total, what)`, `step(done, what)` and `done()`. That is the whole
contract.

The report covers building. Its resolution is the structure of the model. The
measuring pass counts constraints and the writing pass counts nonzeros. A
model with one constraint of one term reports one step and no fraction.

`log=True` requests the log of the solver. Building finishes before a solver
starts, and the report and the log do not 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 the codes of a subset. Nothing stores a column index. A
variable over millions of columns therefore 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 single column space of the
model. `m.n_columns` counts every column declared so far.

## A variable over no dimension

The bracket of a variable lists the dimensions it is declared over. A
variable over no dimension takes no bracket and enters a row on its own. It
is one column: a value-at-risk level, a budget slack, or a bound shared by
every row of a family. `theta[()]` is the same term written out.

A variable over one or more dimensions expresses no term until it is read.
Using one without a bracket raises `TypeError` and reports the reading it
requires. The same rule applies to a parameter, read as `cost[G, T]`, and as
`k` over no dimension.

Comparing a variable expresses a row, and `==` between two variables
expresses a row as well. A list of variables therefore cannot be searched
with `in` or `.index`. Both compare their items and raise `TypeError`
and report the first variable they compare. A dict and a set match on
identity. Store variables in one of them, 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.constraint("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 written 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 correspondence between
a model and its matrix: the column space is a dimension, and 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. Both appear when a `nimblend` array from
inside a model is inspected.

---

# /nimblend/arrays

# nimblend arrays

`nimblend` is the layer below `nimopt`. Its vocabulary is dimensions, labels,
entries and alignment, and it contains no optimization term. A `nimopt` caller
uses these names when reading a solution or inspecting what a model built.

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

Do not assemble one either. `SparseArray.from_canonical` below takes an index
matrix built by the caller, and that is the work of the array layer. No module
of `nimopt` calls it, and a test enforces 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 the interface a caller writes against.

`absence` declares the meaning of a coordinate the array does not have:
`"empty"` that it contributes nothing, `"unknown"` that it was not modeled.
Operators and reductions follow from that declaration. The declaration is
part of the definition of the array, not a hint.

| Member | Returns |
| --- | --- |
| `dims`, `shape`, `nnz` | the dimensions, their sizes, and the number of entries |
| `coords` | the coordinate of each dimension, resolving 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`: a product intersects presence and has at most the entries of
the sparse operand.

Absence and zero remain distinct. A stored `0.0` is a coordinate present with
the value zero. A coordinate the array does not have is a different case.

```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, and this array has four entries
and not two. The zeros are stored, and a stored value is 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.

The storage of presence follows the absence declaration: the two
declarations require opposite behavior from an operator. An `"unknown"` array
marks absence with NaN. NaN propagates through arithmetic at no cost and
requires no storage beside the values. An `"empty"` array stores a boolean
mask: absence is the additive identity there, and substituting it costs less
than marking it.

`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, and
the values reshape with no index built.

```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. A caller with coordinates of its own therefore gives its entries
as labels and resolves no position itself. The columns are read in parallel
and 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; pass columns of equal length
```

`is_canonical` reports whether buffers meet the precondition of
`SparseArray.from_canonical`: entries sorted by ravel key, with no repeat.
Verifying it inside `from_canonical` would cost the ravel that this path
avoids.

```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 dimensions of the two operands alone. A caller reads it
before materialising either operand. A combination therefore reports its
frame while its data is 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 no common dimension to align on, and
`combined_dims` raises. Their combination would be an outer product.

```python raises=ValueError
import nimblend

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

Raises ValueError:

```text
ValueError: frames ('P',) and ('Q',) share no dimension; pass operands that share a dimension
```

## 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: absence is 'unknown' and the array has no value at 3 of 4 coordinates; pass fill=<value> to to_dense()
```

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, and a domain resolves them
for a tuple of dimensions. It reports the multi-indices it has, the position
of each, and the multi-index at a given position. It records which
coordinates are present, and it records no numbering origin.

| 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. An array's `.index` and `.data` are its raw
buffers. Never read them.** `coordinates()` and `labels()` report which
members are present, and `positions_of_coordinates` reports the position of
one. `as_coord` numbers them, and `array` and `identity` return an array over
them. No layer above builds an index matrix. 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, and the query returns `-1`. The other
two are the first and second members of the domain, and `as_coord(100)`
numbers them from 100.

## Crossing a domain with further dimensions

`expand` replicates every member across the full extent of the given
dimensions. That enumerates the coordinates a term can have, before the test
of which ones it does have. The new dimensions are appended, and `transpose`
reads the result in another order. The code of a member is its own code
scaled by the appended extent, plus each position within it. The cross
product is therefore arithmetic on the members, and it builds no index
matrix.

```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 members, and `transpose`
presents them over the dimensions in another order. The members present do
not change.

## A domain returns an array

A caller with one value per member calls the domain, and so does a caller
that pairs each member with its own position along a new dimension. Neither
builds an index matrix. Both calls read above the raw members, as
`coordinates()` and `as_coord()` do.

`array(values)` assigns one value to each member, in their stored order. The
members ascend, 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.]
```

The members of a full domain ascend with the ravel key. That is the order
`values.ravel()` reads a grid in, and 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) requires values of shape (3,); got shape (2,)
```

`identity(into, coord, start)` pairs each member with its own position along
a new dimension, valued 1.0. The position of a member is its rank plus
`start`, the numbering `as_coord(start)` uses. An array built one way and a
coordinate built the other place a member at the same position. `coord` is
the coordinate of the new dimension. It spans the whole extent the positions
are numbered into. That extent is 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 to be numbered raises `ValueError`,
and it writes no position outside itself.

```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 end at position 6, and dimension 'k' has extent 6; pass a smaller start or a larger coord
```

## The three coordinates

A coordinate resolves the position of a label along one dimension. The
dimension determines which of the three is used.

| Coordinate | Is | Used for |
| --- | --- | --- |
| `StoredCoord(labels)` | labels stored as an array | a dimension whose members have labels |
| `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 the position of an entry as its rank among the codes. A
block already in canonical order needs no lookup.

```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 labels it stores, and `"c"`
resolves to position 2. `ProductCoord` ravels a multi-index against the
sizes, and `(0, 2)` is position 2 and `(1, 0)` is position 3. `SubsetCoord`
stores the codes `0` and `4`, the 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 that returns successive slices.

A block computed into a reserved slice never exists as a separate object.
Assembling several blocks therefore stores one copy of the result, and not
one copy per block plus the result. The buffer is the destination the
assembly of a model 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.]]
```

`reserve` returns views of the one allocation. A write into a slice writes
into the returned array.

---

# /nimblend

# nimblend

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

`nimopt` imports it, never the reverse. A model uses `nimblend` in two places:
a solution is returned as a `nimblend` array, and the rows of a constraint are
a `nimblend` domain. `nimblend` is also usable on its own, for labeled sparse
data outside a model.

## 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 modeled.
Division by an absent value raises an error and returns no 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 no storage: `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 the
position of each. Through `array` and `identity` it returns an array over
those members, and a caller assembles no 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` stores the codes in order, and a block already
in canonical order needs no lookup.

```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 stores its
set sizes, the start of its block and, for a subset, the codes of its
members. It stores no integer per member: the column is computed 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 report the position a
member occupies. They differ in whether that position is computed by
arithmetic or by a rank.

---

# /explanation/expressions-are-symbolic

# Expressions are symbolic

A `Term` describes a block of numbers and stores none. It contains 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. That 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
block of the variable, and 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 stores the term list and no block. Adding a constraint
therefore costs the computation of its shape. `n_rows` and `nnz` are known at
declaration, and no coefficient exists 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.constraint(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 stores 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. One expression is live at a
time, and peak memory is set by the largest constraint, not by the sum of
all of them.

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

---

# /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: the array is the matrix.

```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.constraint("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__')`. That is a matrix: 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. The block never exists as a second object.

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

Passing the matrix to a solver returns views. `indices` and `values` are the
buffer, and only `indptr` is built. A model of four million nonzeros passes
its matrix without copying it.

## Why the shape is computed first

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

Where the data of a parameter changes between the two, the constraint
computes one number of coefficients and writes another. The model raises, and
it writes no matrix that differs from the shape it reported.

## What the design costs

Materialising the expression twice costs build time. One expression exists at
a time, and the blocks of the constraints never exist together. Peak memory
is set by the largest constraint, not by the sum of all of them.

The declared shape of a model is also exact before anything is built.
`n_rows`, `n_columns` and `nnz` are exact counts, 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 contains no optimization term.
A function in it referring to a row, a column or a constraint is 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 its matrix is stored in
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 the `nimblend` source
fails on any mention. A second test imports `nimblend` alone and fails where
`nimopt` is imported with it. `nimblend` enforces 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 mentions
one in prose. A package developed on its own fails in its own suite, and not
in the suite of a consumer.

**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 constrained further: the set of
`nimblend` names they import is pinned, and widening it requires an explicit
change to the test.

**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, not by a text search. It detects a read that is not
a subscript, and it permits `dims.index(name)`, a position lookup on a
tuple.

**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 performed one layer too high. That work is part of the interface that
moves 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 are allocated

The strongest boundary test measures allocation. It does not compare the
totals of the two packages.

The cost per nonzero of a build is allocated in `nimblend`: the matrix is an
int32 column and a float64 value per entry, stored in a `nimblend` buffer. A
ninefold increase in the nonzeros grows the `nimblend` share by at least
twelve bytes per added entry.

`nimopt` allocates a near-constant amount 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. The
sets fix those counts, and the density does not.

The test therefore measures the gradient and not the totals. An assertion
that `nimblend` allocates more bytes would measure the workload instead of
the boundary. `nimopt` allocates a vector per column at every density, and a
sparse enough model puts `nimopt` above `nimblend` with no boundary
violation.

## Why the site lives in `nimopt`

`nimopt` imports `nimblend`. A site inside `nimopt` documenting both packages
follows that dependency. A site inside `nimblend` documenting `nimopt` would
invert it, and the `nimblend` tests would then require a model to describe.

---

# /explanation/what-the-numbers-measure

# What the numbers measure

Every figure here is given against a declared baseline. A number without its
denominator is not interpretable, and the choice of denominator changes the
size of a ratio.

## Which package allocates 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.
Rows and columns are fixed at 200 and 20 000, and the density of the
coefficient determines the nonzeros.

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

`nimblend` allocates 16.0 bytes per added nonzero, an int32 column and a
float64 value. The matrix of a model is stored in a `nimblend.EntryBuffer`,
and the CSR conversion returns views of that buffer. `nimopt` moves by
48 bytes across a ninefold change in the matrix. It allocates the per-column
vectors, the bounds and the integrality, and the row bounds. The sets fix
those counts.

**What it does not claim.** It does not claim that `nimblend` allocates more
bytes than `nimopt`. The totals follow from the nonzeros per column of the
model. A model with one entry per column puts `nimopt` above `nimblend` with
no boundary violation. The measure is the gradient, not the totals.

## Two models

`benchmarks/bench_transport.py` ships from plants to warehouses over a
network. Each plant serves a band of nearby warehouses, and the flow variable
spans the arcs, not 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 adjacent
hours through a cyclic `state_of_charge` row, 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 models solve through HiGHS. The two lag rules are visible in the row
counts. Over 168 hours the ramp block has 10 x 167 rows: the first hour has
no predecessor, and its row is not produced. The cyclic `state_of_charge` row
wraps to the last hour and keeps all 2 x 168 of its rows.

## Peak against the matrix

The peak occurs while the model is built, not while it is assembled. On the
400 000-column transport model, the build peaks at 69.72 MB and the assembly
that follows peaks at 68.23 MB. The expression of every constraint is built
once to measure its shape, and the model that produces the matrix remains
live while the matrix is written.

The denominator determines the size of the ratio. 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 constant across the rungs, 8.09x at 4 000 nonzeros and
7.26x at 800 000. The figure of the smallest rung is stable and is not a
first-call artifact: three consecutive measurements in one process give 0.39,
0.38 and 0.38 MB. The excess of the small rung over the large one is the
fixed structures of the model. Those structures 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 occupies. An alignment
therefore costs more than the result it produces. That cost sets the lower
bound of the ratio. It is per operation and transient, and it is not
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
constraints peak at 75.85 MB against a 16.00 MB buffer. That is 2.14 MB
beyond the rows the second constraint adds. A constraint stores its term list
and no block. An expression exists while its shape is measured and again
while it is written, and never between.

**What it does not claim.** It does not claim that peak memory is 7x the
matrix in an absolute sense. It is 7.26x *the CSR matrix* and 3.31x *the full
LP data* on this model. Which of the two applies depends on the comparison
being made.

## 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 passed to both. The measured build runs from an empty
model to the matrix a solver is given: `assemble()` for nimopt, `m.matrices`
for linopy. Each side runs in a process of its own: resident memory includes
the import cost of the library that is loaded. A row is printed once the two
agree on rows, columns, nonzeros and the solved objective.

Resident memory is sampled, not traced. `tracemalloc` records only what
passes through the Python allocator. Two libraries that allocate by different
routes would then be compared on the route and not 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. The
design targets this case.

**Temporal coupling on a dense model.** Ramp limits and a cyclic
`state_of_charge` row 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 build 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, and the design costs time with no
corresponding saving. Every operation aligns by key, while xarray broadcasts
over dense grids and aligns by position. The model is also built twice, once
to measure the shape of each constraint and once to write it. Resident memory
is close on both sides: both store the same dense coefficient grids.

**Integrality.** Making the flow an integer column changes neither 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 solution 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 is
the 400 000-arc transport model, where the solution is unpacked onto the
masked product the build used.

The matrices differ in width. linopy returns a scipy matrix with `int64`
column indices, and the same 800 000 nonzeros occupy 12.80 MB against
9.60 MB for nimopt.

**What it does not claim.** It does not claim that nimopt is faster than
linopy. On a dense temporally coupled model at two million rows it is three
times slower, and the table reports that row with the others. The numbers
support a narrower claim. Where a model is sparse in its variables, not
materialising the grid saves memory and time. Where it is not, the alignment
work costs time and saves nothing.

---

# /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. Those positions 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`, and the message gives the expected
shape 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. 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` contains 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 give six members. `x` therefore occupies six
columns of the coefficient matrix. A column index is computed from the
positions of a member in each set, and nothing stores a column per member. 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 and 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 written 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,
the coefficients and the 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, and an objective takes that form.

## 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: it contains references to `cost` and `x`, not 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.constraint(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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
demand_rows = m.constraint("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 frame of the constraint, giving each row its own value. The
parameter is read at its sets here as it is anywhere else: `supply[P]`, not
`supply`. A parameter without a bracket raises `TypeError` and reports the
reading it requires.

Supply is indexed over `P`, and so are the supply rows. A parameter read over
any other index set raises `ValueError`, and 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.constraint("supply", Sum(W, x[P, W]) <= supply)
```

Raises TypeError:

```text
TypeError: parameter 'supply' is over ('P',) and expresses 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.constraint("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)` and discards the first relation. `nimopt`
raises `TypeError` on the chained form and drops no 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; write each bound in its own constraint
```

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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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. `status` and `feasible` are
readable after any solve. `objective` and `primal` raise `ValueError` where
`feasible` is False. They raise at status `unbounded` and
`unbounded_or_infeasible` whatever `feasible` reports. `bound` and `gap` are
`None` at those two statuses. `dual` raises `ValueError` where `status` is not
`optimal`.

Raising the demand of Berlin to 40 makes total demand 70 against total
supply 55. 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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: status is 'infeasible' and the solver reports no feasible point; 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).

## A solve stopped at a limit

An option in `options()` stops the solver early. A solver stopped at a
limit reports the best point it found, and `feasible` is True for it.
`objective` and `primal` then return that point. `bound` returns what the
solver proved about the optimum, an upper bound under sense `max` and a
lower bound under sense `min`. `gap` returns the relative distance from the
objective to that bound.

A search stopped after one node returns the point the solver found there.
The thread count of the solver determines that point. The example below
reports the properties every such point has.

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

rng = np.random.default_rng(1)
ITEM = Set("item", np.array([f"i{t}" for t in range(40)]))
BIN = Set("bin", np.array([f"b{t}" for t in range(5)]))
weight = Param.from_dense("weight", (BIN, ITEM), rng.uniform(1, 50, (5, 40)))
value = Param.from_dense("value", (ITEM,), rng.uniform(1, 100, 40))
capacity = Param.from_dense("capacity", (BIN,), np.full(5, 306.0))

m = Model("knapsack", sense="max")
x = m.var("x", (ITEM,), integer=True, upper=1.0)
m.constraint("capacity", Sum(ITEM, weight[BIN, ITEM] * x[ITEM]) <= capacity[BIN])
m.set_objective(Sum(ITEM, value[ITEM] * x[ITEM]))

solution = m.solve(options={"node_limit": 1})
print(f"status: {solution.status}")
print(f"feasible: {solution.feasible}")
print(f"the bound is above the objective: {solution.bound > solution.objective}")
print(f"the gap is positive: {solution.gap > 0.0}")
```

Output:

```text
status: solution_limit
feasible: True
the bound is above the objective: True
the gap is positive: True
```

The status identifies the limit the solver stopped at. HiGHS reports a stop
at `node_limit` as `solution_limit`. `gap` is `None` where the solver proved
no bound, and `feasible` is False where it found no point.

## 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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 solution](/tutorial/reading-the-answer).

---

# /tutorial/reading-the-answer

# Reading the solution

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. `dual` also takes a
variable, and returns its reduced costs: the cost of the variable less the
duals of the rows it appears in. A route with no shipment at the optimum
reports the amount by which its cost exceeds those duals.

```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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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. The supply constraint of
Lisbon binds, and each marginal unit is served from Porto at the cost of the
route from Porto.

## 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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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. 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"`. `to_dense()` raises
`ValueError`. It supplies no value for a coordinate the model does not
produce.

```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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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: absence is 'unknown' and the array has no value at 1 of 6 coordinates; pass fill=<value> to to_dense()
```

`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.constraint("supply", Sum(W, x[P, W]) <= supply[P])
m.constraint("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 ships nothing.

## Summary

The tutorial declared index sets and parameters, a decision variable,
expressions, two constraint families and an objective. It 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) present ten
complete formulations.

---

# /guides/at-scale

# Declaring a model at scale

Models with millions of columns are declared the same way as small ones.
A declaration computes 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.constraint("supply", Sum(W, x[P, W]) <= 1.0)
m.constraint("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 contains no matrix. The parameter is the data of the
caller. The model has added two constraints and an objective, each a symbolic
expression recording the variable, the coefficient and the sets that 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. The matrix therefore
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.constraint("supply", Sum(W, x[P, W]) <= 1.0)
m.constraint("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 the rows of a constraint as a slice. A caller locates a
named row block 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()`. Both return
labeled arrays over the sets, not offsets into a vector.

A variable over a subset keeps its values sparse. Reading a solution over a
sparse network therefore builds no grid.

## What a variable costs

A member's column is computed from its multi-index and is not stored. A
variable over millions of columns costs its members, not its columns.
`n_columns` above is twenty thousand, and the variable stores a few numbers:
the start of its block and the sizes of its sets.

## Progress

A large model 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)
```

A solve stopped at `time_limit` reports `feasible` True where the solver
found a point. `objective`, `primal`, `bound` and `gap` then read that
point. See [Solving](/tutorial/solving).

---

# /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. The three periods of each unit take the capacity of that unit.

## A bound has a value at every column

A bound with no value at some member of the variable raises `ValueError`,
and 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' has no value at member ('b',) of variable 'x'; give the bound a value at every member of the variable
```

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 is not over ['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 contains 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.constraint("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, and a linear
model does not require one. 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; pass operands that share a dimension
```

**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`, and the message gives the
coordinate. A quotient of infinity is not passed to a solver. The check
covers a Python number, a NumPy scalar and a coefficient with a zero at any
coordinate. 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 is zero at 1 coordinate(s), first at {'G': 'base', 'T': 1}; remove the zeros or divide by another parameter
```

## 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.constraint("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
```

## Forms that raise

A line of modeling arithmetic produces a linear term, or raises with a
message that gives the form to write instead. There is no third outcome.

| Written | Form to write instead |
| --- | --- |
| `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 is passed to 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 constraint |
| `0.0 <= x[P] <= 1.0` | each bound as its own constraint |
| `x[P] + c[P]` | a coefficient has no row until a variable multiplies it |
| `c[P] * (x[P] + 1.0)` | a coefficient times a constant is one value per row; the constant goes to the right-hand side |
| `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 with a `nimopt` message, not with a bare Python error. The test
suite executes the table and the replacement forms the messages give.

```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 apply the operator and
return an array of expressions. The expression types reject the ufunc, and
the message describes how a coefficient is built.

---

# /guides/conditions

# Conditions on a sum and on a constraint

Two cases require a condition. A constraint may sum over part of the members
of a variable, such as the arcs of a network where the variable is indexed
over the full product. A constraint may also apply to some members of its
frame only, such as a capacity limit on one plant. `where=` covers both
cases. `over=` declares the rows of a constraint 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.constraint("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
gives three of them a coefficient. A variable over a subset would have three
columns from the start. Use a condition where the variable is over the
product and one constraint reads part of it. Declare a subset where the model
never uses the other members.

## Restricting the rows

`m.constraint(..., 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.constraint("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 frame of the constraint raises
`ValueError`, and 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.constraint("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',)
```

## Declaring 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. A row missing one of its
terms would express a constraint that was not written.

`over=domain` declares 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.constraint("capacity", Sum(W, x[P, W]) <= 10.0, over=product((P,)))

print(rows.n_rows)
```

Output:

```text
2
```

## `over` or `where`, not both

`over=` declares the rows and `where=` restricts them. 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.constraint(
    "capacity",
    Sum(W, x[P, W]) <= 10.0,
    where=subset((P,), {"P": np.array(["p0"])}),
    over=product((P,)),
)
```

Raises ValueError:

```text
ValueError: constraint 'capacity' is given over= and where= together; pass one of them
```

A condition on a sum and a condition on the constraint compose. The first
restricts what is summed, and the second restricts 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, such as the stored energy at the first period or 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`, and the reference is indexed over `G` alone. An
initial condition is one row per unit, referencing the column of that unit 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.constraint("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.constraint("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'; read it at a member that set contains
```

---

# /guides/highs-methods

# Interior point and first-order methods

`method=` selects the algorithm HiGHS runs on the assembled matrix. Four
values select one algorithm each, and three further options configure the
interior point and first-order methods. Every setting on this page runs in
the test suite against the bundled models, and each solve returns the optimum
of the model. Which method is fastest or smallest depends on the model.

| 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, the HiGHS interior point method built on a
preconditioned conjugate gradient. `hipo` runs HiPO, an interior point method
built on a direct factorization of the Newton system, parallel across the
elimination tree. `newton_system` selects the augmented system or the normal
equations, and `choose` leaves the selection to the solver. HiPO stores that
factorization in memory and requires more memory than IPX on the same model.
`threads` applies to its iterations and not to the crossover that may follow
them. `pdlp` runs cuPDLP-C, a primal-dual hybrid gradient method that reads
the matrix only through matrix-vector products. `crossover` controls whether
an interior point method passes its solution to the simplex method to obtain
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 support `crossover` and `method` up to `barrier`.
`newton_system`, `pdlp_tol`, `hipo` and `pdlp` are specific to HiGHS. Passing
one of them to another solver raises and reports the option. The
[solvers reference](/reference/solvers) records the same rule.

## 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. A basis, an exact active set and
duals at a vertex require that move. Crossover 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 the same way. A constraint that is tight at the optimum may be
a tolerance away from equality, and a variable at a bound may be a tolerance
inside it. The interior point is sufficient where the solution is read as
quantities and prices and not as a basis.

## PDLP on a GPU

A first-order method computes no factorization. Its memory is the matrix in
two orientations, one for each product, plus working vectors of the row and
column dimensions. That memory grows linearly with the problem, and a problem
too large for a factorization 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 stored
there while presolve runs, and that is the peak on the host. The reduced
problem is transferred to the GPU. It is smaller than the problem as
declared: presolve removes the rows that bound a single column and the
columns it can fix. The GPU therefore needs memory for the reduced matrix
twice, the working vectors, and the buffers of the sparse kernels, and
nothing else. The host keeps the original and the reduced LP throughout, and
postsolve maps the solution back through them.

A HiGHS built without CUDA runs the same method on the CPU. The solver log,
requested with `log=True`, reports the device the method 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 from feasibility. A
tighter tolerance costs more iterations, and each iteration is a pass over
the matrix. The iterations are cheap and numerous, and the tolerance controls
both the solve time and the accuracy of the solution.

HiGHS checks the point PDLP returns against its own `feasibility_tol` and
`optimality_tol` after postsolve. That check is stricter than the criterion
of PDLP. PDLP measures its residuals relative to the scaled problem it
iterates on. 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 fails the HiGHS check reports the status HiGHS calls unknown. `nimopt`
raises on that status and reports it, and it reads no values from that point.
On a model of any size the default `pdlp_tol` is usually too loose for HiGHS
to accept the point. The two agree at a tolerance one to two orders
tighter. Tighten `pdlp_tol` until the point passes, or loosen the feasibility
and optimality tolerances to the accuracy the model requires.

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. The method returns the last row of the log, not the best row.

```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 with HiPO and GPU support

HiGHS 1.15.1 stores the HiPO orderings and its BLAS in a separate library,
`libhighs_extras`. `libhighs` loads that library at run time by name. The
`highspy` wheel on PyPI and the conda-forge package ship without that library
and without CUDA. Given `hipo`, such a HiGHS logs an error and runs simplex.
`nimopt` raises instead and reports the missing library. Given `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. That directory 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 is sufficient.

   ```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 the HiGHS build
   configuration enables HiPO. Pass the prefix of the conda environment
   through the `CMAKE_PREFIX_PATH` environment variable, not through
   `CMAKE_ARGS`. `CMAKE_ARGS` 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 `libhighs` of the wheel searches its own
   directory for the extras library. Copy the built `libhighs_extras.so`
   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. The conda
   environment that provides them must remain installed. The `libhighs` and
   `libcudalin` of the GPU wheel require the `lib` directory of that
   environment in their own run path. `patchelf --set-rpath` sets it on the
   unpacked wheel before the wheel is packed again. The run path of the
   extension module alone does not resolve the CUDA libraries.

   ```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/
   ```

A solve with `method="hipo"` verifies the installation. A HiGHS without the
extras library raises before anything is solved. A HiGHS with the library
returns the optimum.

---

# /guides/lags

# Lags

Time-coupled constraints reference the previous period. A storage balance
relates the stored energy at `t` to the stored energy 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.constraint("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, and 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.constraint("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, and the level at the end of the horizon enters the row of the first
period.

## 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'] and takes the set, not a lag of it; write 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']; write the lag at the variable's reference
```

A lag is an integer number of members. A fractional lag raises `ValueError`
and is not 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/piecewise

# Piecewise-linear curves

A fuel cost that rises in steps, an efficiency that changes with load, and a
revenue that saturates are each a curve through a list of points.
`piecewise` relates an expression `y` to an expression `x` through such
points. `x` is on the curve. `sign` compares `y` with the curve: `==` sets
`y` to the curve, `>=` bounds `y` below by it and `<=` bounds `y` above by it.

The points are two parameters, `x_points` and `y_points`. Each is over some
or all of the sets of `x` and over one breakpoint set. The declaration
generates the variables and the constraints of its method. Each generated
name begins with the declaration's name.

## A cost curve that is not convex

`method="incremental"` is exact for breakpoints that are strictly increasing
or strictly decreasing. It adds one continuous variable and one integer
variable per segment.

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

G = Set("G", np.array(["a", "b"]))
B = Set("B", np.array(["b0", "b1", "b2", "b3"]))
power = Param.from_dense("power", (G, B), [[0, 10, 20, 30], [0, 10, 20, 30]])
cost = Param.from_dense("cost", (G, B), [[0, 5, 30, 35], [0, 20, 30, 40]])

m = Model("curve")
p = m.var("p", (G,))
c = m.var("c", (G,))
m.constraint("demand", Sum(G, p[G]) == 25.0)
m.piecewise(
    "fuel",
    x=p[G],
    x_points=power[G, B],
    y=c[G],
    y_points=cost[G, B],
    sign=">=",
    method="incremental",
)
m.set_objective(Sum(G, c[G]))

s = m.solve()
print(s.objective)
print(s.primal("p").to_dense())
print(list(m.variables))
print(list(m.constraints))
```

Output:

```text
30.0
[10. 15.]
['p', 'c', 'fuel_fill', 'fuel_order']
['demand', 'fuel_x', 'fuel_y', 'fuel_order_bound', 'fuel_fill_order', 'fuel_order_link']
```

The generated variables and constraints are declarations of the model.
`primal`, `dual`, `row` and `explain` read them by their names.

`x` and `y` are expressions. `p[G] + 5.0` is on the curve five units above
`p`, and the incremental method reads that constant. The tangent method
multiplies `x` by the slope of each segment, and a constant in `x` raises
`ValueError`. Subtract it from `x_points` instead.

## A convex curve

`method="tangent"` adds no variable. It adds one row per segment, and two
rows that keep `x` between the first and the last breakpoint. The rows
describe the curve exactly when the points are convex under `>=`, or
concave under `<=`.

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

G = Set("G", np.array(["a"]))
B = Set("B", np.array(["b0", "b1", "b2"]))
power = Param.from_dense("power", (G, B), [[0, 10, 20]])
cost = Param.from_dense("cost", (G, B), [[0, 10, 30]])

m = Model("convex")
p = m.var("p", (G,))
c = m.var("c", (G,))
m.constraint("demand", Sum(G, p[G]) == 15.0)
m.piecewise("fuel", p[G], power[G, B], c[G], cost[G, B], ">=", "tangent")
m.set_objective(Sum(G, c[G]))

s = m.solve()
print(s.objective)
print(m.constraints["fuel_tangent"])
```

Output:

```text
20.0
Constraint('fuel_tangent', ('G', 'fuel_segment'), 2 rows, 4 coefficients)
```

Points that are not convex under `>=` raise `ValueError` and identify the
first entity at fault. The tangent rows of such points describe a different
curve.

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

G = Set("G", np.array(["a"]))
B = Set("B", np.array(["b0", "b1", "b2"]))
power = Param.from_dense("power", (G, B), [[0, 10, 20]])
cost = Param.from_dense("cost", (G, B), [[0, 20, 30]])

m = Model("concave")
p = m.var("p", (G,))
c = m.var("c", (G,))
m.piecewise("fuel", p[G], power[G, B], c[G], cost[G, B], ">=", "tangent")
```

Raises ValueError:

```text
ValueError: piecewise 'fuel' has points that are not convex, required by sign '>=' at {'G': 'a'}; use method 'incremental'
```

## One curve for every entity

`x_points` and `y_points` are over the breakpoint set and over as many of the
sets of `x` as the curves differ along. Points over the breakpoint set alone
give every entity the same curve.

## Entities with fewer breakpoints

A parameter built from a table has entries only where the table lists
them. An entity with fewer breakpoints lists its first ones, and the rest
are absent. An entity with no breakpoint has no generated rows and no
generated columns. Its `x` and `y` are not related by the declaration.

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

G = Set("G", np.array(["a", "b", "c"]))
B = Set("B", np.array(["b0", "b1", "b2"]))
labels = {
    "G": np.array(["a", "a", "a", "b", "b"]),
    "B": np.array(["b0", "b1", "b2", "b0", "b1"]),
}
power = Param.from_long("power", (G, B), labels, [0.0, 10.0, 20.0, 0.0, 15.0])
cost = Param.from_long("cost", (G, B), labels, [0.0, 10.0, 30.0, 0.0, 30.0])

m = Model("short")
p = m.var("p", (G,))
c = m.var("c", (G,))
m.piecewise("fuel", p[G], power[G, B], c[G], cost[G, B], ">=", "incremental")

print(m.variables["fuel_fill"].n_columns)
print(m.constraints["fuel_x"].n_rows)
```

Output:

```text
3
2
```

A breakpoint absent before a present one raises `ValueError`. A table
with an entity of one breakpoint raises `ValueError`.

## A curve that a binary variable switches off

`active=` takes a binary variable over the sets of `x`, or a sum of them.
Where it is 1, `x` is on the curve. Where it is 0, `x` is 0 and `y` is
compared with 0. The curve then starts at its first breakpoint, and a first
breakpoint above 0 is a minimum output. `active=` is supported by
`method="incremental"` only. A variable with bounds outside 0 and 1 and a
scaled variable raise `ValueError`.

A continuous variable raises `ValueError` under the default. A value between
0 and 1 scales every breakpoint, so the curve is met at a fraction of its
first breakpoint and at a fraction of its cost. `relaxed=True` accepts that
variable and declares the scaled curve, which is the linear relaxation of
the switch. `relaxed=True` with no `active=` raises `ValueError`.

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

G = Set("G", np.array(["a"]))
T = Set("T", np.array(["t0", "t1"]))
B = Set("B", np.array(["b0", "b1", "b2"]))
power = Param.from_dense("power", (G, B), [[10, 20, 30]])
cost = Param.from_dense("cost", (G, B), [[10, 15, 30]])
demand = Param.from_dense("demand", (T,), [5.0, 25.0])

m = Model("committed")
p = m.var("p", (G, T))
c = m.var("c", (G, T))
on = m.var("on", (G, T), upper=1.0, integer=True)
spot = m.var("spot", (T,))
m.constraint("balance", Sum(G, p[G, T]) + spot[T] == demand[T])
m.piecewise(
    "fuel",
    x=p[G, T],
    x_points=power[G, B],
    y=c[G, T],
    y_points=cost[G, B],
    sign=">=",
    method="incremental",
    active=on[G, T],
)
m.set_objective(Sum(G, T, c[G, T]) + 1.2 * Sum(T, spot[T]) + 2.0 * Sum(G, T, on[G, T]))

s = m.solve()
print(s.objective)
print(s.primal("on").to_dense())
print(s.primal("p").to_dense())
```

Output:

```text
29.0
[[0. 1.]]
[[ 0. 20.]]
```

In `t0` the demand of 5 is below the minimum output of 10, and the unit is
off.

## A declaration over some of the entities

`where` restricts a declaration to some entities of its breakpoints. It is
a parameter, a tuple of sets or a domain over the sets of `x_points` other
than the breakpoint set, the form of `where` on `Model.constraint`. The
breakpoint checks, the generated columns and the generated rows cover the
entities at its coordinates. `x`, `y` and `active` are compared at those
coordinates only.

Two declarations can share the breakpoints and split the entities. In the
example, `a` has a status variable and points that are not convex, `b` has
convex points, and `c` has no curve. `on` declares the curve of `a` with the
incremental method and `active`. `free` declares the curve of `b` with the
tangent method. A linear term prices the output of `c`. The variable `fuel`
is declared over the entities with a curve, and `p` over every entity.

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

G = Set("G", np.array(["a", "b", "c"]))
T = Set("T", np.array(["t0", "t1"]))
B = Set("B", np.array(["b0", "b1", "b2"]))
points = {
    "G": np.array(["a", "a", "a", "b", "b", "b"]),
    "B": np.array(["b0", "b1", "b2", "b0", "b1", "b2"]),
}
power = Param.from_long("power", (G, B), points, [10.0, 20.0, 30.0, 0.0, 10.0, 20.0])
cost = Param.from_long("cost", (G, B), points, [10.0, 25.0, 30.0, 0.0, 10.0, 30.0])
committed = Param.from_long("committed", (G,), {"G": np.array(["a"])}, [1.0])
flexible = Param.from_long("flexible", (G,), {"G": np.array(["b"])}, [1.0])
price = Param.from_long("price", (G,), {"G": np.array(["c"])}, [3.0])
demand = Param.from_dense("demand", (T,), [5.0, 25.0])
curves = subset(
    (G, T),
    {"G": np.array(["a", "a", "b", "b"]), "T": np.array(["t0", "t1", "t0", "t1"])},
)
status = subset((G, T), {"G": np.array(["a", "a"]), "T": np.array(["t0", "t1"])})

m = Model("split")
p = m.var("p", (G, T))
fuel = m.var("fuel", (G, T), subset=curves)
u = m.var("u", (G, T), subset=status, upper=1.0, integer=True)
m.constraint("balance", Sum(G, p[G, T]) == demand[T])
m.piecewise(
    "on", p[G, T], power[G, B], fuel[G, T], cost[G, B], ">=", "incremental",
    active=u[G, T], where=committed,
)
m.piecewise(
    "free", p[G, T], power[G, B], fuel[G, T], cost[G, B], ">=", "tangent",
    where=flexible,
)
m.set_objective(
    Sum(G, T, fuel[G, T]) + Sum(G, T, price[G] * p[G, T]) + 2.0 * Sum(G, T, u[G, T])
)

s = m.solve()
print(s.objective)
print(m.variables["on_fill"].n_columns, m.constraints["free_tangent"].n_rows)
```

Output:

```text
34.5
4 4
```

`on_fill` has one column per segment of `a` and timestep. `free_tangent`
has one row per segment of `b` and timestep.

## Saving a piecewise declaration

A model file of version 4 contains the declaration under `piecewise`. The
generated variables, constraints and parameters are not written. Loading
the file generates them again.

```python
from nimopt import Definition, Sum

d = Definition("curve", sense="min")
G, B = d.set("G"), d.set("B")
power, cost = d.param("power", (G, B)), d.param("cost", (G, B))
p = d.var("p", (G,))
c = d.var("c", (G,))
d.constraint("demand", Sum(G, p[G]) == 25.0)
d.piecewise("fuel", p[G], power[G, B], c[G], cost[G, B], ">=", "incremental")
d.set_objective(Sum(G, c[G]))

print(d.to_yaml())
```

Output:

```text
version: 4
name: curve
sense: min
sets: [G, B]
parameters:
  power: [G, B]
  cost: [G, B]
variables:
  p:
    sets: [G]
  c:
    sets: [G]
constraints:
  demand:
    relation: Sum(G, p[G]) == 25
piecewise:
  fuel:
    x: p[G]
    x_points: power[G, B]
    y: c[G]
    y_points: cost[G, B]
    sign: '>='
    method: incremental
objective: Sum(G, c[G])
```

`version=3` writes a model in the format a reader of version 3 accepts. The
generated declarations are written in place of the declaration, and the
breakpoints are not written. A model loaded from that file contains the
same rows and no piecewise declaration. A definition has no data to
generate the rows from, and writing it as version 3 raises `ValueError`.

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

G = Set("G", np.array(["a"]))
B = Set("B", np.array(["b0", "b1"]))
power = Param.from_dense("power", (G, B), [[0, 10]])
cost = Param.from_dense("cost", (G, B), [[0, 20]])

m = Model("line")
p = m.var("p", (G,))
c = m.var("c", (G,))
m.piecewise("fuel", p[G], power[G, B], c[G], cost[G, B], ">=", "tangent")

print(m.to_yaml(version=3))
```

Output:

```text
version: 3
name: line
sense: min
sets: [G, fuel_segment]
parameters:
  fuel_slope: [G, fuel_segment]
  fuel_intercept: [G, fuel_segment]
  fuel_x_low: [G]
  fuel_x_high: [G]
variables:
  p:
    sets: [G]
  c:
    sets: [G]
constraints:
  fuel_tangent:
    relation: c[G] - fuel_slope[G, fuel_segment] * p[G] >= fuel_intercept[G, fuel_segment]
  fuel_x_min:
    relation: p[G] >= fuel_x_low[G]
  fuel_x_max:
    relation: p[G] <= fuel_x_high[G]
```

---

# /guides/saving-and-loading

# Saving and loading a model

A definition writes a YAML file whose expressions are written as they are
typed in Python. The file is the canonical form of the model: every derived
coefficient written out, every term with its own sum and sign, and 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.constraint("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: 4
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 schema of the data: every set, and every
parameter with its dimensions. `build` raises `ValueError` for a mapping that
omits one of them.

## 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: 3
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. An
edit is accepted where Python accepts it, and is normalized the same way.
Adding a scalar, reordering terms and reversing a comparison all parse. The
file written back is the canonical form.

```python
from nimopt import loads

edited = """
version: 3
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 that raises in Python raises here with the same message. A construct
outside the expression syntax raises and reports it.

```python raises=ValueError
from nimopt import loads

loads(
    """
version: 3
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': the syntax supports one call; write Sum
```

## 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, and a dense parameter is nested lists. A parameter with values at some
coordinates of its product is a table of the dimensions then `value`. A file
containing 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
```

## Datetime members

A set whose members are `datetime64` or `timedelta64` round trips through both
data sources, in every unit. Inline, such a set is written as a mapping of
`dtype` and `members`: a `datetime64` member as its ISO 8601 string, and a
`timedelta64` member as its integer count of the unit in the `dtype`. A member
fixed in a relation is written as quoted text.

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

T = Set("T", np.array(["2030-01-01T00", "2030-01-01T01"], dtype="datetime64[h]"))
cost = Param.from_dense("cost", (T,), np.array([2.0, 5.0]))
m = Model("dispatch", sense="min")
gen = m.var("gen", (T,), upper=10.0)
m.constraint("start", gen["2030-01-01T00"] == 4.0)
m.constraint("total", Sum(T, gen[T]) >= 6.0)
m.set_objective(Sum(T, cost[T] * gen[T]))
text = m.to_yaml(inline=True)
print(text[text.index("constraints:") :])
print(loads(text).solve().objective)
```

Output:

```text
constraints:
  start:
    relation: gen['2030-01-01T00'] == 4
  total:
    relation: Sum(T, gen[T]) >= 6
objective: Sum(T, cost[T] * gen[T])
data:
  T:
    dtype: datetime64[h]
    members: [2030-01-01T00, 2030-01-01T01]
  cost: [2.0, 5.0]

18.0
```

The fixed member is written as a string and is read back to the same member.
A string is parsed as ISO 8601, and a `datetime.datetime`, a `datetime.date`
and a `datetime64` of another unit are converted. An integer against a
`timedelta64` set is a count of that set's own unit. A conversion that is not
exact raises `ValueError` instead of truncating.

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

T = Set("T", np.array(["2030-01-01T00", "2030-01-01T01"], dtype="datetime64[h]"))
m = Model("dispatch")
gen = m.var("gen", (T,))
gen["2030-01-01T00:30"]
```

Raises ValueError:

```text
ValueError: member '2030-01-01T00:30' does not convert exactly to datetime64[h] at dimension 'T' of variable 'gen'; write a member in the unit of that dimension
```

`Model.row` takes the same text for a datetime dimension, and a `Row` displays
each datetime coordinate in it.

## Data beside the file, for a large model

`save` writes a model's file and its data as an `.npz` beside it, under the
file's stem. The file records the sidecar name. `load` reads both. The file
records no path and no machine name.

```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 contains
data and a `data=` together raises `ValueError`: one model takes one data
source.

## A file that describes its own format

`instructions=True` on `save`, `Definition.to_yaml` and `Model.to_yaml` writes
a comment block at the top of the file. The block is the same in every file.
It describes the format, not the model: the keys and their order, the
defaults, the rules that determine which rows a constraint has, and the
expression syntax. A reader with one file interprets it without this package.

The block is a YAML comment. A file with it loads to the same model as one
without it, and writing the loaded model with the flag gives the same text.

```python
from nimopt import Definition, loads

d = Definition("dispatch", sense="min")
T = d.set("T")
load, gen = d.param("load", (T,)), d.var("gen", (T,))
d.constraint("balance", gen[T] == load[T])
text = d.to_yaml(instructions=True)
print("\n".join(text.splitlines()[:5]))
print(loads(text).to_yaml() == d.to_yaml())
print(loads(text).to_yaml(instructions=True) == text)
```

Output:

```text
# --- Reading this file --------------------------------------------------
# A nimopt model file, format version 4. The keys are written in this
# order, and no other key is accepted: version, name, sense, sets,
# aliases, parameters, variables, constraints, piecewise, objective,
# data. Only version, name and sense are required.
True
True
```

---

# /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, and `x` has three columns.
The assembled matrix is empty: no constraint has been added. A model over a
sparse network costs its arcs, not 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`
resolves no label when positions are already known.

## Constraints over a subset variable

A sum over a subset variable runs over the members the variable has. A row
therefore 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.constraint("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. A bound per
plant therefore applies to each arc of that plant.

```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 incurs a no-load cost. An uncommitted
unit produces nothing. The `capacity` and `minimum` rows are written against
the binary column. This is the only MILP in the corpus.

```text
minimize    Σ_{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. `reference` enumerates every on-off subset per
snapshot and takes the cheapest feasible one. The optimum is 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
minimize    Σ_{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 requires none,
and the corpus writes no coefficient a model does not require.

```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, and the optimum is the merit order per snapshot.
`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, and `absent`
reports no dropped row.

```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 given by the scenario. The shape of the first-stage variable expresses
the information available to it: `cap` is indexed by technology alone, and one
capacity applies to every scenario. `p` and `shed` are indexed by scenario as
well and may differ across it.

```text
minimize    Σ_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. The second and third sums are
therefore an expectation, and the objective is the committed capital plus the
expected cost of the recourse. No row relates the capacity of one scenario to
the capacity of another: the missing dimension expresses
non-anticipativity.

```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 declared over all three dimensions,
and the variable it bounds has one. The rows of every scenario read the same
column, and the capacity is therefore 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 no row couples
one hour to the next. The recourse in each scenario-hour is therefore the
merit order of the built capacity against that demand, with the remainder
unserved at `voll`. The capacity is written as bands: the band of the
cheapest technology, then the next, and last the band between the most
expensive technology and lost load. The total separates into one term per
band. Each term is convex in the level of its band and changes slope only at
a demand. `reference` minimizes the terms one at a time and adds them, and
the result is arithmetic over the inputs, not 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 describes a cost frontier: capital falls as marginal cost rises,
and lost load is more expensive than the most expensive technology. All three
technologies are built. The capacity stacks to 180 and leaves the peak hour
of 225 in the cold scenario short. Shedding 45 costs less than a fourth band
used 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 applies only to data whose technologies form a frontier and
whose bands stack. `reference` raises on any other data, and it returns no
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; pass a capital cost that falls across the merit order
```

---

# /models/fleet

# Fleet

`nimopt.models.fleet` is the same problem as `dispatch`. It is declared as
one variable per unit over the snapshots alone, and the terms of the units
are added into one balance row. The optimum is the same merit order, and the
cost of the declaration differs.

```text
minimize    Σ_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: the number of variables is a property of
the declaration, not 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 functions have different callers. 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
and calls nothing from `nimopt`. A formulation error is therefore 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 lag |
| [`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` is located at bus
`b`. Multiplying the generation by it maps a row over generators to a row
over buses. The coefficient introduces `B`, a dimension no variable has, and
the balance is indexed over the dimensions the lookup defines.

```text
minimize    Σ_{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. 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. A solar unit is bounded by its
hourly availability, and a thermal unit by its rating.

```text
minimize    Σ_{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, and 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. The optimum is the merit order against the
capacities of that hour.

```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 given. `on` is indexed by hour and unit alone, and `p` and `shed` are
indexed by the scenario as well. One commitment applies to every scenario,
and the binary column is therefore taken under uncertainty.

```text
minimize    Σ_{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 the same in every scenario. The second and third sums
are weighted by the probability of a scenario, and 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])
```

No row couples one hour to the next, and 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. That enumeration 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 the model has no sink for
unwanted energy. A unit whose minimum exceeds the demand of the mildest
scenario is therefore not committed: its output would violate the balance row
of that scenario. The first hour demands 38.25 in the mildest scenario. The
minimum of `base` is 40.0 and its fuel is the cheapest of the three, and it
is not committed in that hour. The third hour is the opposite case: the fleet
supplies 260.0 against a coldest demand of 268.75, every unit 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. `base` serves that demand at
the lowest cost. A model given that one number commits `base`, and that
commitment violates the balance row of the mildest scenario, whose
probability is 0.5. The scenario dimension on `demand` excludes the
commitment, and the missing scenario dimension on `on` applies that exclusion
to every scenario.

---

# /models/sector

# Sector

`nimopt.models.sector` has mixed density. The region-technology map is
sparse: a technology exists in some regions and not in others. 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. It is
therefore sparse in one axis and dense in the other.

```text
minimize    Σ_{(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. 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. The capacity bound applies to the sited
pairs alone, and 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`. The row at the first hour therefore references
the last hour, and every hour has a row. The ramp row references `T - 1`. The
first hour has no predecessor, and that row is not produced. The model
exercises both lag rules.

```text
minimize    Σ_{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 and not a bound, and each therefore 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 costs of the fleet span 50.0 to 55.0. Shifting energy through the store
costs more than it saves. The store is therefore idle, and the optimum is the
hourly merit order. A store that cycles requires data with a wider cost
spread, and 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: that lag wraps. The ramp row at the
first hour is not produced: that lag does not wrap.

```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 and not all of them. The
cost parameter has one entry per arc, and the flow variable takes its members
from that parameter. The model therefore has one column per arc, not one per
cell of the plant-warehouse product.

```text
minimize    Σ_{(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 the band of a plant. No supply row binds
therefore, and each warehouse buys from the cheapest plant connected to 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. No plant serves the last
warehouse, and it 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)
```
