Metadata-Version: 2.1
Name: kd-simulator
Version: 0.2.1
Description-Content-Type: text/markdown
Requires-Dist: numpy
Requires-Dist: matplotlib

# Qubit Kirkwood-Dirac Simulator

A Python package for simulating quantum circuits with qubits using the Kirkwood-Dirac (KD) quasiprobability distribution. This simulator implements a hidden variable model that efficiently simulates stabilizer circuits and reveals when quantum computation requires genuine quantum resources.

## Overview

This package implements a classical simulation framework for qubit quantum circuits based on the Kirkwood-Dirac quasiprobability representation. The KD distribution provides a phase-space-like representation of quantum states that becomes a true probability distribution for certain quantum states, enabling efficient classical simulation.

The simulator supports:
- Preparation of arbitrary single-qubit states
- Efficient sampling from KD distributions
- Clifford gates (Hadamard, Pauli X, Pauli Z, CNOT)
- Computational basis measurements
- Fast vectorized operations for statistical simulations
- Multiple visualization tools for measurement outcomes
- Detection of quantum advantage through KD nonpositivity

## Installation

```bash
pip install kd_simulator
```

## Requirements

- Python 3.8+
- NumPy
- Matplotlib

## Quick Start

```python
import numpy as np
from KD_Simulator import (
    make_qubit_rhos,
    compute_per_qubit_kd,
    sample_from_per_qubit_kd,
    Hadamard_vectorized,
    CNOT_vectorized,
    Zmeasurement_vectorized,
    plot_top_k
)

# Initialize random number generator
rng = np.random.default_rng()

# Prepare initial states: two qubits in |0+⟩ state
state_labels = ["0", "+"]
rho_list = make_qubit_rhos(state_labels)

# Compute KD distributions
kd_matrices = compute_per_qubit_kd(rho_list)

# Sample hidden variables
num_samples = 10000
g_all, chi_all = sample_from_per_qubit_kd(kd_matrices, num_samples, rng)

# Apply gates to create Bell state
g_all, chi_all = Hadamard_vectorized(g_all, chi_all)
g_all, chi_all = CNOT_vectorized(g_all, chi_all, 0, 1, num_qubits=2)

# Measure both qubits
outcomes_0, g_all, chi_all = Zmeasurement_vectorized(g_all, chi_all, 0, rng)
outcomes_1, g_all, chi_all = Zmeasurement_vectorized(g_all, chi_all, 1, rng)

# Combine outcomes and count
output_all = np.column_stack([outcomes_0, outcomes_1])
counts = np.bincount([int(''.join(map(str, row)), 2) for row in output_all], 
                      minlength=4)

# Visualize results
plot_top_k(counts, num_qubits=2, k=4)
```

## Theoretical Background

### Kirkwood-Dirac Quasiprobability Distribution

The Kirkwood-Dirac distribution is a quasiprobability representation of quantum states defined with respect to two incompatible observables. For qubits, we consider computational basis states {|g⟩} and Hadamard-basis states {|χ⟩}, where g, χ ∈ {0, 1}.

For a single-qubit density matrix ρ, the KD distribution is:

```
Q(g, χ) = [(ρH†) ∘ Hᵀ]_{g,χ}
```

where:
- H is the Hadamard matrix: H = (1/√2)[[1, 1], [1, -1]]
- ∘ denotes element-wise (Hadamard) product
- The result is a real-valued 2×2 matrix representing K(g, χ)

### Key Properties

- **Normalization**: ΣK(g, χ) = 1 (marginals sum to 1)
- **Quasiprobability**: K(g, χ) can be negative for quantum states
- **KD-positive states**: When K(g, χ) ≥ 0 for all (g, χ), the state is "classical" and admits efficient simulation
- **Resource for quantum advantage**: KD nonpositivity (negative values) is a necessary resource for quantum computational advantage

### Classical Simulation via Hidden Variables

When all KD values are non-negative, K(g, χ) can be interpreted as a classical joint probability distribution over the binary variables (g, χ). The simulation workflow:

1. **Initialization**: Compute K(g, χ) for each qubit's initial state
2. **Sampling**: Draw (g, χ) pairs according to K as a probability distribution
3. **Gate evolution**: Update (g, χ) deterministically via linear transformations (mod 2)
4. **Measurement**: Extract outcome from g with stochastic update to χ

