"""
SDP SOS Certificate Solver — Python I/O shim.

Formulates sum-of-squares (SOS) polynomial optimization problems as
semidefinite programs (SDP), calls an external solver, rationalizes
the floating-point solution to exact rationals, and emits:
  1. A Lean 4 data file (GoormaghtighCert.lean) with concrete coefficients
  2. A JSON receipt for the audit trail

Architecture (per AGENTS.md programming choice flow):
  - Python owns I/O: calling solver, formatting output
  - Lean owns decisions: verifying the certificate (SDPVerify.lean)
  - Float at the solver boundary ONLY, immediately rationalized

Solver: CVXPY + SCS (open-source). Falls back to stub mode if not installed.

TODO(lean-port): The polynomial formulation logic in _build_moment_matrix
  could eventually move to Lean once a Lean SDP interface exists.

Usage:
  python3 sdp_sos_solver.py --problem test_simple --degree 2
  python3 sdp_sos_solver.py --problem goormaghtigh --degree 4 --level 1 \\
      --output-lean path/to/GoormaghtighCert.lean \\
      --output-json path/to/sdp_certificate.json
"""

from __future__ import annotations

import argparse
import hashlib
import json
import math
import os
import sys
import time
from fractions import Fraction
from itertools import combinations_with_replacement
from typing import Any

# --- Lazy import: solver only needed at the API boundary ---
try:
    import cvxpy as cp
    import numpy as np
    HAS_CVXPY = True
except ImportError:
    HAS_CVXPY = False

# --- Q16_16 constants (per AGENTS.md) ---
Q16_SCALE = 65536  # 2^16 = 0x00010000


# ================================================================
# §0  MONOMIAL BASIS GENERATION
# ================================================================

def monomial_basis(n_vars: int, max_degree: int) -> list[tuple[int, ...]]:
    """Generate all monomials in n_vars variables up to max_degree.

    Returns list of exponent tuples, sorted lexicographically.
    Example: monomial_basis(2, 2) = [(0,0), (0,1), (0,2), (1,0), (1,1), (2,0)]
    """
    basis: list[tuple[int, ...]] = []
    # Generate all exponent vectors with sum <= max_degree
    _generate_monomials(n_vars, max_degree, 0, (), basis)
    basis.sort()
    return basis


def _generate_monomials(
    n_vars: int, remaining_deg: int, var_idx: int,
    current: tuple[int, ...], result: list[tuple[int, ...]]
) -> None:
    """Recursively generate monomial exponent vectors."""
    if var_idx == n_vars:
        result.append(current)
        return
    for d in range(remaining_deg + 1):
        _generate_monomials(
            n_vars, remaining_deg - d, var_idx + 1,
            current + (d,), result
        )


# ================================================================
# §1  SPARSE POLYNOMIAL REPRESENTATION
# ================================================================

class SparsePoly:
    """Sparse polynomial with Fraction coefficients.

    Internally: dict mapping exponent tuples to Fraction coefficients.
    """
    def __init__(self, terms: dict[tuple[int, ...], Fraction] | None = None):
        self.terms: dict[tuple[int, ...], Fraction] = terms or {}

    @staticmethod
    def from_list(n_vars: int, data: list[tuple[Fraction, list[int]]]) -> "SparsePoly":
        """Build from [(coeff, [exponents])] list."""
        terms: dict[tuple[int, ...], Fraction] = {}
        for coeff, expon in data:
            key = tuple(expon + [0] * (n_vars - len(expon)))
            terms[key] = terms.get(key, Fraction(0)) + coeff
        return SparsePoly({k: v for k, v in terms.items() if v != 0})

    def __add__(self, other: "SparsePoly") -> "SparsePoly":
        result = dict(self.terms)
        for k, v in other.terms.items():
            result[k] = result.get(k, Fraction(0)) + v
        return SparsePoly({k: v for k, v in result.items() if v != 0})

    def __mul__(self, other: "SparsePoly") -> "SparsePoly":
        result: dict[tuple[int, ...], Fraction] = {}
        for ka, va in self.terms.items():
            for kb, vb in other.terms.items():
                key = tuple(a + b for a, b in zip(ka, kb))
                result[key] = result.get(key, Fraction(0)) + va * vb
        return SparsePoly({k: v for k, v in result.items() if v != 0})

    def square(self) -> "SparsePoly":
        return self * self

    def scale(self, c: Fraction) -> "SparsePoly":
        return SparsePoly({k: c * v for k, v in self.terms.items() if c * v != 0})

    def to_lean_list(self) -> list[tuple[str, list[int]]]:
        """Convert to Lean-compatible [(coeff_string, exponents)] format."""
        result = []
        for expon, coeff in sorted(self.terms.items()):
            if coeff != 0:
                result.append((fraction_to_lean(coeff), list(expon)))
        return result


