exp(duck): the LP file rendered in SQL is the same file, and a worse one

bench/duck_lp.py

@@ -0,0 +1,163 @@
+"""The ``lp_file`` sink, rendered inside duckdb — #399's experiment.
+
+The experiment #399 names: build the LP text in SQL and ``COPY`` it out, so the
+model never crosses into polars on the way to a file. What the ordinary writer
+streams out of frames this streams out of tables, and the bytes must be
+identical — `tests/test_duck_lp_text.py` compares them over the whole ports
+corpus.
+
+In `bench/` rather than under `writers/` because the answer it produced was
+no: it is heavier and slower than the writer it would replace, and the one
+shape that would bound it is the one duckdb will not stream. Measurement code
+for a question already answered, kept so the answer can be re-checked.
+
+**The float rendering is the whole risk.** LP text is mostly numbers, and the
+ordinary writer spells them with polars' float cast. duckdb's differs in two
+places and nowhere else, measured over 141,840 doubles including 20,000 drawn
+from the one decade they disagree about: it pads a single-digit negative
+exponent (``1e-07`` for ``1e-7``), and it switches to scientific notation one
+decade earlier, at ``1e-5``. :func:`_number` repairs both.
+"""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+    import duckdb
+
+#: A double as the ordinary writer spells it. The two repairs are duckdb's
+#: only departures from polars' cast; every other value renders identically.
+_NUMBER = r"""
+    CASE
+        WHEN {v} = 0 OR NOT isfinite({v}) THEN CAST({v} AS VARCHAR)
+        WHEN abs({v}) >= 1e-5 AND abs({v}) < 1e-4 THEN
+            CASE WHEN {v} < 0 THEN '-0.0000' ELSE '0.0000' END ||
+            replace(split_part(CAST(abs({v}) AS VARCHAR), 'e', 1), '.', '')
+        ELSE regexp_replace(CAST({v} AS VARCHAR), 'e-0(\d)$', 'e-\1')
+    END
+"""
+
+
+def _number(value: str) -> str:
+    """*value* as LP text, matching polars' float cast."""
+    return _NUMBER.format(v=value)
+
+
+def _term(coeff: str, col: str) -> str:
+    """One ``+1.5 x7`` term, sign always explicit.
+
+    Zero is spelled rather than rendered for the reason the ordinary writer
+    spells it: ``-0.0`` takes the ``+`` arm and would give ``+-0.0``, which no
+    LP parser accepts.
+    """
+    sign = f"CASE WHEN {coeff} >= 0 THEN '+' ELSE '' END"
+    magnitude = f"CASE WHEN {coeff} = 0 THEN '0.0' ELSE {_number(coeff)} END"
+    return f"{sign} || {magnitude} || ' x' || CAST({col} AS BIGINT)::VARCHAR"
+
+
+def _bound(value: str, infinite: str) -> str:
+    return f"CASE WHEN isinf({value}) THEN '{infinite}' ELSE {_number(value)} END"
+
+
+def _line(section: int, order: str, text: str) -> str:
+    """One SELECT contributing lines to the file, keyed for the global order."""
+    return f'SELECT {section} AS s, {order}::BIGINT AS k, 0::BIGINT AS j, {text} AS line'
+
+
+def write(
+    con: duckdb.DuckDBPyConnection,
+    path: Path,
+    *,
+    col_tables: list[str],
+    obj_tables: list[str],
+    row_tables: list[str],
+    matrix_tables: list[str],
+    vtype_runs: list[tuple[str, int]],
+    sense: str,
+    constant: float,
+) -> None:
+    """Render the model at *con* as LP text at *path*, in one statement.
+
+    Every section is a branch of one union ordered by ``(section, key)``, so
+    the file is one stream duckdb sorts and writes. That sort is the point of
+    the experiment: it is the only part of a write whose size follows the
+    model, so whether the whole write stays inside ``memory_limit`` is whether
+    duckdb can spill it.
+    """
+    cols = _union(col_tables, 'SELECT lb, ub FROM {t}')
+    objective = _union(obj_tables, 'SELECT col, coeff FROM {t}')
+    rows = _union(row_tables, 'SELECT row, sense, rhs FROM {t}')
+    matrix = _union(matrix_tables, 'SELECT row, col, coeff FROM {t}')
+
+    numbered = f'SELECT *, (row_number() OVER () - 1) AS col FROM ({cols})'
+    vtype = _vtype_case(vtype_runs)
+
+    parts = [
+        _line(0, '0', f"'{'min' if sense == 'min' else 'max'}'"),
+        _line(0, '1', "''"),
+        _line(0, '2', "'obj:'"),
+    ]
+    if constant:
+        parts.append(_line(0, '3', f"printf('%+.17g', {constant!r}::DOUBLE)"))
+    parts += [
+        f'SELECT 1, col::BIGINT, 0::BIGINT, {_term("coeff", "col")} FROM ({objective})',
+        _line(2, '0', "''"),
+        _line(2, '1', "'s.t.'"),
+        _line(2, '2', "''"),
+        # a constraint's lines are keyed (row, slot): header, placeholder, one
+        # per term at its column, sense. Two integers rather than the ordinary
+        # writer's single product, which it has to bound per chunk to keep from
+        # overflowing; ordering on the pair cannot overflow at any model size.
+        f"SELECT 3, row::BIGINT, 0::BIGINT, 'c' || CAST(row AS BIGINT)::VARCHAR || ':' FROM ({rows})",
+        f"""SELECT 3, r.row::BIGINT, 1::BIGINT, '+0 x0' FROM ({rows}) AS r
+            WHERE NOT EXISTS (SELECT 1 FROM ({matrix}) AS m WHERE m.row = r.row)""",
+        f'SELECT 3, row::BIGINT, col::BIGINT + 2, {_term("coeff", "col")} FROM ({matrix})',
+        f"""SELECT 3, row::BIGINT, 9223372036854775807::BIGINT,
+                   replace(sense, '==', '=') || ' ' || {_number('rhs')} FROM ({rows})""",
+        _line(4, '0', "''"),
+        _line(4, '1', "'bounds'"),
+        f"""SELECT 5, col::BIGINT, 0::BIGINT,
+                   {_bound('lb', '-infinity')} || ' <= x' || CAST(col AS BIGINT)::VARCHAR
+                   || ' <= ' || {_bound('ub', '+infinity')} FROM ({numbered})""",
+    ]
+    for section, kind, keyword in ((6, 'binary', 'binary'), (7, 'integer', 'general')):
+        if not any(t == kind for t, _ in vtype_runs):
+            continue
+        parts += [
+            _line(section, '-2', "''"),
+            _line(section, '-1', f"'{keyword}'"),
+            f"SELECT {section}, col::BIGINT, 0::BIGINT, 'x' || CAST(col AS BIGINT)::VARCHAR "
+            f"FROM ({numbered}) WHERE {vtype} = '{kind}'",
+        ]
+    parts += [_line(8, '0', "''"), _line(8, '1', "'end'")]
+
+    union = ' UNION ALL '.join(f'({p})' for p in parts)
+    con.execute(
+        f"COPY (SELECT line FROM ({union}) ORDER BY s, k, j) TO '{path}' "
+        "(FORMAT CSV, HEADER false, QUOTE '', ESCAPE '', DELIMITER '')"
+    )
+
+
+def _union(tables: list[str], template: str) -> str:
+    """The declarations' shares as one relation, or an empty one."""
+    if not tables:
+        return template.format(t='(SELECT NULL) AS empty') + ' WHERE false'
+    return ' UNION ALL '.join(template.format(t=f'"{t}"') for t in tables)
+
+
+def _vtype_case(runs: list[tuple[str, int]]) -> str:
+    """A column's variable type, from the run lengths rather than a stored word.
+
+    The engine holds one word per *declaration*, not per column, for the reason
+    `cols` states: the frame's dtype is an Enum and the runs already say how
+    long each is. Rebuilt here as a boundary over ``col`` so the word still
+    never exists once per row.
+    """
+    at, arms = 0, []
+    for kind, height in runs:
+        arms.append(f"WHEN col < {at + height} THEN '{kind}'")
+        at += height
+    return 'CASE ' + ' '.join(arms) + " ELSE 'continuous' END" if arms else "'continuous'"