This approach efficiently simulates circuits where KD distributions remain non-negative throughout, revealing when quantum resources are genuinely required.

## Core Functions

### State Preparation

**`make_qubit_rhos(state_labels)`**

Creates density matrices for common single-qubit states.

- **Parameters:**
  - `state_labels`: List of state labels, one per qubit
    - `"0"`: |0⟩ state
    - `"1"`: |1⟩ state
    - `"+"`: |+⟩ state (equal superposition)
    - `"-"`: |−⟩ state (X-basis)
    - `"r"`: |R⟩ state (right circular)
    - `"l"`: |L⟩ state (left circular)
- **Returns:** List of density matrices (2×2 complex arrays)

**`compute_per_qubit_kd(rho_list)`**

Computes the Kirkwood-Dirac distribution for each qubit independently.

- **Parameters:**
  - `rho_list`: List of single-qubit density matrices
- **Returns:** List of 2×2 KD matrices (real-valued)
- **Raises:** `ValueError` if the state has negative KD values (non-simulatable)

This function verifies that:
- The KD matrix is real (checks imaginary components)
- All KD values are non-negative (required for classical simulation)

### Hidden Variable Sampling

**`sample_from_per_qubit_kd(kd_matrices, num_samples, rng)`**

Samples hidden variables (g, χ) for all qubits across many runs.

- **Parameters:**
  - `kd_matrices`: List of KD matrices from `compute_per_qubit_kd()`
  - `num_samples`: Number of sampling runs
  - `rng`: NumPy random number generator (`np.random.default_rng()`)
- **Returns:** Tuple (g_all, chi_all), each of shape (num_samples, num_qubits)

Each qubit's (g, χ) is sampled independently from its marginal KD distribution.

### Quantum Gates

All gates operate on batched hidden variables and return updated copies. Gates transform (g, χ) via deterministic linear operations over ℤ₂.

**Single-Qubit Gates:**

- `Hadamard_vectorized(g_all, chi_all)`: Hadamard on all qubits (swaps g and χ)
- `PauliX_vectorized(g_all, chi_all, qubit)`: Pauli X gate (bit flip on g)
- `PauliZ_vectorized(g_all, chi_all, qubit)`: Pauli Z gate (bit flip on χ)

**Two-Qubit Gates:**

- `CNOT_vectorized(g_all, chi_all, control_qubit, target_qubit, num_qubits)`: Controlled-NOT gate

**Parameters:**
- `g_all, chi_all`: Hidden variable arrays (num_samples, num_qubits)
- `qubit`: Target qubit index
- `control_qubit, target_qubit`: Control and target indices for CNOT
- `num_qubits`: Total number of qubits in the circuit

**Returns:** Updated (g_all, chi_all) arrays

### Measurement

**`Zmeasurement_vectorized(g_all, chi_all, qubit, rng)`**

Performs computational basis (Z-basis) measurement on a specified qubit.

- **Parameters:**
  - `g_all, chi_all`: Hidden variable arrays (num_samples, num_qubits)
  - `qubit`: Index of qubit to measure
  - `rng`: NumPy random number generator
- **Returns:** Tuple (outcomes, g_all, chi_all)
  - `outcomes`: Measurement results (num_samples,) - values are 0 or 1
  - `g_all`: Unchanged g values
  - `chi_all`: Updated χ values (stochastically flipped with probability 0.5)

The measurement outcome is deterministically given by g[qubit], with a random post-measurement update to χ[qubit].

### Visualization Functions

**`plot_top_k(counts, num_qubits, k=20)`**

Displays the k most probable measurement outcomes as a bar chart.

**`plot_hamming_weight(counts, num_qubits)`**

Shows the distribution of measurement outcomes grouped by Hamming weight (number of 1s).

**`plot_single_qubit_marginals(counts, num_qubits)`**

Visualizes the marginal probability distribution for each individual qubit as a heatmap.

**`plot_sample_heatmap(output_all)`**

Creates a heatmap showing sampled bitstrings across shots and qubits.

**`plot_qubit_correlation(output_all)`**

Displays the correlation matrix between all pairs of qubits.

