refactor: three modules in the engine, one per box in the diagram

ARCHITECTURE.md

@@ -31,10 +31,11 @@ flowchart TB
     subgraph REL["Relational lane — streaming · memory-bounded · linopy-free"]
         direction TB
         LOWER["lowering.py"] --> PLAN["logical plan<br/>(relational/plan.py)"]
-        DR[("data<br/>parquet paths / Arrow tables")] --> EXEC
-        PLAN --> EXEC["executor.py<br/>tidy tables in file-backed duckdb<br/>under memory_limit"]
-        EXEC --> LPS["lp_file sink<br/>portability, debugging<br/>(mps planned)"]
-        EXEC --> DIRECT["solver_direct sink<br/>COO batches → highspy → HiGHS"]
+        PLAN --> COMP["compiler.py<br/>plan → SQL text<br/>pure: no connection"]
+        DR[("data<br/>parquet paths / pandas")] --> EXEC
+        COMP --> EXEC["executor.py<br/>tidy tables in file-backed duckdb<br/>under memory_limit"]
+        EXEC --> LPS["lp_file sink (sinks.py)<br/>portability, debugging<br/>(mps planned)"]
+        EXEC --> DIRECT["solver_direct sink (sinks.py)<br/>COO batches → highspy → HiGHS"]
         DIRECT --> SOL["solution tables<br/>(label join, never dense)"]
     end
 
