perf(duckdb): the four frames have one schema, and both engines produce it

bench/duckdb-spike.md

@@ -31,11 +31,12 @@ uv run python -m bench.run --arms lpspec duckdb \
 uv run python -m bench.report bench/results/duckdb-spike.jsonl --arms lpspec duckdb
 ```
 
-**Expect the in-tree engine to be slower than the recorded numbers**, and say so
-rather than quietly re-baselining: the port has no counterpart to
-`labels._factored` (§9), so masked models fall back to the counted path and pay
-an ordered window the shipped engine avoided. The committed results file is the
-*old* engine and is kept as the measurement the argument was built on.
+**The committed results file is the *old* engine**, kept as the measurement the
+argument was built on rather than re-baselined. The in-tree engine has since
+grown the `_factored` label path §9 recorded as missing — masked models no
+longer pay an ordered window over the whole coordinate product — which took the
+build on `dispatch/l` from 1.84 s to 0.69 s and on `nodal/l` from 1.09 s to
+0.53 s, against 0.59 s and 0.29 s on polars.
 
 ---
 
@@ -335,8 +336,11 @@ Where the tax actually falls:
   column lists on both sides, alias qualification and null-safe
   `IS NOT DISTINCT FROM` predicates where polars gets all three from
   column-name semantics. Every fragment composition pays it.
-- **`labels._factored`, 37 lines, no counterpart.** The one place polars is
-  *algorithmically* ahead rather than incidentally, and real work to port.
+- **`labels._factored`, 37 lines.** Was the one place polars was
+  *algorithmically* ahead rather than incidentally. Now ported, and the port
+  moved the *choice* — which of the three routes a mask allows — up into
+  `plan.free_prefix`, where it is a question about the plan that both engines
+  ask once rather than two implementations that could drift.
 - SQL **wins** on `_compare` (15 → 4) and on the semi-joins, where
   `WHERE EXISTS` says plainly what `how='semi'` needs a reader to know is not
   a join.
@@ -369,10 +373,10 @@ read-back are written against `ModelTables` and the label frames, so
 ~2,300-line figure priced a replacement; a *second* engine is a compiler and an
 assembler.
 
-**Still open, and the reason to re-read §7 before relying on it:**
-`labels._factored` has no duckdb counterpart, so masked models fall back to the
-counted path and pay an ordered window. The committed results predate the
-in-tree port and do not include that cost.
+**Read §7 as provenance, not as the current cost.** The committed results
+predate the in-tree port, and the port has been measured and improved since —
+the factored label path alone halves the build on every masked model. Re-run
+the command at the top before quoting a ratio.
 
 **Recommended next step, if this is live:** decide the write-path question
 first. If the answer is no, this closes for good and

docs/ARCHITECTURE.md

@@ -348,7 +348,9 @@ in the lane is order-free, which is what lets the query planner rearrange it.
   depending on how much of the product survives the mask — arithmetic, factored,
   counted — and they must agree integer for integer. That agreement is why
   labelling is a module with stated inputs rather than three methods among
-  twenty: nothing else about a build can move an index.
+  twenty: nothing else about a build can move an index. **Which** of the three
+  is a property of the plan, not of an engine, so both engines ask
+  `plan.free_prefix` and only the executions are written twice.
 - The same order comes **back**: `primal` / `dual` / `to_parquet` sort on the
   label before handing rows over, the order the LP sink writes. Stated at the
   read rather than assumed, because a `where` decides which rows survive and a
@@ -413,7 +415,7 @@ than discovered at solve time.
 | `relational/engines/polars/labels.py` | which coordinate gets which solver index; three routes to one number, which must agree |
 | `relational/engines/polars/executor.py` | assemble the model frames from the bound data |
 | `relational/result.py` | what a solve returned: status, objective, and the label joins that read values back |
-| `relational/sinks/tables.py` | what every sink reads and no more — the four frames plus the batching scalars; what an engine produces |
+| `relational/sinks/tables.py` | what every sink reads and no more — the four frames, their dtypes, and the batching scalars; what an engine produces |
 | `relational/sinks/` | how a built model leaves: `lp_file`, `solver_direct` (one module each, [README](https://github.com/FBumann/lpspec/blob/main/src/lpspec/relational/sinks/README.md)) |
 | `linopy/__init__.py` | opt-in shim: `build` / `extend` on a `linopy.Model` |
 | `linopy/loader.py` | data coercion to `xr.Dataset`, master coords |

src/lpspec/relational/engines/duck/compiler.py

@@ -56,8 +56,14 @@ def lit(value: object) -> str:
         return 'TRUE' if value else 'FALSE'
     if value is None:
         return 'NULL'
-    if isinstance(value, float) and not math.isfinite(value):
-        return "'Infinity'::DOUBLE" if value > 0 else "'-Infinity'::DOUBLE"
+    if isinstance(value, float):
+        if not math.isfinite(value):
+            return "'Infinity'::DOUBLE" if value > 0 else "'-Infinity'::DOUBLE"
+        # `::DOUBLE`, or SQL reads `0.1` as DECIMAL(2,1) — a different number
+        # from the double the plan holds, and a different one from what the
+        # polars engine multiplies by. Decimal also propagates: a coefficient
+        # built from one stays decimal all the way into `obj`.
+        return f'{value!r}::DOUBLE'
     return repr(value)
 
 
@@ -346,7 +352,7 @@ def _variable_fragment(self, name: str) -> TermFragment:
         # syntax error rather than an empty projection
         prefix = ''.join(f'{q(d)}, ' for d in dims)
         table = q(self.variables[name])
-        rel = Rel(f'SELECT {prefix}var_label, 1.0 AS coeff FROM {table}', (*dims, 'var_label', 'coeff'))
+        rel = Rel(f'SELECT {prefix}var_label, 1.0::DOUBLE AS coeff FROM {table}', (*dims, 'var_label', 'coeff'))
         masked = self.program.variable(name).where is not None
         presence = Rel(f'SELECT {", ".join(q(d) for d in dims) or q(UNIT)} FROM {table}', dims) if masked else None
         return TermFragment(

src/lpspec/relational/engines/duck/executor.py

@@ -33,10 +33,10 @@
 if TYPE_CHECKING:
     from collections.abc import Iterator, Sequence
 
-_COLS = ('col', 'lb', 'ub', 'vtype')
-_ROWS = ('row', 'sense', 'rhs')
-_MATRIX = ('row', 'col', 'coeff')
-_OBJ = ('col', 'coeff')
+#: The four frames a sink reads, and their dtypes — both stated by
+#: `sinks/tables.py`, which is what reads them. Shared with the polars engine
+#: for the reason the frames are: a sink cannot see which engine filled them.
+_COLS, _OBJ, _ROWS, _MATRIX = sinks.COLS, sinks.OBJ, sinks.ROWS, sinks.MATRIX
 
 
 class _Labels(Mapping[str, 'pl.LazyFrame']):
@@ -91,6 +91,9 @@ def __init__(self) -> None:
         self._program: plan.Program | None = None
         self._var_tables: dict[str, str] = {}
         self._row_tables: dict[str, str] = {}
+        #: the dims each declared name is read through — what `plan.free_prefix`
+        #: needs to know which label path a mask allows
+        self._name_dims: dict[str, tuple[str, ...]] = {}
         self._var_labels = _Labels(self._con, self._var_tables)
         self._row_labels = _Labels(self._con, self._row_tables)
         self._n_cols = 0
@@ -115,10 +118,30 @@ def _view(self, sql: str, prefix: str) -> str:
         return name
 
     def _register(self, name: str, frame: pl.DataFrame) -> str:
-        # the frame itself, not `to_arrow()`: duckdb reads polars natively, and
-        # the round-trip through a pyarrow table is what used to drag pandas
-        # into a runtime that declares it a bridge out and not a dependency
-        self._con.register(name, frame)
+        """Copy *frame* into a table of its own, and return the table's name.
+
+        **Copied, not scanned in place.** A registered frame stays a Python
+        object: duckdb reads it through the buffer protocol from whichever
+        worker thread the scan lands on, and those threads need the GIL — which
+        the thread that called `execute` is holding while it waits. On a plan
+        with several scans of registered frames feeding one pipeline, that
+        deadlocks outright: `transport/l` reproduces it as a build that sits at
+        0% CPU indefinitely, while the same query over copies returns in 0.3 s.
+
+        Copying is not free — the frame and the table are both resident for the
+        length of one statement — but it buys a build that cannot stall, and
+        duckdb's own storage for every scan after the first.
+
+        The frame itself, not `to_arrow()`: duckdb reads polars natively, and
+        the round-trip through a pyarrow table is what used to drag pandas into
+        a runtime that declares it a bridge out and not a dependency.
+        """
+        source = f'__source {name}__'
+        self._con.register(source, frame)
+        try:
+            self._con.execute(f'CREATE TABLE {q(name)} AS SELECT * FROM {q(source)}')
+        finally:
+            self._con.unregister(source)
         return name
 
     def build(self, program: plan.Program, sources: Mapping[str, Any]) -> None:
@@ -129,6 +152,7 @@ def build(self, program: plan.Program, sources: Mapping[str, Any]) -> None:
         is what a caller's data has to survive whoever builds the model.
         """
         self._program = program
+        self._name_dims = plan.name_dims(program)
         bound = bind(program, sources)
         dims = {n: self._register(f'dim_{n}', f.collect()) for n, f in bound.dimensions.items()}
         params = {n: self._register(f'par_{n}', f.collect()) for n, f in bound.parameters.items()}
@@ -170,29 +194,30 @@ def _label_frame(
     ) -> tuple[str, int]:
         """The masked coord product of *dims* with a dense *label* from *start*.
 
-        Two paths, not the polars side's three. **Unmasked**, a row's label is
-        its position in the product — arithmetic on the ordinals, no sort and
-        nothing to count. **Otherwise** it is counted, which costs the ordered
-        window. The `_factored` middle path is an optimisation rather than a
-        semantics, and is left out here on purpose: it is the one place the
-        polars engine is *algorithmically* ahead, so including a half version
-        of it would flatter this port.
+        The same three paths the polars engine has, chosen by the same
+        function. **Unmasked**, a row's label is its position in the product —
+        arithmetic on the ordinals, no sort and nothing to count. **Masked but
+        factoring**, the survivors are a rectangle and the label is arithmetic
+        again over a ranked survivor set (:meth:`_factored`). **Otherwise** it
+        is counted, which costs the ordered window over the whole product.
+
+        `plan.free_prefix` decides between them for both engines, because a
+        label is the solver's own column index: two engines choosing routes
+        independently is how they would come to build different models.
         """
-        rel = self._q.frame(dims, where)
-        if where is None and not restrictions:
-            # row-major: the leading dim's ordinal times the width of the rest
-            terms: list[str] = []
-            stride = 1
-            for d in reversed(dims):
-                terms.append(f'{q(_ordinal(d))} * {stride}')
-                stride *= self._q.cardinality[d]
-            position = ' + '.join(reversed(terms)) if terms else '0'
-            cols = ', '.join(q(d) for d in dims) or q(UNIT)
-            sql = f'SELECT {cols}, ({position} + {start})::BIGINT AS {q(label)} FROM {rel.alias("p")}'
-            name = self._view(sql, 'lbl')
-            return name, start + math.prod(self._q.cardinality[d] for d in dims)
-
-        carrier = rel
+        if not restrictions:
+            if where is None:
+                rows = math.prod(self._q.cardinality[d] for d in dims)
+                cols = ', '.join(q(d) for d in dims) or q(UNIT)
+                rel = self._q.frame(dims, None)
+                sql = f'SELECT {cols}, ({self._row_major(dims, start)})::BIGINT AS {q(label)} FROM {rel.alias("p")}'
+                return self._view(sql, 'lbl'), start + rows
+
+            free = plan.free_prefix(dims, plan.predicate_dims(where, self._name_dims))
+            if free:
+                return self._factored(dims, free, where, label, start)
+
+        carrier = self._q.frame(dims, where)
         for on, presence in restrictions:
             keep = ', '.join(f'l.{q(c)}' for c in carrier.columns)
             on_sql = ' AND '.join(f'l.{q(d)} IS NOT DISTINCT FROM r.{q(d)}' for d in on)
@@ -202,17 +227,94 @@ def _label_frame(
                 f'WHERE {on_sql})',
                 carrier.columns,
             )
+        return self._counted(carrier, dims, label, start)
+
+    def _row_major(self, dims: tuple[str, ...], start: int, alias: str = '') -> str:
+        """Row-major position in *dims*' coordinate product, offset by *start*.
+
+        The trailing dim has stride 1 and every other is the product of the
+        cardinalities to its right, so the position is a dot product against
+        the ordinals the frame already carries — no ordering imposed, because
+        the answer does not depend on the order rows arrive in. The polars twin
+        is `Labeller.row_major`, and both arithmetic paths reach a label
+        through it for the reason the label itself is written once.
+        """
+        prefix = f'{alias}.' if alias else ''
+        terms: list[str] = []
+        stride = 1
+        for d in reversed(dims):
+            terms.append(f'{prefix}{q(_ordinal(d))} * {stride}')
+            stride *= self._q.cardinality[d]
+        return ' + '.join([*reversed(terms), str(start)])
+
+    def _counted(self, carrier: Rel, dims: tuple[str, ...], label: str, start: int) -> tuple[str, int]:
+        """Rank the surviving rows of *carrier*: the ordered window, and a count.
+
+        The general answer and the expensive one — a global sort of the whole
+        product, which is what the two arithmetic paths exist to avoid.
+        """
         order = ', '.join(q(_ordinal(d)) for d in dims) or '1'
         cols = ', '.join(q(d) for d in dims) or q(UNIT)
         sql = (
             f'SELECT {cols}, (ROW_NUMBER() OVER (ORDER BY {order}) - 1 + {start})::BIGINT AS {q(label)} '
             f'FROM {carrier.alias("c")}'
         )
         name = self._view(sql, 'lbl')
-        got = self._con.execute(f'SELECT count(*) FROM {q(name)}').fetchone()
+        return name, start + self._height(name)
+
+    def _factored(
+        self,
+        dims: tuple[str, ...],
+        free: int,
+        where: plan.Predicate,
+        label: str,
+        start: int,
+    ) -> tuple[str, int]:
+        """Labels for a mask that reads none of the first *free* dims.
+
+        A mask that cannot see the leading dims removes the same coordinates
+        under every one of their values, so the survivors are a rectangle: the
+        full product of the leading dims against one surviving set. Ranking
+        that set costs a window over the *set*, not over the product — on
+        `dispatch` it ranks 100 generators instead of 10M
+        ``(snapshot, generator)`` pairs, and the window is what a global sort
+        made the dominant cost of the build.
+
+        The label is then arithmetic again, through the same
+        :meth:`_row_major` the unmasked path uses: row-major over the leading
+        dims, times the width of the surviving set, plus a survivor's rank
+        within it. That is the same number the window would have counted,
+        because for each leading coordinate the same survivors appear in the
+        same order — which is what `tests/test_engine_parity.py` checks by
+        comparing the built model against the other engine's.
+        """
+        head, kept = dims[:free], dims[free:]
+        ranked = self._view(
+            f'SELECT {", ".join(q(d) for d in kept)}, '
+            f'(ROW_NUMBER() OVER (ORDER BY {", ".join(q(_ordinal(d)) for d in kept)}) - 1)::BIGINT AS "__rank" '
+            f'FROM {self._q.frame(kept, where).alias("k")}',
+            'srv',
+        )
+        width = self._height(ranked)
+        if width == 0:
+            # nothing survived anywhere, so there is no rectangle to describe.
+            # The counted path returns the right columns and dtypes for free.
+            return self._counted(self._q.frame(dims, where), dims, label, start)
+
+        product = self._q.frame(head, None)
+        picked = ', '.join([*(f'h.{q(d)}' for d in head), *(f's.{q(d)}' for d in kept)])
+        position = self._row_major(head, 0, alias='h')
+        sql = (
+            f'SELECT {picked}, (({position}) * {width} + s."__rank" + {start})::BIGINT AS {q(label)} '
+            f'FROM {product.alias("h")} CROSS JOIN {q(ranked)} AS s'
+        )
+        rows = math.prod(self._q.cardinality[d] for d in head) * width
+        return self._view(sql, 'lbl'), start + rows
+
+    def _height(self, table: str) -> int:
+        got = self._con.execute(f'SELECT count(*) FROM {q(table)}').fetchone()
         assert got is not None
-        height = got[0]
-        return name, start + height
+        return int(got[0])
 
     # -- declarations ---------------------------------------------------
 
@@ -223,11 +325,11 @@ def _build_variable(self, v: plan.VariableDeclaration) -> pl.DataFrame:
         labelled = Rel(f'SELECT * FROM {q(name)}', (*v.dims, 'var_label'))
 
         bounded = self._q.bounds(labelled, v)
-        sql = (
-            f'SELECT var_label AS col, lb::DOUBLE AS lb, ub::DOUBLE AS ub, '
-            f'{lit(v.variable_type)} AS vtype FROM {bounded.alias("b")}'
-        )
-        cols = self._fetch(sql)
+        sql = f'SELECT var_label AS col, lb::DOUBLE AS lb, ub::DOUBLE AS ub FROM {bounded.alias("b")}'
+        # `vtype` is attached here rather than selected as a literal in SQL:
+        # one word per column is one *copy* of that word per row over the wire,
+        # and the frame's stated dtype is an Enum holding four bytes.
+        cols = self._fetch(sql).with_columns(pl.lit(v.variable_type, dtype=sinks.VTYPE).alias('vtype'))
         bad = cols.filter(pl.col('lb').is_null() | pl.col('ub').is_null()).height
         if bad:
             raise DataError(null_bounds_message(v.name, bad))
@@ -286,9 +388,10 @@ def _build_constraint(self, c: plan.ConstraintDeclaration) -> tuple[pl.DataFrame
             return rows, self._fetch(stacked)
         # `sum` over `(row, col)` is the terminal aggregate — where duplicates
         # from Sum and GroupSum, which project rather than aggregate, collapse.
-        return rows, self._fetch(
-            f'SELECT row, col, sum(coeff) AS coeff FROM ({stacked}) GROUP BY row, col ORDER BY row, col'
-        )
+        # Unordered: every sink sorts the matrix into the order it needs
+        # (`lp_file` by `(row, col)`, `solver_direct` by `row`), so an ORDER BY
+        # here is a second sort of the largest frame in the model for nothing.
+        return rows, self._fetch(f'SELECT row, col, sum(coeff) AS coeff FROM ({stacked}) GROUP BY row, col')
 
     def _build_objective(self, o: plan.ObjectiveDeclaration) -> pl.DataFrame | None:
         """The objective as ``(col, coeff)``, or ``None`` if it has no terms."""
@@ -342,9 +445,7 @@ def close(self) -> None:
 def _stack(frames: list[pl.DataFrame], columns: tuple[str, ...]) -> pl.DataFrame:
     kept = [f for f in frames if f.height]
     if not kept:
-        return pl.DataFrame(
-            schema={c: (pl.String if c in ('sense', 'vtype') else pl.Float64) for c in columns}
-        ).with_columns([pl.col(c).cast(pl.Int64) for c in columns if c in ('row', 'col')])
+        return pl.DataFrame(schema={name: sinks.DTYPES[name] for name in columns})
     return pl.concat([f.select(columns) for f in kept])