**Common Parameters:**
- `counts`: Array of length 2^num_qubits with outcome counts
- `output_all`: Array of shape (num_samples, num_qubits) with measurement outcomes
- `num_qubits`: Number of qubits in the circuit
- `k`: Number of top outcomes to display

<!-- ## Advanced Usage -->

<!-- ### Creating and Verifying Bell States

```python
import numpy as np
from qubit_kd_simulator import *

rng = np.random.default_rng(seed=42)

# Prepare |00⟩
rho_list = make_qubit_rhos(["0", "0"])
kd_matrices = compute_per_qubit_kd(rho_list)

num_samples = 50000
g_all, chi_all = sample_from_per_qubit_kd(kd_matrices, num_samples, rng)

# Create Bell state |Φ+⟩ = (|00⟩ + |11⟩)/√2
g_all, chi_all = Hadamard_vectorized(g_all, chi_all)
g_all, chi_all = CNOT_vectorized(g_all, chi_all, 0, 1, num_qubits=2)

# Measure both qubits
outcomes_0, g_all, chi_all = Zmeasurement_vectorized(g_all, chi_all, 0, rng)
outcomes_1, g_all, chi_all = Zmeasurement_vectorized(g_all, chi_all, 1, rng)

output_all = np.column_stack([outcomes_0, outcomes_1])

# Verify perfect correlation
plot_qubit_correlation(output_all)
# Should show strong correlation between qubits 0 and 1

# Check outcome distribution
counts = np.bincount([int(''.join(map(str, row)), 2) for row in output_all], 
                      minlength=4)
print(f"P(00) = {counts[0]/num_samples:.3f}")
print(f"P(11) = {counts[3]/num_samples:.3f}")
# Should each be approximately 0.5
``` -->

<!-- ### Multi-Qubit GHZ State

```python
# Prepare n qubits in |0⟩
n_qubits = 5
rho_list = make_qubit_rhos(["0"] * n_qubits)
kd_matrices = compute_per_qubit_kd(rho_list)

g_all, chi_all = sample_from_per_qubit_kd(kd_matrices, 10000, rng)

# Create GHZ state: (|00...0⟩ + |11...1⟩)/√2
g_all, chi_all = Hadamard_vectorized(g_all, chi_all)
for i in range(n_qubits - 1):
    g_all, chi_all = CNOT_vectorized(g_all, chi_all, 0, i+1, n_qubits)

# Measure all qubits
outcomes = []
for q in range(n_qubits):
    result, g_all, chi_all = Zmeasurement_vectorized(g_all, chi_all, q, rng)
    outcomes.append(result)

output_all = np.column_stack(outcomes)
counts = np.bincount([sum(row * (2**i) for i, bit in enumerate(row)) 
                       for row in output_all], minlength=2**n_qubits)

# Visualize GHZ correlations
plot_hamming_weight(counts, n_qubits)
# Should show peaks at Hamming weight 0 and n_qubits
```

### Testing Non-Classical States

```python
# Attempt to simulate a non-stabilizer state
# T gate creates a "magic state" with negative KD values

theta = np.pi / 4
cos_t = np.cos(theta / 2)
sin_t = np.sin(theta / 2)

# |ψ⟩ = cos(θ/2)|0⟩ + e^(iπ/4)sin(θ/2)|1⟩
psi = np.array([cos_t, np.exp(1j * np.pi / 4) * sin_t])
rho_magic = np.outer(psi, psi.conj())

try:
    kd = compute_per_qubit_kd([rho_magic])
    print("State is KD-positive - can be simulated classically")
except ValueError as e:
    print(f"State has negative KD values - requires quantum resources!")
    print(f"Error: {e}")
```

### Checking State Simulatability

```python
# Analyze which states are KD-positive (classically simulatable)
test_states = {
    "|0⟩": make_qubit_rhos(["0"])[0],
    "|+⟩": make_qubit_rhos(["+"])[0],
    "|R⟩": make_qubit_rhos(["r"])[0],
}

for name, rho in test_states.items():
    try:
        kd = compute_per_qubit_kd([rho])[0]
        min_val = np.min(kd)
        print(f"{name}: KD-positive (min = {min_val:.6f}) ✓ simulatable")
    except ValueError:
        print(f"{name}: Has negative KD values ✗ not simulatable")
``` -->