# ================================================================
# §2  PROBLEM FORMULATION
# ================================================================

def build_test_simple() -> dict[str, Any]:
    """Test problem: verify x₀² + x₁² ≥ 0 (trivial SOS)."""
    return {
        "name": "test_simple",
        "n_vars": 2,
        "target": SparsePoly.from_list(2, [
            (Fraction(1), [2, 0]),
            (Fraction(1), [0, 2]),
        ]),
        "constraints": [],
        "description": "x₀² + x₁² ≥ 0 (trivial SOS)",
    }


def build_test_weighted() -> dict[str, Any]:
    """Test problem: verify x₀² + x₀ ≥ 0 on {x₀ ≥ 0}."""
    return {
        "name": "test_weighted",
        "n_vars": 1,
        "target": SparsePoly.from_list(1, [
            (Fraction(1), [2]),
            (Fraction(1), [1]),
        ]),
        "constraints": [
            SparsePoly.from_list(1, [(Fraction(1), [1])]),  # g₀ = x₀
        ],
        "description": "x₀² + x₀ ≥ 0 on {x₀ ≥ 0}",
    }


def build_goormaghtigh_fixed(m: int = 3, n: int = 3) -> dict[str, Any]:
    """Goormaghtigh collision polynomial for fixed (m, n).

    For fixed m, n:
      R(x, m) = 1 + x + x² + ... + x^(m-1)
      R(y, n) = 1 + y + y² + ... + y^(n-1)
      p = (R(x,m) - R(y,n))²

    Variables: x = v₀, y = v₁ (2 variables for fixed m,n).
    Domain: K = {x ≥ 91, y ≥ 2}.
    """
    # Build R(x, m) = Σ_{k=0}^{m-1} x^k
    rx_terms: dict[tuple[int, ...], Fraction] = {}
    for k in range(m):
        rx_terms[(k, 0)] = Fraction(1)
    rx = SparsePoly(rx_terms)

    # Build R(y, n) = Σ_{k=0}^{n-1} y^k
    ry_terms: dict[tuple[int, ...], Fraction] = {}
    for k in range(n):
        ry_terms[(0, k)] = Fraction(1)
    ry = SparsePoly(ry_terms)

    # p = (R(x,m) - R(y,n))²
    diff = rx + ry.scale(Fraction(-1))
    target = diff.square()

    # Constraints: x ≥ 91, y ≥ 2
    constraints = [
        SparsePoly.from_list(2, [(Fraction(1), [1, 0]), (Fraction(-91), [0, 0])]),
        SparsePoly.from_list(2, [(Fraction(1), [0, 1]), (Fraction(-2), [0, 0])]),
    ]

    return {
        "name": f"goormaghtigh_m{m}_n{n}",
        "n_vars": 2,
        "target": target,
        "constraints": constraints,
        "description": (
            f"Goormaghtigh collision (R(x,{m}) - R(y,{n}))² ≥ 0 "
            f"on {{x ≥ 91, y ≥ 2}}"
        ),
    }


# ================================================================
# §3  SDP SOLVER INTERFACE
# ================================================================

