qlbm D1Q3 examples

[richhl]

We are researching using qlbm. Would like to simulate 1D Burgers flow using D1Q3. Do you have any example of building 1d lattice with a set of 3 discrete speeds and periodic boundary conditions?

Thanks in advance. We have found examples, unit tests and tutorials for 2d, 3d but its being hard figuring out how to do D1Q3.

[gcalin]

Thank you for your interest.

Which encoding are you exploring? If you are interested in the linear encoding used in QLGA (noncompressed, one qubit per velocity channel), you can use the example in the LQLGA simulation notebook as a guideline.

If you are interested in the amplitude-based encoding (logarithmically many grid qubits + logarithmically many velocity qubits), we do not yet have examples for that in the regular DdQq discretizations or 1D. Is this something you are specifically interested in? If so, let me know.

[richhl]

Thx for the hints. About the encoding I've started with LQLGA but... is it possible to encode a sinusoidal initial state of velocities? If so, how?

Maybe I need weights for encoding velocities and this lead to amplitude-based encoding. am I wrong? If so, then we have to use amplitude based encondng and i'm interested in it :)

[gcalin]

Had to think about this one for a little bit, as it is not how I typically think about the linear encoding. But it is indeed possible. We don't have any initial conditions like this the code yet, but I wrote some (messy) code that uses Ry gates to encode the sin wave into the \ket{1} state of each velocity channel. Hope this is useful!

import math

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister, transpile
from qiskit_aer import AerSimulator

sns.set_theme()

def piecewise_sin(Ng: int, k: int = 1) -> np.ndarray:
    centers = [(b + 0.5) for b in range(Ng)]
    band_vals = [math.sin(2 * math.pi * k * xc / Ng) for xc in centers]
    v = np.zeros(Ng)
    for x in range(Ng):
        b = min(int((x * Ng) // Ng), Ng - 1)
        v[x] = band_vals[b]
    return v

def ry_angle_from_prob(p: float) -> float:
    p = min(1.0, max(0.0, p))
    return 2.0 * math.asin(math.sqrt(p))

def d1q3_qubit_index(site: int, channel: str) -> int:
    base = 3 * site
    if channel == 'rest':
        return base + 0
    if channel == 'right':
        return base + 1
    if channel == 'left':
        return base + 2

def marginal_probs_from_counts(counts: dict[str, int], nbits: int) -> np.ndarray:
    total = sum(counts.values())
    probs = np.zeros(nbits, dtype=float)
    for bitstring, c in counts.items():
        s = bitstring.replace(' ', '')  # just in case
        s_rev = s[::-1]
        for i, ch in enumerate(s_rev):
            if ch == '1':
                probs[i] += c
    return probs / max(1, total)

def build_d1q3_density_sinusoid_circuit(
    Ng: int = 5,
    k: int = 1,
    rho_base: float = 0.30,
    rho_amp: float = 0.25,
    frac_rest: float = 0.20,
):

    # piecewise sinusoid in [-1,1]
    sin_pw = piecewise_sin(Ng, k=k)

    # target per-site density
    rho = np.clip(rho_base + rho_amp * sin_pw, 0.0, 1.0)

    # Zero-momentum split across channels
    frac_pm = (1.0 - frac_rest) / 2.0
    p_rest = frac_rest * rho
    p_right = frac_pm * rho
    p_left = frac_pm * rho

    # Build circuiit with 3*Ng qubits
    qreg = QuantumRegister(3 * Ng, "q")
    creg = ClassicalRegister(3 * Ng, "c")
    qc = QuantumCircuit(qreg, creg, name="D1Q3_piecewise_density")

    # Apply per-qubit Ry to set P(\ket{1}) = p_channel(site)
    for x in range(Ng):
        for ch, p in (("rest", p_rest[x]), ("right", p_right[x]), ("left", p_left[x])):
            qi = d1q3_qubit_index(x, ch)
            theta = ry_angle_from_prob(float(p))
            if abs(theta) > 1e-12:
                qc.ry(theta, qreg[qi])

    # Measure all qubits (for sampling the initial distribution)
    qc.measure(qreg, creg)

    metadata = {
        "Ng": Ng,
        "k": k,
        "rho_target": rho,
        "p_rest_target": p_rest,
        "p_right_target": p_right,
        "p_left_target": p_left,
    }
    return qc, metadata

def run_and_plot(qc: QuantumCircuit, meta: dict, shots: int):
    sim = AerSimulator(method="statevector")
    tqc = transpile(qc, sim, optimization_level=0)
    res = sim.run(tqc, shots=shots).result()
    counts = res.get_counts(0)

    nbits = qc.num_qubits
    probs_1 = marginal_probs_from_counts(counts, nbits)

    Ng = meta["Ng"]

    # extract per-site, per-channel from marginals
    p_rest_meas = np.array([probs_1[d1q3_qubit_index(x, "rest")] for x in range(Ng)])
    p_right_meas = np.array([probs_1[d1q3_qubit_index(x, "right")] for x in range(Ng)])
    p_left_meas  = np.array([probs_1[d1q3_qubit_index(x, "left")]  for x in range(Ng)])
    rho_meas = p_rest_meas + p_right_meas + p_left_meas

    # targets
    rho_t = meta["rho_target"]
    p0_t = meta["p_rest_target"]
    pr_t = meta["p_right_target"]
    pl_t = meta["p_left_target"]

    x = np.arange(Ng)

    # Density Plot
    plt.figure(figsize=(7, 4.2))
    plt.plot(x, rho_meas, marker='o', linewidth=1.5, label="Measured ρ(x)")
    plt.plot(x, rho_t, marker='x', linestyle='--', linewidth=1.0, label="Target ρ(x)")
    plt.title("D1Q3 piecewise-sinusoid density (total)")
    plt.xlabel("Lattice site x")
    plt.ylabel("Density")
    plt.legend()
    plt.tight_layout()
    plt.show()

    # Per Channel Plot
    plt.figure(figsize=(7, 4.2))
    plt.plot(x, p_rest_meas, marker='o', linewidth=1.2, label="Measured p_rest")
    plt.plot(x, p0_t, marker='x', linestyle='--', linewidth=1.0, label="Target p_rest")

    plt.plot(x, p_right_meas, marker='o', linewidth=1.2, label="Measured p_right")
    plt.plot(x, pr_t, marker='x', linestyle='--', linewidth=1.0, label="Target p_right")

    plt.plot(x, p_left_meas, marker='o', linewidth=1.2, label="Measured p_left")
    plt.plot(x, pl_t, marker='x', linestyle='--', linewidth=1.0, label="Target p_left")

    plt.title("D1Q3 piecewise-sinusoid per-channel occupancies")
    plt.xlabel("Lattice site x")
    plt.ylabel("Occupancy probability")
    plt.legend()
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    qc, meta = build_d1q3_density_sinusoid_circuit(
        Ng=5,
        k=1,
        rho_base=0.30,
        rho_amp=0.25,
        frac_rest=0.20,
    )

    print(qc.name, "with", qc.num_qubits, "qubits and", qc.num_clbits, "clbits")
    print("meta:", {k: meta[k] for k in ("Ng","k")})

    run_and_plot(qc, meta, shots=2**16)