Random Objects and Testing

def rand_ket(N: int, density: float = 1, dims: list = None) -> Qobj:
    """
    Generate random state vector.
    
    Parameters:
    - N: Dimension of Hilbert space
    - density: Sparsity parameter (1 = dense, <1 = sparse)
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random normalized state vector |ψ⟩
    """

def rand_dm(N: int, density: float = 0.75, dims: list = None) -> Qobj:
    """
    Generate random density matrix.
    
    Parameters:
    - N: Dimension of Hilbert space
    - density: Sparsity parameter
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random positive semidefinite density matrix
    """

def rand_herm(N: int, density: float = 0.75, dims: list = None) -> Qobj:
    """
    Generate random Hermitian operator.
    
    Parameters:
    - N: Matrix dimension
    - density: Sparsity parameter
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random Hermitian operator H = H†
    """

def rand_unitary(N: int, density: float = 0.75, dims: list = None) -> Qobj:
    """
    Generate random unitary operator.
    
    Parameters:
    - N: Matrix dimension
    - density: Sparsity parameter
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random unitary operator U†U = UU† = I
    """

def rand_stochastic(N: int, density: float = 0.75, kind: str = 'left') -> Qobj:
    """
    Generate random stochastic matrix.
    
    Parameters:
    - N: Matrix dimension
    - density: Sparsity parameter
    - kind: 'left', 'right', or 'doubly' stochastic
    
    Returns:
    - Qobj: Random stochastic matrix
    """

def rand_kraus_map(N: int, dims: list = None) -> list:
    """
    Generate random quantum channel as Kraus operators.
    
    Parameters:
    - N: Dimension of Hilbert space
    - dims: Dimensions for tensor structure
    
    Returns:
    - list: List of Kraus operators {Kᵢ} with ∑ᵢ Kᵢ†Kᵢ = I
    """

def rand_super(N: int, dims: list = None) -> Qobj:
    """
    Generate random superoperator.
    
    Parameters:
    - N: Dimension of underlying Hilbert space
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random completely positive superoperator
    """

def rand_super_bcsz(N: int, enforce_tp: bool = True, rank: int = None, dims: list = None) -> Qobj:
    """
    Generate random superoperator using BCSZ construction.
    
    Parameters:
    - N: Dimension of Hilbert space
    - enforce_tp: Enforce trace preservation
    - rank: Rank of Choi matrix (default: N²)
    - dims: Dimensions for tensor structure
    
    Returns:
    - Qobj: Random CPTP or CP superoperator
    """

import qutip as qt
import numpy as np

# Random state vectors
N = 10
psi_random = qt.rand_ket(N)
print(f"Random state norm: {psi_random.norm():.6f}")  # Should be 1
print(f"Random state shape: {psi_random.shape}")

# Sparse random state
psi_sparse = qt.rand_ket(N, density=0.3)  # 30% non-zero elements
print(f"Sparse state sparsity: {psi_sparse.data.nnz / (N * 1):.2f}")

# Random density matrices
rho_random = qt.rand_dm(N)
print(f"Random DM trace: {rho_random.tr():.6f}")  # Should be 1
print(f"Random DM positivity: {all(rho_random.eigenenergies() >= -1e-12)}")

# Check density matrix properties
eigenvals = rho_random.eigenenergies()
print(f"DM eigenvalues range: [{eigenvals.min():.6f}, {eigenvals.max():.6f}]")
print(f"DM purity: {(rho_random**2).tr():.3f}")

# Random Hermitian operators
H_random = qt.rand_herm(N)
print(f"Hermitian check: {(H_random - H_random.dag()).norm():.2e}")  # Should be ~0

# Random unitary operators
U_random = qt.rand_unitary(N)
unitarity_check = (U_random * U_random.dag() - qt.qeye(N)).norm()
print(f"Unitarity check: {unitarity_check:.2e}")  # Should be ~0

# Verify unitary preserves norms
psi_test = qt.rand_ket(N)
psi_evolved = U_random * psi_test
print(f"Norm preservation: {psi_test.norm():.6f} -> {psi_evolved.norm():.6f}")

# Random stochastic matrices
S_left = qt.rand_stochastic(N, kind='left')  # Columns sum to 1
S_right = qt.rand_stochastic(N, kind='right')  # Rows sum to 1
S_doubly = qt.rand_stochastic(N, kind='doubly')  # Both rows and columns sum to 1

# Check stochasticity
col_sums = np.sum(S_left.full(), axis=0)
row_sums = np.sum(S_right.full(), axis=1)
print(f"Left stochastic column sums: {col_sums[:3]}...")  # Should be all 1s
print(f"Right stochastic row sums: {row_sums[:3]}...")    # Should be all 1s

# Random quantum channels
kraus_ops = qt.rand_kraus_map(4)  # 4-dimensional system
print(f"Number of Kraus operators: {len(kraus_ops)}")

# Verify completeness relation: ∑ᵢ Kᵢ†Kᵢ = I
completeness = sum(K.dag() * K for K in kraus_ops)
identity_check = (completeness - qt.qeye(4)).norm()
print(f"Kraus completeness check: {identity_check:.2e}")

# Apply random channel to state
rho_in = qt.rand_dm(4)
rho_out = sum(K * rho_in * K.dag() for K in kraus_ops)
print(f"Channel preserves trace: {rho_in.tr():.6f} -> {rho_out.tr():.6f}")

# Random superoperator
S = qt.rand_super(3)  # 3x3 system -> 9x9 superoperator
print(f"Superoperator shape: {S.shape}")

