Struct DensityMatrix 

pub struct DensityMatrix {
    pub num_qubits: usize,
    /* private fields */
}

Quantum density matrix representing mixed or pure quantum states.

A density matrix describes the statistical state of an N-qubit quantum system. Unlike a statevector which can only represent pure states, a density matrix can represent mixed states (ensembles of pure states).

Memory Layout

The data vector uses a contiguous memory layout representing a flattened $2^N \times 2^N$ matrix. To optimize simulation performance, the simulator employs a $2N$-qubit isomorphism:

• The matrix is treated as a statevector of $2N$ qubits.
• The “ket” side (Left U) acts on the upper $N$ qubits (indices $N$ to $2N-1$).
• The “bra” side (Right $U^\dagger$) acts on the lower $N$ qubits (indices 0 to $N-1$).

Example

use cqlib_core::qis::DensityMatrix;

// Create a 1-qubit density matrix in state |0><0|
let mut dm = DensityMatrix::new(1);
assert_eq!(dm.num_qubits, 1);

// Apply Hadamard gate -> |+><+|
dm.apply_h(0);

// Probabilities should be 0.5 for both |0> and |1>
let probs = dm.probabilities();
assert!((probs[0] - 0.5).abs() < 1e-10);
assert!((probs[1] - 0.5).abs() < 1e-10);

Fields

num_qubits: usize

Number of qubits in the system ($N$).

Implementations

impl DensityMatrix

pub fn data(&self) -> &[Complex64]

Returns the flattened row-major density matrix data.

pub fn new(num_qubits: usize) -> Self

Creates a new density matrix initialized to the pure state $|0\dots 0\rangle\langle 0\dots 0|$.

Arguments

• num_qubits - Number of qubits in the system.

Panics

Panics if the dimension is unrepresentable or allocation fails. Use Self::try_new to receive an error instead.

pub fn try_new(num_qubits: usize) -> Result<Self, QisError>

Creates the ground state, returning an error for an unrepresentable dimension or a failed allocation.

Errors

Returns QisError::InvalidParameterValue for an unrepresentable dimension, or QisError::UnsupportedOperation if allocation fails.

pub fn maximally_mixed(num_qubits: usize) -> Self

Creates the maximally mixed state $I / 2^N$.

The returned density matrix has equal probability on every computational-basis state and zero off-diagonal entries. A zero-qubit system is represented by the one-dimensional density matrix [1].

Like every dense density-matrix constructor, this allocates $O(4^N)$ complex values.

Panics

Panics if the dimension is unrepresentable or allocation fails. Use Self::try_maximally_mixed to receive an error instead.

pub fn try_maximally_mixed(num_qubits: usize) -> Result<Self, QisError>

Fallible variant of Self::maximally_mixed.

Errors

Returns QisError::InvalidParameterValue for an unrepresentable dimension, or QisError::UnsupportedOperation if allocation fails.

pub fn from_state( num_qubits: usize, initial_state: Vec<Complex64>, ) -> Result<Self, QisError>

Creates a density matrix from an initial statevector (pure state).

Internally computes the outer product $\rho = |\psi\rangle\langle\psi|$.

Arguments

• num_qubits - Number of qubits in the system.
• initial_state - Vector of $2^N$ complex amplitudes representing a normalized pure state.

Errors

Returns QisError::InvalidStateDimension if the initial_state length is incorrect, or QisError::NotNormalized if it is not normalized.

pub fn from_density_matrix_state( num_qubits: usize, dm_state: Vec<Complex64>, ) -> Result<Self, QisError>

Creates a density matrix directly from a flattened $2^N \times 2^N$ matrix.

Validates all physical constraints: Hermiticity, positive semidefiniteness, and unit trace.

Arguments

• num_qubits - Number of qubits in the system.
• dm_state - Vector of $4^N$ complex values representing the density matrix.

Errors

Returns QisError::InvalidStateDimension if the matrix length is incorrect. Returns QisError::NotHermitian if the matrix is not Hermitian. Returns QisError::NotPositiveSemidefinite if the matrix has negative eigenvalues. Returns QisError::NotNormalized if the trace is not equal to 1.