def solve_sos_sdp(
    problem: dict[str, Any],
    degree: int,
    level: int = 0,
    solver: str = "SCS",
) -> dict[str, Any]:
    """Formulate and solve the SOS SDP problem.

    Args:
        problem: Problem dict from build_* functions.
        degree: Maximum degree for SOS components (half the polynomial degree).
        level: Putinar hierarchy level (0 = pure SOS, 1+ = with constraints).
        solver: CVXPY solver name ("SCS", "MOSEK", "SDPA").

    Returns:
        Dict with:
        - sos_components: List[SparsePoly]  (rationalized)
        - weighted_pairs: List[(SparsePoly, SparsePoly)]
        - status: str
        - solve_time: float
        - solver: str
    """
    if not HAS_CVXPY:
        print("WARNING: cvxpy not installed. Using stub certificate.", file=sys.stderr)
        return _stub_solve(problem, degree, level)

    target = problem["target"]
    n_vars = problem["n_vars"]
    constraints_polys = problem["constraints"]

    # Monomial basis for SOS components of the given degree
    basis = monomial_basis(n_vars, degree)
    basis_size = len(basis)

    # All monomials that appear in the expansion (up to 2*degree)
    full_basis = monomial_basis(n_vars, 2 * degree)

    # Build the moment matrix variable Q (PSD)
    Q = cp.Variable((basis_size, basis_size), symmetric=True)

    # Constraint: Q ≽ 0 (positive semidefinite)
    psd_constraints = [Q >> 0]

    # Build coefficient matching constraints
    # The SOS polynomial is: p_sos(x) = basis(x)ᵀ Q basis(x)
    # = Σᵢⱼ Qᵢⱼ · x^(αᵢ + αⱼ)
    # We need: coeff of x^γ in p_sos = coeff of x^γ in target

    coeff_constraints = []
    for gamma in full_basis:
        # Sum all Qᵢⱼ where αᵢ + αⱼ = γ
        lhs_terms = []
        for i, alpha_i in enumerate(basis):
            for j, alpha_j in enumerate(basis):
                if tuple(a + b for a, b in zip(alpha_i, alpha_j)) == gamma:
                    lhs_terms.append((i, j))

        if not lhs_terms:
            # This monomial doesn't appear in the SOS expansion
            # It should also not appear in the target
            target_coeff = float(target.terms.get(gamma, Fraction(0)))
            if abs(target_coeff) > 1e-12:
                print(f"WARNING: monomial {gamma} in target but not in SOS "
                      f"basis (degree too low?)", file=sys.stderr)
            continue

        lhs = sum(Q[i, j] for i, j in lhs_terms)
        target_coeff = float(target.terms.get(gamma, Fraction(0)))
        coeff_constraints.append(lhs == target_coeff)

    # If level > 0, add weighted SOS components for each constraint
    weighted_Qs = []
    if level > 0 and constraints_polys:
        weighted_basis = monomial_basis(n_vars, degree - 1)
        w_basis_size = len(weighted_basis)
        for g_poly in constraints_polys:
            Qw = cp.Variable((w_basis_size, w_basis_size), symmetric=True)
            psd_constraints.append(Qw >> 0)
            weighted_Qs.append(Qw)
            # TODO(lean-port): The coefficient matching for weighted terms
            # is more complex — need to convolve with constraint polynomial

    # Solve
    prob = cp.Problem(cp.Minimize(0), psd_constraints + coeff_constraints)
    t0 = time.time()
    try:
        prob.solve(solver=solver, verbose=False)
    except cp.SolverError as e:
        return {
            "sos_components": [],
            "weighted_pairs": [],
            "status": f"SOLVER_ERROR: {e}",
            "solve_time": time.time() - t0,
            "solver": solver,
        }
    solve_time = time.time() - t0

    if prob.status not in ("optimal", "optimal_inaccurate"):
        return {
            "sos_components": [],
            "weighted_pairs": [],
            "status": prob.status,
            "solve_time": solve_time,
            "solver": solver,
        }

    # Extract SOS components from Q via eigendecomposition
    Q_val = Q.value
    sos_components = _extract_sos_from_gram(Q_val, basis, n_vars)

    # Extract weighted components
    weighted_pairs = []
    for idx, Qw in enumerate(weighted_Qs):
        Qw_val = Qw.value
        weighted_basis = monomial_basis(n_vars, degree - 1)
        w_components = _extract_sos_from_gram(Qw_val, weighted_basis, n_vars)
        # For each weighted SOS component sⱼ, pair with constraint gⱼ
        for s_comp in w_components:
            weighted_pairs.append((s_comp, constraints_polys[idx]))

    return {
        "sos_components": sos_components,
        "weighted_pairs": weighted_pairs,
        "status": prob.status,
        "solve_time": solve_time,
        "solver": solver,
    }