## Performance Considerations

- **Vectorized operations**: All gate operations are fully vectorized over samples for maximum efficiency
- **Memory scaling**: O(num_samples × num_qubits) for hidden variable storage
- **Time complexity**: O(num_gates × num_samples × num_qubits) for circuit simulation
- **Optimal sample sizes**: 10,000-100,000 samples provide good statistical accuracy
- **Qubit scaling**: Efficient for 10-20 qubits with reasonable sample counts

## Limitations

- **KD-positive states only**: Cannot simulate states with negative KD values (magic states, T-gate outputs)
- **Clifford gates only**: Supports H, X, Z, CNOT; does not support T gate or arbitrary rotations
- **Computational basis measurements**: Only Z-basis measurements are directly supported
- **Product state initialization**: Initial states must be tensor products of single-qubit states
- **Positivity requirement**: The simulator will raise an error if it encounters negative KD values

## Mathematical Details

### KD Distribution Formula

For a single-qubit density matrix ρ, the Kirkwood-Dirac distribution is:

```
K(g, χ) = ⟨g|ρ|χ⟩⟨χ|g⟩
```

In the computational basis {|0⟩, |1⟩} and Hadamard basis {|+⟩, |−⟩}, this becomes:

```
K = (ρH†) ∘ Hᵀ
```

where ∘ is element-wise multiplication. The result is a 2×2 real matrix satisfying ΣK(g,χ) = 1.

### Gate Update Rules (mod 2 arithmetic)

**Hadamard:** (g, χ) → (χ, g)

**Pauli X on qubit q:** gq → gq ⊕ 1 (flip g bit)

**Pauli Z on qubit q:** χq → χq ⊕ 1 (flip χ bit)

**CNOT from control c to target t:**
- g' = g · Act where Act = I + |t⟩⟨c|
- χ' = χ · Bct where Bct = I + |c⟩⟨t|

All operations are linear transformations over ℤ₂ (binary field).

### Measurement Protocol

Measuring qubit q in the computational basis:
1. **Outcome**: mq = gq (deterministic)
2. **Post-measurement update**: With probability 1/2, flip χq

This reflects the contextual nature of quantum measurements in the KD framework.

## Quantum Advantage and KD Nonpositivity

This simulator implements results from recent research showing that KD nonpositivity is a necessary resource for quantum computational advantage:

- **Stabilizer circuits** preserve KD-positivity and can be efficiently simulated classically
- **Magic states** (e.g., T-gate outputs) introduce negative KD values
- **Quantum speedup** requires operations that produce negative KD values
<!-- - **Resource monotone**: KD nonpositivity cannot increase under free operations -->

When the simulator encounters negative KD values, it signals that the quantum circuit requires genuine quantum resources beyond classical hidden variable models.

## Relation to Other Frameworks

- **Wigner functions**: The KD distribution is analogous to Wigner functions but defined over discrete phase space for qubits
- **Stabilizer formalism**: Stabilizer states are precisely the KD-positive states
<!-- - **Resource theories**: KD nonpositivity is a monotone for quantum computational resources -->
<!-- - **Contextuality**: Negative KD values indicate quantum contextuality -->

## Contributing

Contributions are welcome! Please open an issue to discuss changes or submit a Pull Request.


## Citation

If you use this simulator in your research, please cite:

```bibtex
@software{qubit_kd_simulator,
  title={Qubit Kirkwood-Dirac Simulator},
  author={Rishi Goel},
  year={2025},
  url={https://github.com/rishigoel2003/KD_Simulator}
}
```

And please cite the foundational paper on KD distributions as a computational resource:

```bibtex
@article{thio2025kirkwood,
  title={Kirkwood-Dirac Nonpositivity is a Necessary Resource for Quantum Computing},
  author={Thio, Jonathan J and Yang, Songqinghao and De Bi{\`e}vre, Stephan and Barnes, Crispin HW and Arvidsson-Shukur, David RM},
  journal={arXiv preprint arXiv:2506.08092},
  year={2025}
}
```


## Contact

For questions and support, please open an issue on GitHub or contact [rishigoel25@gmail.com].