pub fn probabilities(&self) -> Vec<f64>

Computes the measurement probability distribution over all computational basis states.

Extracts the diagonal elements of the density matrix, which represent the probabilities $P(|i\rangle) = \rho_{ii}$.

Returns

Vector of probabilities summing to 1.0.

pub fn from_circuit(circuit: &Circuit) -> Result<Self, QisError>

Constructs a density matrix by simulating a quantum circuit.

Executes the circuit gates sequentially to evolve the initial $|0\dots 0\rangle\langle 0\dots 0|$ state.

Arguments

• circuit - The quantum circuit to simulate.

Returns

• Ok(DensityMatrix) - The resulting density matrix after execution.
• Err(QisError) - If the circuit contains unsupported operations.

pub fn apply_circuit(&mut self, circuit: &Circuit) -> Result<(), QisError>

Applies a quantum circuit to this density matrix in-place.

Once bound to circuit qubit IDs, subsequent circuits must have the same IDs in the same order. The mapping is committed only on success; gates completed before an error are not rolled back.

As in from_circuit, terminal measurements are ignored. A measured qubit used by a later gate or Reset is rejected before the state is changed. Independent gates, repeated measurements, Barrier, Delay, and Store are allowed.

The circuit is first decomposed into basic gates via Circuit::decompose, then each operation is applied sequentially to the current density matrix.

Arguments

• circuit - The circuit to apply. Must have the same number of qubits as this density matrix.

Returns

• Ok(()) - If all operations were applied successfully
• Err(QisError) - If the circuit cannot be simulated

Errors

• QisError::InvalidStateDimension if circuit.num_qubits() != self.num_qubits
• QisError::UnsupportedOperation if the circuit contains classical control flow
• QisError::CircuitError if a gate lacks a matrix representation, contains unresolved symbolic parameters, or a referenced qubit is not found

Supported Instructions

• Standard single and multi-qubit gates
• Controlled gates (CX, CY, CZ, CRX, CRY, CRZ)
• Multi-controlled gates (CCX / Toffoli)
• Unitary gates with an explicit matrix representation
• Reset
• Barriers and delays (ignored)
• Measurement output declarations from Circuit::measure* (ignored)

pub fn apply_standard_gate( &mut self, gate: StandardGate, qubits: &[usize], params: &[f64], ) -> Result<(), QisError>

pub fn zeros(num_qubits: usize) -> Self

Creates a new density matrix filled entirely with zeros.

This is not a valid physical state (trace = 0) but is useful as an accumulator during operations like Kraus channel application.

Arguments

• num_qubits - Number of qubits in the system.

Panics

Panics if the dimension is unrepresentable or allocation fails. Use Self::try_zeros to receive an error instead.

pub fn try_zeros(num_qubits: usize) -> Result<Self, QisError>

Fallible variant of Self::zeros.

Errors

Returns QisError::InvalidParameterValue for an unrepresentable dimension, or QisError::UnsupportedOperation if allocation fails.

pub fn trace(&self) -> Complex64

Computes the trace of the density matrix.

For any valid physical state, the trace must equal 1.0.

Returns

The trace (sum of diagonal elements) as a complex number.

pub fn is_hermitian(&self, tol: f64) -> bool

Checks if the density matrix is Hermitian (self-adjoint) within a tolerance.

A valid density matrix must satisfy ρ = ρ†, i.e., ρ_ij = ρ_ji*.

Arguments

• tol - Tolerance for floating-point comparison (e.g., 1e-10).

Returns

true if the matrix is Hermitian within the specified tolerance.

pub fn is_positive_semidefinite_approx(&self, tol: f64) -> bool

Checks if the density matrix is positive semidefinite (approximately).

Uses faer’s self-adjoint eigenvalue decomposition to check that all eigenvalues satisfy λᵢ ≥ -tol.