def _extract_sos_from_gram(
    Q: "np.ndarray",
    basis: list[tuple[int, ...]],
    n_vars: int,
    tol: float = 1e-8,
) -> list[SparsePoly]:
    """Extract SOS components from a Gram matrix via eigendecomposition.

    Q = Σ λᵢ vᵢvᵢᵀ where λᵢ ≥ 0 (PSD).
    Each eigenvector vᵢ with λᵢ > tol gives an SOS component:
      qᵢ(x) = √λᵢ · Σⱼ vᵢⱼ · x^αⱼ

    Returns list of SparsePoly (rationalized coefficients).
    """
    eigenvalues, eigenvectors = np.linalg.eigh(Q)
    components = []

    for k in range(len(eigenvalues)):
        lam = eigenvalues[k]
        if lam < tol:
            continue
        sqrt_lam = math.sqrt(float(lam))
        v = eigenvectors[:, k]

        # Build the polynomial: qₖ(x) = √λₖ · Σⱼ vⱼ · x^αⱼ
        terms: dict[tuple[int, ...], Fraction] = {}
        for j, alpha in enumerate(basis):
            coeff_float = sqrt_lam * float(v[j])
            if abs(coeff_float) < tol:
                continue
            # Float → rational at the boundary (immediately rationalized)
            coeff_rat = Fraction(coeff_float).limit_denominator(10**8)
            key = tuple(alpha)
            terms[key] = terms.get(key, Fraction(0)) + coeff_rat

        if terms:
            components.append(SparsePoly(
                {k: v for k, v in terms.items() if v != 0}
            ))

    return components


def _stub_solve(
    problem: dict[str, Any], degree: int, level: int
) -> dict[str, Any]:
    """Stub solver for when CVXPY is not installed.

    Returns a hand-constructed certificate for known test problems.
    """
    name = problem["name"]
    n_vars = problem["n_vars"]

    if name == "test_simple":
        # x₀² + x₁² = x₀² + x₁²  (trivial: q₀ = x₀, q₁ = x₁)
        return {
            "sos_components": [
                SparsePoly.from_list(2, [(Fraction(1), [1, 0])]),
                SparsePoly.from_list(2, [(Fraction(1), [0, 1])]),
            ],
            "weighted_pairs": [],
            "status": "stub_optimal",
            "solve_time": 0.0,
            "solver": "stub",
        }
    elif name == "test_weighted":
        # x₀² + x₀ on {x₀ ≥ 0}: sos=[x₀], weighted=[(1, x₀)]
        return {
            "sos_components": [
                SparsePoly.from_list(1, [(Fraction(1), [1])]),
            ],
            "weighted_pairs": [
                (SparsePoly.from_list(1, [(Fraction(1), [0])]),  # s₀ = 1
                 SparsePoly.from_list(1, [(Fraction(1), [1])])), # g₀ = x₀
            ],
            "status": "stub_optimal",
            "solve_time": 0.0,
            "solver": "stub",
        }
    else:
        return {
            "sos_components": [],
            "weighted_pairs": [],
            "status": f"stub_unsupported:{name}",
            "solve_time": 0.0,
            "solver": "stub",
        }


# ================================================================
# §4  RATIONALIZATION
# ================================================================

def fraction_to_lean(f: Fraction) -> str:
    """Convert a Fraction to a Lean ℚ literal string.

    Examples:
      Fraction(3, 1)  → "3"
      Fraction(1, 2)  → "(1/2)"
      Fraction(-5, 3) → "(-5/3)"
    """
    if f.denominator == 1:
        return str(f.numerator)
    else:
        return f"({f.numerator}/{f.denominator})"


def rationalize_float(x: float, max_denom: int = 10**8) -> Fraction:
    """Convert a float to the closest rational with bounded denominator.

    This is the BOUNDARY CONVERSION: float (solver output) → ℚ (exact).
    Per AGENTS.md: Float at external boundary only, immediately bracketed.
    """
    return Fraction(x).limit_denominator(max_denom)


# ================================================================
# §5  LEAN CERTIFICATE EMITTER
# ================================================================

