Fix binary encoding: use ceil(log2(n_k)) qubits per bucket instead of ceil(log2(n_k+1))

[Copilot]

The previous implementation allocated ceil(log2(n_k + 1)) qubits per bucket, implicitly reserving an index for a "null route". Since the problem structure guarantees exactly one route is selected per bucket, the correct formula per docs/main_hubo.tex §5.2 (Eq. 10) is ceil(log2(n_k)) — no null slot needed.

Core fix

• Qubit count (__init__): max(1, ceil(log2(n_k))) — saves 1 qubit for any bucket where n_k is a power of 2
• _decode: idx = v % n_k; route = bucket_columns[k][idx] — direct modular map into [0, n_k), no null-route guard
• create_circuit: same v % n_k decode in both cost_vals and cov_vals build loops

Test updates

For the 7-column test problem (bucket 0: 3 routes, bucket 1: 4 routes):

• N_qubits: 4 → was 5 (bucket 1 now needs ceil(log2(4))=2 qubits, not 3)
• Bitstrings for known feasible solutions updated: "0010" (c0+c4), "1001" (c1+c5), "0111" (c2+c6)
• Infeasible test case replaced: "0000" decodes to c0+c3 (a real but order-violating route pair), replacing the former all-zeros "null route" which no longer has meaning

Original prompt

Implement BucketExactCover problem class with IO support

Mathematical specification

The full mathematical derivation is in docs/main_hubo.tex on the exact-cover-hubo branch. Read this file carefully before implementing. Key sections:

• §5.1: Partitioning columns into boat-wise groups (bucket definition)
• §5.2: Binary encoding of route selection — use Strategy B: Modular wrapping for overflow indices
• §5.3: Constructing the HUBO (multilinear polynomial interpolation, cost function, order coverage polynomials, order constraint penalties, penalization factor, scaling §5.3.6)
• §5.4: Mapping to a Hamiltonian and the QAOA circuit (substitution q_m = (1 - Z_m)/2, multi-qubit Z-rotation via CNOT-ladder + RZ)
• §5.5: QUBO vs HUBO comparison (Table 2)

Existing code to read first

• qaoa/problems/base_problem.py — Problem base class (inherit from this directly)
• qaoa/problems/qubo_problem.py — QUBO class (reference for circuit construction style)
• qaoa/problems/exactcover_problem.py — ExactCover class (the drop-in replacement target; mirror its interface)
• qaoa/utils/qaoaIO.py — IO infrastructure to extend
• unittests/test_problems.py — existing test patterns to follow

---

New file: qaoa/problems/bucketexactcover_problem.py

Implement class BucketExactCover(Problem).

Constructor signature (note: num_buckets is required, not optional):

def __init__(
    self,
    columns,              # np.ndarray shape (N_rows, N_routes)
    num_buckets,          # int: number of top rows used for bucket assignment (REQUIRED)
    weights=None,         # np.ndarray shape (N_routes,) or None (defaults to zeros)
    penalty_factor=None,  # float or None (auto-computed if None)
    allow_infeasible=False,
    scale_problem=False,
) -> None

Bucket extraction:

• columns[:num_buckets, :] is the boat-assignment block.
• Column r belongs to bucket k if columns[k, r] == 1.
• Each column belongs to exactly one bucket (exactly one 1 in its top num_buckets entries).
• Columns that belong to no bucket (all zeros in the top num_buckets rows) are ignored — they cannot be selected in any feasible solution.
• Store the list of valid column indices per bucket in self._bucket_columns[k].

Binary encoding per bucket:

• Bucket k has n_k = len(self._bucket_columns[k]) columns.
• Assign b_k = ceil(log2(n_k + 1)) qubits (minimum 1 qubit even if n_k == 1).
• A b_k-bit binary string z encodes integer v = sum_j 2^j * z_j.
• Map to column selection via Strategy B (modular wrapping): idx = v % (n_k + 1). If idx == 0 → no column selected (null route). If idx ∈ [1, n_k] → global column index self._bucket_columns[k][idx - 1] is selected.
• Total qubits: N_qubits = sum_k b_k. Qubits are laid out bucket-by-bucket: first b_0 qubits for bucket 0, then b_1 for bucket 1, etc.

cost(string) -> float:

• Decode the bitstring → per-bucket selections → route indicator vector x (length = number of valid columns, indexed in order of self._valid_columns).
• Compute: -(weights_used @ x + penalty_factor * sum((Y_orders @ x - 1)**2)) where Y_orders = columns[num_buckets:, :][:, valid_column_indices].
• When scale_problem=True, weights and penalty are already scaled in self.weights and self.penalty_factor.
• The cost convention matches ExactCover.cost(): negated (QAOA maximizes).

unscaled_cost(string) -> float:

• Same as cost() but always uses self._original_weights and self._original_penalty_factor.

isFeasible(string) -> bool:

• Decode → x. Return True iff Y_orders @ x == ones(N_orders) (use np.allclose with atol=1e-7), or if self.allow_infeasible.

create_circuit():

1. Build the HUBO polynomial using multilinear interpolation (Eq. 20–22 in docs/main_hubo.tex). For each bucket k and each function value slot (all 2^b_k binary strings of length b_k), evaluate the cost and penalty contributions at that point using the modular wrapping decoding.
2. Collect all polynomial terms by expanding products of qubit variables.
3. Substitute q_m = (1 - Z_m)/2 for each qubit (Eq. 33). Collect Pauli-Z monomials and their coefficients beta_T.
4. Skip terms where abs(beta_T) < 1e-10.
5. For each non-zero monomial over qubit index set T of degree d:
   • If d == 0: skip (constant term, no gate needed).
   • If d == 1: add RZ(2 * gamma * beta_T, qubit).
   • If d > 1: implement using CNOT-ladder + RZ + reverse CNOT-ladder. Specifically: CNOT from T[0] to T[1], CNOT from T[1] to T[2], ..., CNOT from T[d-2] to T[d-1], then RZ(2 * gamma * beta_T, T[d-1]), then reverse CNOTs.
6. Use a single Qiskit Parameter("x_gamma") for gamma (same name as in QUBO.create_circuit()).
7. The circuit returned is only the cost unitary (phase separator). Initial state (|+>^n) and mixer are handled by the QAOA framework.
8. Store the circuit in self.circuit and retur...

This pull request was created from Copilot chat.

---

📱 Kick off Copilot coding agent tasks wherever you are with GitHub Mobile, available on iOS and Android.