For Hermitian matrices the eigenvalue check is definitive. For non-Hermitian inputs the decomposition may fail or produce meaningless results, in which case false is returned.

Edge cases

• Matrices containing NaN or Inf return false without panicking.
• Negative eigenvalues smaller in magnitude than tol are treated as numerical noise and pass the check.

Arguments

• tol - Tolerance for floating-point comparison.

Returns

true if all eigenvalues are ≥ -tol, false otherwise.

pub fn validate_physical(&self, tol: f64) -> Result<(), QisError>

Validates all physical constraints of the density matrix.

Checks:

1. Hermiticity: ρ = ρ†
2. Positive semidefiniteness: All eigenvalues >= 0
3. Unit trace: Tr(ρ) = 1

Arguments

• tol - Tolerance for floating-point comparisons (e.g., 1e-10).

Returns

• Ok(()) if all constraints are satisfied.
• Err(QisError::NotHermitian) if not Hermitian.
• Err(QisError::NotPositiveSemidefinite) if not positive semidefinite.
• Err(QisError::NotNormalized) if trace is not 1.

pub fn apply_x(&mut self, qubit: usize) -> Result<(), QisError>

Applies the Pauli-X (NOT) gate to the specified qubit.

Errors

Returns QisError::IndexOutOfBounds if qubit >= num_qubits.

pub fn apply_y(&mut self, qubit: usize) -> Result<(), QisError>

Applies the Pauli-Y gate to the specified qubit.

pub fn apply_z(&mut self, qubit: usize) -> Result<(), QisError>

Applies the Pauli-Z gate to the specified qubit.

pub fn apply_h(&mut self, qubit: usize) -> Result<(), QisError>

Applies the Hadamard gate to the specified qubit.

pub fn apply_u( &mut self, qubit: usize, theta: f64, phi: f64, lambda: f64, ) -> Result<(), QisError>

Applies the generic single-qubit U gate with parameters theta, phi, and lam.