def emit_lean_certificate(
    result: dict[str, Any],
    problem: dict[str, Any],
    degree: int,
    level: int,
    output_path: str,
) -> str:
    """Generate a Lean 4 data file with the SOS certificate.

    Returns the SHA256 of the generated file.
    """
    n_vars = problem["n_vars"]
    target = problem["target"]
    sos_components = result["sos_components"]
    weighted_pairs = result["weighted_pairs"]

    lines = [
        "/-",
        "  Machine-generated SOS certificate.",
        f"  Problem: {problem['name']}",
        f"  Description: {problem['description']}",
        f"  Solver: {result['solver']}",
        f"  Status: {result['status']}",
        f"  Solve time: {result['solve_time']:.4f}s",
        f"  Generated: {time.strftime('%Y-%m-%dT%H:%M:%SZ', time.gmtime())}",
        "",
        "  DO NOT EDIT — regenerate with:",
        f"    python3 sdp_sos_solver.py --problem {problem['name']} "
        f"--degree {degree} --level {level}",
        "-/",
        "",
        "import Semantics.SDPVerify",
        "",
        "open Semantics.SDPVerify",
        "",
        f"namespace Semantics.SDPCert.{_lean_name(problem['name'])}",
        "",
    ]

    # Emit SOS components
    lines.append(f"/-- SOS certificate for: {problem['description']} -/")
    lines.append(f"def certificate : SDPCertificate {n_vars} :=")
    lines.append(f"  mkCertificate {n_vars}")

    # SOS components
    lines.append("    -- SOS components (qᵢ polynomials)")
    lines.append("    [" + ",\n     ".join(
        _poly_to_lean(q, n_vars) for q in sos_components
    ) + "]")

    # Weighted pairs
    lines.append("    -- Weighted pairs (sⱼ, gⱼ)")
    if weighted_pairs:
        lines.append("    [" + ",\n     ".join(
            f"({_poly_to_lean(s, n_vars)}, {_poly_to_lean(g, n_vars)})"
            for s, g in weighted_pairs
        ) + "]")
    else:
        lines.append("    []")

    # Target polynomial
    lines.append("    -- Target polynomial p")
    lines.append("    " + _poly_to_lean_flat(target, n_vars))

    # Degree and level
    lines.append(f"    -- degree = {degree}, level = {level}")
    lines.append(f"    {degree} {level}")
    lines.append("")

    # Verification witness
    lines.append("/-- #eval witness: certificate is valid. -/")
    lines.append("#eval verifyCertificate certificate  -- Expected: true")
    lines.append("")
    lines.append(f"end Semantics.SDPCert.{_lean_name(problem['name'])}")
    lines.append("")

    content = "\n".join(lines)

    os.makedirs(os.path.dirname(output_path) or ".", exist_ok=True)
    with open(output_path, "w") as f:
        f.write(content)

    sha = hashlib.sha256(content.encode()).hexdigest()
    print(f"Wrote Lean certificate to {output_path} (SHA256: {sha})",
          file=sys.stderr)
    return sha


def _lean_name(name: str) -> str:
    """Convert a problem name to a Lean-valid namespace component."""
    return "".join(c if c.isalnum() else "_" for c in name).capitalize()


def _poly_to_lean(p: SparsePoly, n_vars: int) -> str:
    """Convert a SparsePoly to a Lean list literal."""
    terms = p.to_lean_list()
    if not terms:
        return "[]"
    parts = []
    for coeff_str, expon in terms:
        parts.append(f"({coeff_str}, {list(expon)})")
    return "[" + ", ".join(parts) + "]"


def _poly_to_lean_flat(p: SparsePoly, n_vars: int) -> str:
    """Convert a SparsePoly to a Lean list literal on a single line."""
    return _poly_to_lean(p, n_vars)


# ================================================================
# §6  JSON RECEIPT EMITTER
# ================================================================