# Test on random density matrix
rho_test = qt.rand_dm(3)
rho_vec = qt.operator_to_vector(rho_test)
rho_out_vec = S * rho_vec
rho_out = qt.vector_to_operator(rho_out_vec)
print(f"Superoperator trace preservation: {rho_test.tr():.6f} -> {rho_out.tr():.6f}")

# BCSZ random superoperator (guaranteed CPTP)
S_bcsz = qt.rand_super_bcsz(3, enforce_tp=True)
rho_bcsz_vec = S_bcsz * rho_vec
rho_bcsz = qt.vector_to_operator(rho_bcsz_vec)
print(f"BCSZ trace preservation: {rho_bcsz.tr():.6f}")

# Benchmarking with random objects
def benchmark_expectation_values(N, num_trials=100):
    """Benchmark expectation value calculations."""
    import time
    
    # Generate random objects
    operators = [qt.rand_herm(N) for _ in range(10)]
    states = [qt.rand_ket(N) for _ in range(num_trials)]
    
    start_time = time.time()
    for psi in states:
        for op in operators:
            _ = qt.expect(op, psi)
    end_time = time.time()
    
    total_calculations = num_trials * len(operators)
    avg_time = (end_time - start_time) / total_calculations
    print(f"Average expectation value time (N={N}): {avg_time*1000:.3f} ms")

# Run benchmarks
for N in [10, 50, 100]:
    benchmark_expectation_values(N, num_trials=20)

# Random object properties analysis
def analyze_random_properties(N=20, num_samples=100):
    """Analyze statistical properties of random objects."""
    
    # Random Hermitian eigenvalue statistics
    eigenvals_all = []
    for _ in range(num_samples):
        H = qt.rand_herm(N)
        eigenvals_all.extend(H.eigenenergies())
    
    eigenvals_all = np.array(eigenvals_all)
    print(f"Random Hermitian eigenvalues:")
    print(f"  Mean: {eigenvals_all.mean():.3f}")
    print(f"  Std:  {eigenvals_all.std():.3f}")
    print(f"  Range: [{eigenvals_all.min():.3f}, {eigenvals_all.max():.3f}]")
    
    # Random density matrix purity statistics
    purities = []
    for _ in range(num_samples):
        rho = qt.rand_dm(N)
        purities.append((rho**2).tr())
    
    purities = np.array(purities)
    print(f"Random density matrix purities:")
    print(f"  Mean: {purities.mean():.3f}")
    print(f"  Std:  {purities.std():.3f}")
    print(f"  Range: [{purities.min():.3f}, {purities.max():.3f}]")
    print(f"  (Maximum purity = 1.0, minimum = 1/{N} = {1/N:.3f})")

analyze_random_properties()

# Tensor product random objects
# Multi-qubit random states
dims = [2, 2, 2]  # 3-qubit system
psi_3qubit = qt.rand_ket(8, dims=dims)
print(f"3-qubit random state dims: {psi_3qubit.dims}")

# Random two-qubit gate
U_2qubit = qt.rand_unitary(4, dims=[[2,2], [2,2]])
print(f"2-qubit gate dims: {U_2qubit.dims}")

# Test gate on product states
psi_00 = qt.tensor(qt.basis(2,0), qt.basis(2,0))
psi_entangled = U_2qubit * psi_00

# Measure entanglement
rho_total = psi_entangled * psi_entangled.dag()
rho_A = qt.ptrace(rho_total, 0)
entanglement = qt.entropy_vn(rho_A)
print(f"Entanglement generated by random gate: {entanglement:.3f}")

# Random quantum circuit simulation
def random_circuit(n_qubits, depth, gate_prob=0.7):
    """Generate random quantum circuit."""
    N = 2**n_qubits
    U_total = qt.qeye([2]*n_qubits)
    
    for layer in range(depth):
        # Random single-qubit gates
        for qubit in range(n_qubits):
            if np.random.random() < gate_prob:
                # Random single-qubit unitary
                U_local = qt.rand_unitary(2)
                
                # Embed in full space
                gate_list = [qt.qeye(2) for _ in range(n_qubits)]
                gate_list[qubit] = U_local
                U_layer = qt.tensor(*gate_list)
                U_total = U_layer * U_total
        
        # Random two-qubit gates
        for qubit in range(n_qubits-1):
            if np.random.random() < gate_prob/2:  # Less frequent
                # Random CNOT-like gate
                if np.random.random() < 0.5:
                    gate = qt.cnot(N=n_qubits, control=qubit, target=qubit+1)
                else:
                    gate = qt.cz_gate(N=n_qubits, control=qubit, target=qubit+1)
                U_total = gate * U_total
    
    return U_total

# Test random circuit
n_qubits = 3
circuit_depth = 5
U_circuit = random_circuit(n_qubits, circuit_depth)

# Apply to initial state
initial_state = qt.tensor(*[qt.basis(2,0) for _ in range(n_qubits)])  # |000⟩
final_state = U_circuit * initial_state

print(f"Random circuit unitarity: {(U_circuit * U_circuit.dag() - qt.qeye(2**n_qubits)).norm():.2e}")
print(f"Final state norm: {final_state.norm():.6f}")

# Measure final state in computational basis
probs = [(abs(final_state[i,0])**2) for i in range(2**n_qubits)]
print(f"Computational basis probabilities: {[f'{p:.3f}' for p in probs]}")

# All random object functions return Qobj instances with appropriate
# type ('ket', 'oper', 'super') and specified dimensions