QUBO

qubo_utils.py

from typing import Dict, Tuple, List, Optional

import numpy as np
import dimod

def choose_penalties(S: np.ndarray, Pu: Optional[float], Pc: Optional[float]) -> Tuple[float, float]:
    """
    Heuristic penalty chooser: scale penalties to dominate typical |S|.
    If user-specified, pass through.
    """
    sscale = float(np.max(np.abs(S))) if S.size else 1.0
    base = max(1.0, sscale)
    # conservative: a bit larger than base to ensure constraints hold
    Pu_val = float(Pu) if Pu is not None else 2.0 * base
    Pc_val = float(Pc) if Pc is not None else 2.0 * base
    return Pu_val, Pc_val

def build_qubo_naive(S: np.ndarray, C: List[int], Pu: Optional[float]=None, Pc: Optional[float]=None):
    """
    Build a naive QUBO dict (i,j)->Qij as in the user's definition (one-hot per row, capacity per column).
    Returns (Qdict, varmap) with varmap[k] = (i,a).
    """
    S = np.asarray(S, dtype=float)
    C = np.asarray(C, dtype=int)
    N, M = S.shape
    assert C.size == M
    mask = np.ones_like(S, dtype=bool)

    Pu_val, Pc_val = choose_penalties(S, Pu, Pc)

    varmap: List[Tuple[int,int]] = []
    idx_of = -np.ones_like(S, dtype=int)
    for i in range(N):
        for a in range(M):
            if mask[i, a]:
                idx = len(varmap)
                idx_of[i, a] = idx
                varmap.append((i, a))
    nvar = len(varmap)
    if nvar == 0:
        raise ValueError("No variables after mask.")

    Qdict: Dict[Tuple[int,int], float] = {}

    def addQ(p: int, q: int, w: float):
        key = (p, q) if p <= q else (q, p)
        Qdict[key] = Qdict.get(key, 0.0) + float(w)

    # Objective: diagonal -S_{i,a}
    for k, (i, a) in enumerate(varmap):
        addQ(k, k, -float(S[i, a]))

    # Row one-hot penalties: (sum_a x_{i,a} - 1)^2
    for i in range(N):
        cols = [a for a in range(M) if mask[i, a]]
        for a in cols:
            k = idx_of[i, a]
            addQ(k, k, -Pu_val)  # linear part
        for u in range(len(cols)):
            k1 = idx_of[i, cols[u]]
            for v in range(u+1, len(cols)):
                k2 = idx_of[i, cols[v]]
                addQ(k1, k2, 2.0 * Pu_val)

    # Column capacities: (sum_i x_{i,a} - C[a])^2
    for a in range(M):
        rows = [i for i in range(N) if mask[i, a]]
        lin = Pc_val * (1.0 - 2.0 * float(C[a]))
        for i in rows:
            k = idx_of[i, a]
            addQ(k, k, lin)
        for u in range(len(rows)):
            k1 = idx_of[rows[u], a]
            for v in range(u+1, len(rows)):
                k2 = idx_of[rows[v], a]
                addQ(k1, k2, 2.0 * Pc_val)

    return Qdict, varmap, Pu_val, Pc_val

def qubo_to_bqm(Qdict: Dict[Tuple[int,int], float]) -> dimod.BinaryQuadraticModel:
    return dimod.BinaryQuadraticModel.from_qubo(Qdict)

def bqm_problem_graph(bqm: dimod.BinaryQuadraticModel):
    """Return networkx.Graph of the quadratic interactions (edges where bias != 0)."""
    import networkx as nx
    G = nx.Graph()
    for u in bqm.variables:
        G.add_node(int(u))
    for (u, v), bias in bqm.quadratic.items():
        if bias != 0.0:
            G.add_edge(int(u), int(v))
    return G

run_qubo_embedding_minimize.py

import argparse, os, json, numpy as np, pandas as pd
import networkx as nx

import dimod

from qubo_utils import build_qubo_naive, qubo_to_bqm, bqm_problem_graph
from run_embedding_minimize_qubits import build_hardware_graph, find_best_embedding

try:
    from dwave.preprocessing import fix_variables
    HAVE_PRE = True
except Exception:
    HAVE_PRE = False

def load_matrix(path: str) -> np.ndarray:
    ext = os.path.splitext(path)[1].lower()
    if ext in [".npy", ".npz"]:
        arr = np.load(path)
        if isinstance(arr, np.lib.npyio.NpzFile):
            # pick first array
            arr = arr[list(arr.files)[0]]
        return np.asarray(arr, dtype=float)
    elif ext in [".csv", ".tsv"]:
        df = pd.read_csv(path) if ext==".csv" else pd.read_csv(path, sep="\t")
        return df.values.astype(float)
    else:
        raise ValueError("Unsupported format for S. Use .npy/.npz/.csv/.tsv")