def emit_json_receipt(
    result: dict[str, Any],
    problem: dict[str, Any],
    degree: int,
    level: int,
    lean_sha: str,
    output_path: str,
) -> None:
    """Write a machine-readable JSON receipt for the SDP computation."""
    receipt = {
        "schema": "sdp_sos_certificate_v1",
        "problem": {
            "name": problem["name"],
            "description": problem["description"],
            "n_vars": problem["n_vars"],
            "n_target_terms": len(problem["target"].terms),
            "n_constraints": len(problem.get("constraints", [])),
        },
        "solver": {
            "name": result["solver"],
            "status": result["status"],
            "solve_time_seconds": result["solve_time"],
        },
        "certificate": {
            "n_sos_components": len(result["sos_components"]),
            "n_weighted_pairs": len(result["weighted_pairs"]),
            "degree": degree,
            "putinar_level": level,
            "lean_file_sha256": lean_sha,
        },
        "claim_boundary": (
            "sdp-computed;python-rationalized;lean-verified"
        ),
        "generated_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
    }

    os.makedirs(os.path.dirname(output_path) or ".", exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(receipt, f, indent=2)
    print(f"Wrote JSON receipt to {output_path}", file=sys.stderr)


# ================================================================
# §7  CLI ENTRY POINT
# ================================================================

PROBLEMS = {
    "test_simple": build_test_simple,
    "test_weighted": build_test_weighted,
}


def main() -> None:
    parser = argparse.ArgumentParser(
        description="SDP SOS Certificate Solver",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=__doc__,
    )
    parser.add_argument(
        "--problem", required=True,
        help="Problem name (test_simple, test_weighted, goormaghtigh_mM_nN)",
    )
    parser.add_argument(
        "--degree", type=int, default=2,
        help="Maximum degree for SOS components (default: 2)",
    )
    parser.add_argument(
        "--level", type=int, default=0,
        help="Putinar hierarchy level (default: 0 = pure SOS)",
    )
    parser.add_argument(
        "--solver", default="SCS",
        help="CVXPY solver name (default: SCS)",
    )
    parser.add_argument(
        "--output-lean",
        help="Path to write the Lean certificate file",
    )
    parser.add_argument(
        "--output-json",
        help="Path to write the JSON receipt file",
    )
    parser.add_argument(
        "--m", type=int, default=3,
        help="Goormaghtigh parameter m (default: 3)",
    )
    parser.add_argument(
        "--n", type=int, default=3,
        help="Goormaghtigh parameter n (default: 3)",
    )

    args = parser.parse_args()

    # Build the problem
    if args.problem in PROBLEMS:
        problem = PROBLEMS[args.problem]()
    elif args.problem.startswith("goormaghtigh"):
        problem = build_goormaghtigh_fixed(args.m, args.n)
    else:
        print(f"ERROR: Unknown problem '{args.problem}'", file=sys.stderr)
        print(f"Available: {list(PROBLEMS.keys())} + goormaghtigh",
              file=sys.stderr)
        sys.exit(1)

    print(f"Problem: {problem['name']}", file=sys.stderr)
    print(f"  {problem['description']}", file=sys.stderr)
    print(f"  {len(problem['target'].terms)} target terms", file=sys.stderr)
    print(f"  {len(problem.get('constraints', []))} constraints",
          file=sys.stderr)
    print(f"  Degree: {args.degree}, Level: {args.level}", file=sys.stderr)

    # Solve
    result = solve_sos_sdp(
        problem, args.degree, args.level, args.solver
    )

    print(f"Status: {result['status']}", file=sys.stderr)
    print(f"Solve time: {result['solve_time']:.4f}s", file=sys.stderr)
    print(f"SOS components: {len(result['sos_components'])}", file=sys.stderr)
    print(f"Weighted pairs: {len(result['weighted_pairs'])}", file=sys.stderr)

    # Verify locally (Python-side sanity check)
    if result["sos_components"]:
        recon = SparsePoly()
        for q in result["sos_components"]:
            recon = recon + q.square()
        for s, g in result["weighted_pairs"]:
            recon = recon + (s * g)

        # Compare with target
        diff = recon + problem["target"].scale(Fraction(-1))
        max_residual = max(
            (abs(v) for v in diff.terms.values()), default=Fraction(0)
        )
        print(f"Max residual: {float(max_residual):.2e}", file=sys.stderr)

    # Emit outputs
    lean_sha = ""
    if args.output_lean:
        lean_sha = emit_lean_certificate(
            result, problem, args.degree, args.level, args.output_lean
        )
    if args.output_json:
        emit_json_receipt(
            result, problem, args.degree, args.level,
            lean_sha, args.output_json
        )

    # Summary
    if result["status"] in ("optimal", "optimal_inaccurate", "stub_optimal"):
        print("SUCCESS: SOS certificate computed.", file=sys.stderr)
        sys.exit(0)
    else:
        print(f"FAILED: solver status = {result['status']}", file=sys.stderr)
        sys.exit(1)


if __name__ == "__main__":
    main()