pub fn apply_rx(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the single-qubit rotation about the X-axis by angle theta.

pub fn apply_ry(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the single-qubit rotation about the Y-axis by angle theta.

pub fn apply_rz(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the single-qubit rotation about the Z-axis by angle theta.

pub fn apply_phase(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies a phase shift of theta to the specified qubit.

pub fn apply_s(&mut self, qubit: usize) -> Result<(), QisError>

Applies the S (Phase) gate to the specified qubit.

pub fn apply_sdg(&mut self, qubit: usize) -> Result<(), QisError>

Applies the inverse S (SDG) gate to the specified qubit.

pub fn apply_t(&mut self, qubit: usize) -> Result<(), QisError>

Applies the T (Pi/8) gate to the specified qubit.

pub fn apply_tdg(&mut self, qubit: usize) -> Result<(), QisError>

Applies the inverse T (TDG) gate to the specified qubit.

pub fn apply_x2p(&mut self, qubit: usize) -> Result<(), QisError>

Applies a +Pi/2 rotation about the X-axis.

pub fn apply_x2m(&mut self, qubit: usize) -> Result<(), QisError>

Applies a -Pi/2 rotation about the X-axis.

pub fn apply_y2p(&mut self, qubit: usize) -> Result<(), QisError>

Applies a +Pi/2 rotation about the Y-axis.

pub fn apply_y2m(&mut self, qubit: usize) -> Result<(), QisError>

Applies a -Pi/2 rotation about the Y-axis.

pub fn apply_rxy( &mut self, qubit: usize, theta: f64, phi: f64, ) -> Result<(), QisError>

Applies the parameterized RXY rotation gate.

pub fn apply_xy2p(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the XY2P gate.

pub fn apply_xy2m(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the XY2M gate.

pub fn apply_xy(&mut self, qubit: usize, theta: f64) -> Result<(), QisError>

Applies the XY gate.

Matrix: [[0, -i·e^(-iθ)], [-i·e^(iθ), 0]]

pub fn apply_gphase(&mut self, _phi: f64)

Applies a global phase (has no observable effect on a density matrix).

pub fn apply_cx( &mut self, control: usize, target: usize, ) -> Result<(), QisError>

Applies the Controlled-X (CNOT) gate.

pub fn apply_cz(&mut self, q0: usize, q1: usize) -> Result<(), QisError>

Applies the Controlled-Z gate.

pub fn apply_single_qubit_gate( &mut self, qubit: usize, matrix: [[Complex64; 2]; 2], ) -> Result<(), QisError>

Applies an arbitrary 2x2 unitary matrix to a single qubit.

pub fn apply_two_qubit_gate( &mut self, q0: usize, q1: usize, matrix: [[Complex64; 4]; 4], ) -> Result<(), QisError>

Applies an arbitrary 4x4 unitary matrix to two qubits.

pub fn apply_swap(&mut self, q0: usize, q1: usize) -> Result<(), QisError>

Applies the SWAP gate between two qubits.

pub fn apply_cy( &mut self, control: usize, target: usize, ) -> Result<(), QisError>

Applies the Controlled-Y gate.

pub fn apply_ccx( &mut self, c0: usize, c1: usize, target: usize, ) -> Result<(), QisError>

Applies the Toffoli (Controlled-Controlled-X) gate.

pub fn apply_rxx( &mut self, q0: usize, q1: usize, theta: f64, ) -> Result<(), QisError>

Applies the parameterized RXX (Ising XX) gate.

pub fn apply_ryy( &mut self, q0: usize, q1: usize, theta: f64, ) -> Result<(), QisError>

Applies the parameterized RYY (Ising YY) gate.

pub fn apply_rzz( &mut self, q0: usize, q1: usize, theta: f64, ) -> Result<(), QisError>

Applies the parameterized RZZ (Ising ZZ) gate.

pub fn apply_rzx( &mut self, q0: usize, q1: usize, theta: f64, ) -> Result<(), QisError>

Applies the parameterized RZX gate.

pub fn apply_crx( &mut self, control: usize, target: usize, theta: f64, ) -> Result<(), QisError>

Applies the Controlled-RX gate.

pub fn apply_cry( &mut self, control: usize, target: usize, theta: f64, ) -> Result<(), QisError>

Applies the Controlled-RY gate.

pub fn apply_crz( &mut self, control: usize, target: usize, theta: f64, ) -> Result<(), QisError>

Applies the Controlled-RZ gate.

pub fn apply_fsim( &mut self, q0: usize, q1: usize, theta: f64, phi: f64, ) -> Result<(), QisError>

Applies the Fermionic Simulation (FSIM) gate.

pub fn apply_unitary_gate( &mut self, qubits: &[usize], matrix: &Array2<Complex64>, ) -> Result<(), QisError>

Applies an arbitrary n-qubit unitary gate.

The evolution is given by $\rho \to U \rho U^\dagger$.

Arguments

• qs - Slice of qubit indices the gate acts on.
• mat - The unitary matrix as a $2^n \times 2^n$ ndarray.

pub fn apply_kraus( &mut self, ops: &[Vec<Complex64>], qs: &[usize], ) -> Result<(), QisError>

Applies a general quantum channel specified by Kraus operators.

The evolution of the density matrix is given by $\rho \to \sum_k K_k \rho K_k^\dagger$, where $\sum_k K_k^\dagger K_k = I$ for a trace-preserving channel.

Arguments

• ops - A slice of Kraus operators, where each operator is represented as a flattened vector of Complex64.
• qs - The target qubit indices the channel acts upon.

pub fn partial_trace(&self, keep: &[usize]) -> Result<Self, QisError>

Computes the partial trace over a set of qubits.

Reduces the N-qubit system to a smaller subsystem containing only the specified keep qubits by tracing out all other qubits.

Arguments

• keep - A slice of qubit indices to keep in the resulting reduced density matrix.

Returns

A new DensityMatrix representing the subsystem, with num_qubits = keep.len().

pub fn expectation(&self, h: &dyn Observable) -> Result<f64, QisError>

Computes the expectation value of a Hamiltonian observable.

Calculates Tr(ρ * H) for the current density matrix ρ and a given Hamiltonian H. The Hamiltonian is represented as a sum of Pauli strings with coefficients.

Arguments

• h - The Hamiltonian observable.

Returns

The expectation value as a real number (f64), or a QisError if the qubit counts do not match.

Errors

Returns QisError::CircuitError if the observable acts on a different number of qubits than the density matrix.

pub fn measure(&mut self, qubit: usize) -> Result<bool, QisError>

Measures qubit in the Z basis, collapsing the density matrix in-place.

Returns true for outcome |1⟩ and false for outcome |0⟩. The density matrix is projected into the subspace consistent with the outcome and renormalized: ρ' = Π_b ρ Π_b / Tr(Π_b ρ).

pub fn measure_with_rng( &mut self, qubit: usize, rng: &mut impl Rng, ) -> Result<bool, QisError>

Measures one qubit using the supplied random stream.

pub fn reset(&mut self, qubit: usize) -> Result<(), QisError>

Resets qubit to the |0⟩ state.

Applies the deterministic reset channel with Kraus operators |0⟩⟨0| and |0⟩⟨1|. The outcome is discarded, preserving the reduced density matrix of the other qubits, which may be mixed for an entangled input.

pub fn measure_all(&mut self) -> Outcome

Measures all qubits sequentially, returning a bit-packed Outcome.

Equivalent to calling measure(q) for each qubit q in order 0..num_qubits. The density matrix is fully collapsed after this call. Use Outcome::is_one(q) to read qubit q’s result.

pub fn measure_all_with_rng(&mut self, rng: &mut impl Rng) -> Outcome

Measures all qubits in storage order using the supplied random stream.

pub fn sample_shots(&self, shots: usize) -> Vec<Outcome>

Samples shots independent measurement outcomes in parallel.

Each Rayon worker thread reuses a single pre-allocated clone of the density matrix (via reset_from), avoiding per-shot heap allocation. Returns a Vec<Outcome> of bit-packed results.

Example

use cqlib_core::qis::DensityMatrix;

let mut dm = DensityMatrix::new(2);
dm.apply_h(0).unwrap();
dm.apply_cx(0, 1).unwrap(); // |Φ⁺⟩ Bell state

let shots = dm.sample_shots(200);
// All outcomes must be |00⟩ or |11⟩
assert!(shots.iter().all(|v| v.is_one(0) == v.is_one(1)));

pub fn sample_shots_with_seed( &self, shots: usize, seed: Option<u64>, ) -> Vec<Outcome>

Samples without changing this state. An explicit seed assigns a random stream to each shot before parallel execution, independent of scheduling.

pub fn sample( &self, measurement: &Measurement, shots: usize, ) -> Result<ExecutionResult, QisError>

Samples the state using a circuit Measurement as the output contract.

measurement.qubits()[i] becomes bit i in each Outcome. The returned ExecutionResult reuses the device result shape, with counts projected to only the qubits named by the measurement.

pub fn sample_with_seed( &self, measurement: &Measurement, shots: usize, seed: Option<u64>, ) -> Result<ExecutionResult, QisError>

Projects seeded full-state samples onto the requested measurement.

pub fn probs( &self, measurement: &Measurement, ) -> Result<HashMap<Outcome, f64>, QisError>

Returns the marginal probability distribution selected by a circuit Measurement.

This method does not execute circuit IR. It treats the Measurement as a self-contained output contract and aggregates computational-basis probabilities over unmeasured qubits.

Trait Implementations

impl AddAssign for DensityMatrix

fn add_assign(&mut self, rhs: Self)

Performs the += operation.

impl Clone for DensityMatrix

fn clone(&self) -> Self

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for DensityMatrix

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl StateVisualizationSource for DensityMatrix

fn num_qubits(&self) -> usize

Number of qubits represented by the state.

fn density_matrix_data(&self) -> Result<Vec<Complex64>, VisualizationError>

Return row-major density-matrix data for visualization.