@@ -195,6 +196,14 @@ sink, is [ROADMAP Track 4](ROADMAP.md#track-4--sink-capabilities).
 
 ## The relational lane
 
+**Three modules, one per box above.** `compiler.py` turns plan nodes into SQL
+and holds no connection; `executor.py` owns the database and fills the tables;
+`sinks.py` drains them. The split is what makes the admissibility test below
+something you can perform rather than reason about — build a `SqlCompiler`,
+hand it a node, read the `SELECT` (`tests/test_compiler.py` does exactly that,
+with no engine installed). It is also why a new sink is a function in one file
+instead of another method on the executor.
+
 **Tidy tables.** Parameters are `(dims…, value)`; a variable frame is
 `(dims…, var_label)`, one row per *existing* variable; a linear expression is
 `(frame dims…, var_label, coeff)` plus a constant part; constraint rows are
@@ -274,7 +283,9 @@ than the plan level, which is
 | `errors.py` | the exception hierarchy; the one module the engine may import |
 | `relational/plan.py` | frozen logical-plan dataclasses |
 | `relational/arrow.py` | the Arrow boundary — caller tables in, via the PyCapsule protocol |
-| `relational/executor.py` | duckdb execution + `lp_file` / `solver_direct` sinks |
+| `relational/compiler.py` | plan → SQL text; pure, no connection |
+| `relational/executor.py` | duckdb: bind sources, label, assemble the tables |
+| `relational/sinks/` | how a built model leaves: `lp_file`, `solver_direct` (one module each, [README](linopy_yaml/relational/sinks/README.md)) |
 | `compat/__init__.py` | opt-in shim: `build` / `extend` on a `linopy.Model` |
 | `compat/loader.py` | data coercion to `xr.Dataset`, master coords |
 | `compat/builder.py` | eager backend: core AST → `linopy.Model` |

linopy_yaml/relational/sinks/README.md

@@ -0,0 +1,69 @@
+# Sinks
+
+How a built model leaves the engine. These are the two boxes downstream of the
+executor in [ARCHITECTURE.md](../../../ARCHITECTURE.md)'s pipeline.
+
+| Module | Sink | Needs |
+|---|---|---|
+| `tables.py` | — | the contract every sink reads |
+| `lp_file.py` | `lp_file` — LP text | nothing beyond duckdb |
+| `highs.py` | `solver_direct` — COO batches → HiGHS | `highspy` |
+
+## The contract
+
+A sink takes a `ModelTables` and nothing else: a connection holding `cols`
+(col, lb, ub, vtype), `obj` (col, coeff), `rows` (row, sense, rhs) and `A`
+(row, col, coeff), plus the counts it chunks by and the objective's sense and
+constant — those last two live outside the tables because a constant has no
+column to attach to.
+
+A sink never learns how the tables were filled, and the executor never learns
+how they are drained. That is the point: adding `mps` is a new module here, not
+another method on `DuckdbExecutor`.
+
+## Adding a sink
+
+1. New module, named for the sink. Take `ModelTables`, return whatever the sink
+   naturally returns (`None` for a writer, `(status, objective)` for a solver).
+2. Re-export it from `__init__.py`.
+3. Thin delegation on `DuckdbExecutor` — three lines, no logic.
+4. If it needs an optional dependency, import it **inside the function**. The
+   module boundary is the fence; the lazy import is what keeps importing this
+   package free for callers who will never use that sink.
+5. Stream. Nothing here may materialise the model — hard rule 4. Aggregate
+   inside duckdb, or hand the solver batches.
+
+## Why one module per sink
+
+Because that is where the fences are. `highspy` is an optional dependency of
+`solver_direct` alone, and a caller that only writes LP files should not import
+it. Splitting by *kind* — all writers together, all solvers together — would
+put two optional imports in one module and a function that branches on which
+solver you meant.
+
+When `mps` lands it may well belong beside `lp_file` (both are chunked `COPY`
+of `printf`'d rows into part files, concatenated bytewise; they would share
+`_cat`, the chunking and the float formatting). That is a decision to take with
+the code in hand, not now — a `text.py` holding one function today would be a
+guess about a sink that does not exist.
+
+## When Track 4 lands
+
+[Track 4](../../../ROADMAP.md#track-4--sink-capabilities) gives each sink a
+declared capability table so `check(model, sink=...)` can answer "will this
+sink take it". Two notes for whoever writes it:
+
+- The capability table belongs next to the sink it describes, in that sink's
+  module. `__init__.py` collects them; it does not own them.
+- Make the set of sinks **closed**, like `helpers.BUILTINS` — not an open
+  `register_sink()`. An installed plugin that can change the answer to
+  `ly.check(model, sink=...)` is hard rule 5's failure mode one level down.
+
+## Known issue
+
+Neither file sink emits in a stable order: several `COPY` statements have no
+`ORDER BY` and `preserve_insertion_order=false` is set on the connection, so
+two runs of the same model produce byte-different files with identical content
+— [#109](https://github.com/FBumann/linopy-yaml/issues/109). `solver_direct` is
+unaffected; every read it makes is explicitly ordered, because `searchsorted`
+requires it.

linopy_yaml/relational/sinks/__init__.py

@@ -0,0 +1,22 @@
+"""Sinks: how a built model leaves the engine.
+
+ARCHITECTURE.md's pipeline draws two boxes downstream of the executor. This
+package is those boxes — one module each, because that is where the fences
+are: ``highspy`` is an optional dependency of ``solver_direct`` alone, and a
+caller that only writes LP files should never import it.
+
+A sink reads :class:`ModelTables` and nothing else. Neither sink knows how the
+tables were filled, and the executor does not know how they are drained, which
+is what makes the planned ``mps`` sink a new module here rather than another
+method on the executor.
+"""
+
+from linopy_yaml.relational.sinks.highs import solve_direct
+from linopy_yaml.relational.sinks.lp_file import write_lp_file
+from linopy_yaml.relational.sinks.tables import ModelTables
+
+__all__ = [
+    'ModelTables',
+    'solve_direct',
+    'write_lp_file',
+]

linopy_yaml/relational/sinks/highs.py

@@ -0,0 +1,94 @@
+"""The ``solver_direct`` sink: COO batches straight into HiGHS.
+
+No float→text→parse round trip — that is the whole reason this exists beside
+:mod:`~linopy_yaml.relational.sinks.lp_file`. Columns arrive as arrow batches,
+rows as numpy slices of ``A``, and the full model never lands in one array.
+
+``highspy`` is imported inside the function rather than at module scope: it is
+an optional dependency, and importing this module must stay free for callers
+that only ever write LP files. The module boundary is the fence; the lazy
+import is what keeps the fence cheap.
+"""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from linopy_yaml.relational.sinks.tables import ModelTables
+
+
+def solve_direct(model: ModelTables, batch_rows: int = 100_000) -> tuple[str, float]:
+    """Stream the model into HiGHS and solve it. Returns ``(status, objective)``.
+
+    Leaves the primal values in a ``sol`` table on the connection, so reading
+    results back stays a label join like every other read — the caller owns
+    the mapping from solver column index to coordinates.
+    """
+    import highspy
+    import numpy as np
+
+    con = model.connection
+    inf = highspy.kHighsInf
+    h = highspy.Highs()
+    h.setOptionValue('output_flag', False)
+
+    empty_i = np.empty(0, dtype=np.int32)
+    empty_f = np.empty(0, dtype=np.float64)
+    reader = con.execute(
+        'SELECT c.col, c.lb, c.ub, c.vtype, COALESCE(o.coeff, 0) AS cost '
+        'FROM cols c LEFT JOIN obj o USING (col) ORDER BY c.col'
+    ).to_arrow_reader(batch_rows)
+    for batch in reader:
+        d = batch.to_pydict()
+        lb = np.nan_to_num(np.asarray(d['lb'], dtype=np.float64), neginf=-inf, posinf=inf)
+        ub = np.nan_to_num(np.asarray(d['ub'], dtype=np.float64), neginf=-inf, posinf=inf)
+        cost = np.asarray(d['cost'], dtype=np.float64)
+        h.addCols(len(cost), cost, lb, ub, 0, empty_i, empty_i, empty_f)
+        variable_type = np.asarray(d['vtype'])
+        noncontinuous = np.flatnonzero(variable_type != 'continuous')
+        if len(noncontinuous):
+            cols_idx = np.asarray(d['col'], dtype=np.int32)[noncontinuous]
+            integrality = np.full(len(noncontinuous), int(highspy.HighsVarType.kInteger), dtype=np.uint8)
+            h.changeColsIntegrality(len(noncontinuous), cols_idx, integrality)
+
+    for lo, hi in model.row_chunks(batch_rows):
+        rows = con.execute(
+            f'SELECT row, sense, rhs FROM rows WHERE row >= {lo} AND row < {hi} ORDER BY row'
+        ).fetchnumpy()
+        a = con.execute(f'SELECT row, col, coeff FROM A WHERE row >= {lo} AND row < {hi} ORDER BY row').fetchnumpy()
+        rhs = np.asarray(rows['rhs'], dtype=np.float64)
+        sense = rows['sense']
+        rlb = np.where(sense == '<=', -inf, rhs)
+        rub = np.where(sense == '>=', inf, rhs)
+        starts = np.searchsorted(np.asarray(a['row']), np.asarray(rows['row'])).astype(np.int32)
+        h.addRows(
+            len(rhs),
+            rlb,
+            rub,
+            len(a['col']),
+            starts,
+            np.asarray(a['col'], dtype=np.int32),
+            np.asarray(a['coeff'], dtype=np.float64),
+        )
+
+    if model.objective_sense == 'max':
+        h.changeObjectiveSense(highspy.ObjSense.kMaximize)
+    h.run()
+
+    status = str(h.getModelStatus()).rsplit('.', 1)[-1].removeprefix('k')
+    objective = h.getInfo().objective_function_value + model.objective_constant
+
+    import pyarrow as pa
+
+    primal = pa.table(
+        {
+            'col': pa.array(np.arange(model.column_count, dtype=np.int64)),
+            'value': pa.array(np.asarray(h.getSolution().col_value, dtype=np.float64)),
+        }
+    )
+    con.execute('DROP TABLE IF EXISTS sol')
+    con.register('sol_src', primal)
+    con.execute('CREATE TABLE sol AS SELECT * FROM sol_src')
+    con.unregister('sol_src')
+    return status, objective

linopy_yaml/relational/sinks/lp_file.py

@@ -0,0 +1,102 @@
+"""The ``lp_file`` sink: the model as LP text.
+
+Portability, debugging, and the differential oracle. Every section is produced
+by a duckdb ``COPY`` into a part file and the parts are concatenated bytewise,
+so the LP text never exists in this process's memory either — only the file
+handle does.
+
+The one hand-managed chunk in the engine that is not label assignment lives
+here: string aggregates do not spill, so the constraint text is emitted in
+fixed row ranges. That costs nothing in a debugging sink.
+
+Block order is not stable run to run (#109); the content is.
+"""
+
+from __future__ import annotations
+
+import shutil
+from pathlib import Path
+from typing import TYPE_CHECKING, Any
+
+if TYPE_CHECKING:
+    from linopy_yaml.relational.sinks.tables import ModelTables
+
+#: Raw lines, no CSV quoting to undo.
+_COPY_OPTS = "(FORMAT csv, HEADER false, QUOTE '', ESCAPE '')"
+
+
+def write_lp_file(model: ModelTables, path: str | Path) -> None:
+    """Write the model as LP text.
+
+    Every section is produced by a duckdb ``COPY`` into a part file, then the
+    parts are concatenated bytewise — so the LP text never exists in this
+    process's memory either, only the file handle does.
+    """
+    path = Path(path)
+    parts = model.workdir / 'lp_parts'
+    parts.mkdir(exist_ok=True)
+    con = model.connection
+
+    con.execute(f"COPY (SELECT printf('%+.17g x%d', coeff, col) FROM obj) TO '{parts / 'obj'}' {_COPY_OPTS}")
+
+    nnz = model.scalar('SELECT count(*) FROM A')
+    avg = max(1, nnz // max(1, model.row_count))
+    con_parts = []
+    for i, (lo, hi) in enumerate(model.row_chunks(max(1, model.chunk_rows // avg))):
+        part = parts / f'cons.{i}'
+        con_parts.append(part)
+        con.execute(
+            f"""
+            COPY (
+                SELECT printf('c%d:', r.row) || chr(10)
+                       || COALESCE(string_agg(printf('%+.17g x%d', a.coeff, a.col), chr(10)), '+0 x0')
+                       || chr(10)
+                       || printf('%s %.17g', CASE r.sense WHEN '==' THEN '=' ELSE r.sense END, r.rhs)
+                FROM rows r LEFT JOIN A a USING (row)
+                WHERE r.row >= {lo} AND r.row < {hi}
+                GROUP BY r.row, r.sense, r.rhs
+            ) TO '{part}' {_COPY_OPTS}
+            """
+        )
+
+    con.execute(
+        f"""
+        COPY (
+            SELECT CASE WHEN lb = '-infinity'::DOUBLE THEN '-infinity' ELSE printf('%.17g', lb) END
+                   || printf(' <= x%d <= ', col)
+                   || CASE WHEN ub = 'infinity'::DOUBLE THEN '+infinity' ELSE printf('%.17g', ub) END
+            FROM cols
+        ) TO '{parts / 'bounds'}' {_COPY_OPTS}
+        """
+    )
+
+    integrality_sections = []
+    for variable_type, keyword in (('binary', 'binary'), ('integer', 'general')):
+        if model.scalar(f"SELECT count(*) FROM cols WHERE vtype = '{variable_type}'"):
+            part = parts / keyword
+            integrality_sections.append((keyword, part))
+            con.execute(
+                f"COPY (SELECT printf('x%d', col) FROM cols WHERE vtype = '{variable_type}') TO '{part}' {_COPY_OPTS}"
+            )
+
+    sense = b'min' if model.objective_sense == 'min' else b'max'
+    with open(path, 'wb') as f:
+        f.write(sense + b'\n\nobj:\n')
+        if model.objective_constant:
+            f.write(f'{model.objective_constant:+.17g}\n'.encode())
+        _cat(f, parts / 'obj')
+        f.write(b'\ns.t.\n\n')
+        for part in con_parts:
+            _cat(f, part)
+        f.write(b'\nbounds\n')
+        _cat(f, parts / 'bounds')
+        for keyword, part in integrality_sections:
+            f.write(f'\n{keyword}\n'.encode())
+            _cat(f, part)
+        f.write(b'\nend\n')
+    shutil.rmtree(parts)
+
+
+def _cat(f: Any, part: Path) -> None:
+    with open(part, 'rb') as src:
+        shutil.copyfileobj(src, f)

linopy_yaml/relational/sinks/tables.py

@@ -0,0 +1,48 @@
+"""What every sink reads, and nothing more.
+
+The contract between the executor and the sinks: four tables in a connection,
+plus the handful of scalars a writer needs to size its own chunking. A sink
+that needs a fifth thing states it here, where both sides can see it.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Any
+
+if TYPE_CHECKING:
+    from collections.abc import Iterator
+    from pathlib import Path
+
+
+@dataclass(frozen=True)
+class ModelTables:
+    """The built model, as a sink sees it.
+
+    ``connection`` holds four tables — ``cols`` (col, lb, ub, vtype), ``obj``
+    (col, coeff), ``rows`` (row, sense, rhs) and ``A`` (row, col, coeff). The
+    scalars alongside are the ones a sink cannot cheaply recover: the counts
+    it chunks by, and the objective's sense and constant, which live outside
+    the tables because a constant has no column to attach to.
+    """
+
+    connection: Any
+    workdir: Path
+    chunk_rows: int
+    column_count: int
+    row_count: int
+    objective_sense: str
+    objective_constant: float
+
+    def scalar(self, sql: str) -> Any:
+        row = self.connection.execute(sql).fetchone()
+        assert row is not None
+        return row[0]
+
+    def row_chunks(self, per_chunk: int) -> Iterator[tuple[int, int]]:
+        """``(lo, hi)`` half-open row ranges covering the constraint matrix."""
+        for lo in range(0, max(self.row_count, 1), per_chunk):
+            hi = min(lo + per_chunk, self.row_count)
+            if hi <= lo:
+                return
+            yield lo, hi