tests/conftest.py

@@ -89,6 +89,21 @@ def engine_internals(pytestconfig: pytest.Config) -> None:
         pytest.skip(f'reaches polars-engine internals; this run is on {name!r}')
 
 
+@pytest.fixture
+def duck_internals(pytestconfig: pytest.Config) -> None:
+    """Skip unless this run is on the duckdb engine.
+
+    The mirror of :func:`engine_internals`, for the few tests that reach inside
+    *that* engine. Same reasoning: an engine does not owe another engine its
+    internals, only the same model.
+    """
+    from lpspec.relational import engines
+
+    name = pytestconfig.getoption('--engine') or engines.DEFAULT_ENGINE
+    if name != 'duckdb':
+        pytest.skip(f'reaches duckdb-engine internals; this run is on {name!r}')
+
+
 @pytest.fixture(autouse=True, scope='session')
 def _engine_under_test(pytestconfig: pytest.Config) -> Iterator[None]:
     """Point `lps.build`'s default at `--engine`, for the whole session.

tests/test_duck_lp_text.py

@@ -0,0 +1,96 @@
+"""The experiment's evidence: LP text rendered in SQL is the same LP text.
+
+#399 asks whether the LP file can be written from inside duckdb, so the model
+never crosses into polars on the way out. Whether that is *worth* doing is a
+measurement and lives in the PR. Whether it is even the same file is this, and
+it has to hold before any measurement means anything.
+"""
+
+from __future__ import annotations
+
+import json
+import sys
+from pathlib import Path
+from typing import Any
+
+import polars as pl
+import pytest
+
+import lpspec as lps
+from lpspec.relational.engines.duck.executor import DuckExecutor
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+from bench import duck_lp
+
+PORTS = Path(__file__).resolve().parent.parent / 'examples' / 'ports'
+REFERENCES: dict[str, dict[str, Any]] = json.loads((PORTS / 'references.json').read_text())
+
+
+def _sources(name: str) -> dict[str, Any]:
+    data = json.loads((PORTS / 'data' / f'{name}.json').read_text())
+    return {k: pl.DataFrame(v) if isinstance(v, dict) else v for k, v in data.items()}
+
+
+def _render_in_sql(ex: DuckExecutor, path: Path) -> None:
+    tables = ex._tables()
+    duck_lp.write(
+        ex._con,
+        path,
+        col_tables=list(ex._col_tables),
+        obj_tables=list(ex._obj_tables),
+        row_tables=list(ex._row_shares),
+        matrix_tables=[name for name, _, _ in ex._matrix_spans],
+        vtype_runs=[(vtype, height) for _, vtype, height in ex._col_runs],
+        sense=tables.objective_sense,
+        constant=tables.objective_constant,
+    )
+
+
+@pytest.mark.parametrize('name', sorted(REFERENCES))
+def test_the_sql_writer_renders_the_same_bytes(name: str, tmp_path: Path, duck_internals: None) -> None:
+    """Every port, both writers, compared as bytes.
+
+    The corpus is the point: it carries the sections a small model does not —
+    ``binary`` and ``general`` for the four MILPs, termless rows, infinite
+    bounds, an objective constant, and coefficients across enough magnitudes to
+    reach both places duckdb's float rendering departs from polars'.
+    """
+    program = str(PORTS / f'{name}.yaml')
+    with lps.build(program, _sources(name)) as ex:
+        assert isinstance(ex, DuckExecutor), 'this run is not on the duckdb engine'
+        ex.write(tmp_path / 'ordinary.lp')
+        _render_in_sql(ex, tmp_path / 'sql.lp')
+    ordinary = (tmp_path / 'ordinary.lp').read_bytes()
+    assert (tmp_path / 'sql.lp').read_bytes() == ordinary, f'{name}: the SQL writer rendered different bytes'
+
+
+@pytest.mark.parametrize(
+    'value',
+    [
+        pytest.param(1e-7, id='padded-negative-exponent'),
+        pytest.param(2.5e-8, id='padded-negative-exponent-with-mantissa'),
+        pytest.param(1.2345678901234568e-5, id='the-decade-polars-keeps-positional'),
+        pytest.param(-9.99e-5, id='the-same-decade-negative'),
+        pytest.param(1e16, id='the-upper-switch-both-agree-on'),
+        pytest.param(-0.0, id='negative-zero'),
+    ],
+)
+def test_the_number_repairs_are_the_ones_that_matter(value: float, tmp_path: Path, duck_internals: None) -> None:
+    """The two departures :mod:`duck_lp` repairs, each as its own case.
+
+    A coefficient is where they reach the file, so the check is the rendered
+    LP rather than the SQL: a repair that stopped working would show up here as
+    a term a solver reads as a different number.
+    """
+    model = {
+        'dimensions': {'i': {'dtype': 'int', 'values': [0]}},
+        'parameters': {'c': {'dims': ['i']}},
+        'variables': {'x': {'foreach': ['i'], 'bounds': {'lower': 0, 'upper': 1}}},
+        'constraints': {'k': {'foreach': ['i'], 'expression': 'x * c >= 0'}},
+        'objectives': {'o': {'sense': 'minimize', 'expression': 'sum(x, over=i)'}},
+    }
+    sources = {'c': pl.DataFrame({'i': [0], 'value': [value]})}
+    with lps.build(model, sources) as ex:
+        ex.write(tmp_path / 'ordinary.lp')
+        _render_in_sql(ex, tmp_path / 'sql.lp')
+    assert (tmp_path / 'sql.lp').read_bytes() == (tmp_path / 'ordinary.lp').read_bytes()