def load_capacities(path: str) -> np.ndarray:
    ext = os.path.splitext(path)[1].lower()
    if ext in [".npy", ".npz"]:
        arr = np.load(path)
        if isinstance(arr, np.lib.npyio.NpzFile):
            arr = arr[list(arr.files)[0]]
        return np.asarray(arr, dtype=int)
    elif ext in [".csv", ".tsv", ".txt"]:
        data = np.loadtxt(path, dtype=int, delimiter="," if ext==".csv" else None)
        return np.asarray(data, dtype=int).ravel()
    else:
        raise ValueError("Unsupported format for C. Use .npy/.npz/.csv/.tsv/.txt")

def main():
    ap = argparse.ArgumentParser(description="QUBO -> BQM -> problem graph -> embedding qubit minimization")
    ap.add_argument("--S", type=str, required=True, help="Path to S matrix (.npy/.npz/.csv/.tsv)")
    ap.add_argument("--C", type=str, required=True, help="Path to capacities vector (.npy/.npz/.csv/.tsv/.txt)")
    ap.add_argument("--Pu", type=float, default=None, help="Row (one-hot) penalty; if omitted, choose automatically")
    ap.add_argument("--Pc", type=float, default=None, help="Column (capacity) penalty; if omitted, choose automatically")
    ap.add_argument("--hardware", choices=["pegasus", "zephyr"], default="pegasus")
    ap.add_argument("--scale", type=int, default=16)
    ap.add_argument("--attempts", type=int, default=64)
    ap.add_argument("--patience-grid", type=int, nargs="+", default=[5,10,15])
    ap.add_argument("--seed", type=int, default=13)
    ap.add_argument("--init", choices=["layout", "none"], default="layout")
    ap.add_argument("--layout", choices=["spring", "spectral", "kamada_kawai"], default="spring")
    ap.add_argument("--refine-topk", type=float, default=0.1)
    ap.add_argument("--time-budget-sec", type=float, default=180)
    ap.add_argument("--presolve", action="store_true", help="Apply variable fixing preprocessing (roof duality etc.)")
    ap.add_argument("--out", type=str, required=True)

    args = ap.parse_args()
    os.makedirs(args.out, exist_ok=True)

    # 1) Build QUBO and BQM
    S = load_matrix(args.S)
    C = load_capacities(args.C)
    Qdict, varmap, Pu_val, Pc_val = build_qubo_naive(S, C, Pu=args.Pu, Pc=args.Pc)
    bqm = qubo_to_bqm(Qdict)

    # 2) Optional presolve
    presolve_info = {}
    if args.presolve and HAVE_PRE:
        fixed, bqm_reduced = fix_variables(bqm)
        presolve_info = {
            "fixed_count": int(len(fixed)),
            "remaining_variables": int(len(bqm_reduced.variables))
        }
        bqm = bqm_reduced
    elif args.presolve and not HAVE_PRE:
        presolve_info = {"note": "dwave-preprocessing not installed at runtime"}

    # 3) Problem graph
    H = bqm_problem_graph(bqm)
    nx.write_gpickle(H, os.path.join(args.out, "problem_graph.gpickle"))

    # 4) Build hardware
    T = build_hardware_graph(args.hardware, args.scale)

    # 5) Embedding minimize (multi-start + refine)
    best_emb, summary, df = find_best_embedding(
        H, T,
        attempts=args.attempts,
        patience_grid=tuple(args.patience_grid),
        init=args.init,
        layout_kind=args.layout,
        seed=args.seed,
        topk_refine_ratio=args.refine_topk,
        time_budget_sec=args.time_budget_sec
    )

    # 6) Save artifacts
    df.to_csv(os.path.join(args.out, "attempts_summary.csv"), index=False)
    with open(os.path.join(args.out, "best_embedding.json"), "w") as f:
        json.dump({str(int(k)): list(map(int, v)) for k, v in best_emb.items()}, f, indent=2)
    meta = {
        "Pu": Pu_val, "Pc": Pc_val,
        "presolve": presolve_info,
        "best_metrics": summary.get("best_metrics"),
        "attempts": int(summary.get("attempts", 0)),
        "successes": int(summary.get("successes", 0))
    }
    with open(os.path.join(args.out, "meta.json"), "w") as f:
        json.dump(meta, f, indent=2)

    print("Done. Best:", summary.get("best_metrics"))
    print("Artifacts in:", args.out)

if __name__ == "__main__":
    main()

QUBO-specific: minimize qubits for this problem

Given your QUBO definition (rows: one-hot; cols: capacity), place your S (N×M) and C (length M) on disk and run:

# S.npy: float matrix; C.npy: int vector (capacities)
python src/run_qubo_embedding_minimize.py \
  --S S.npy --C C.npy \
  --hardware pegasus --scale 16 \
  --attempts 64 --patience-grid 5 10 15 \
  --init layout --layout spring \
  --refine-topk 0.1 \
  --presolve \
  --time-budget-sec 180 \
  --out out/my_qubo_embed