diff --git a/Cargo.lock b/Cargo.lock index a23d7d0f..438edc76 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -165,7 +165,7 @@ checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" [[package]] name = "binding-c" -version = "0.1.0-beta.5" +version = "0.1.0-beta.6" dependencies = [ "bitvec", "cbindgen", @@ -178,7 +178,7 @@ dependencies = [ [[package]] name = "binding-python" -version = "2.0.0-beta5" +version = "2.0.0-beta6" dependencies = [ "chrono", "cqlib-core", @@ -465,14 +465,14 @@ dependencies = [ [[package]] name = "cqlib" -version = "0.1.0-beta.5" +version = "0.1.0-beta.6" dependencies = [ "cqlib-core", ] [[package]] name = "cqlib-core" -version = "0.1.0-beta.5" +version = "0.1.0-beta.6" dependencies = [ "approx", "bitvec", diff --git a/Cargo.toml b/Cargo.toml index 12d91b1d..32a64e5c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,7 +8,7 @@ members = [ resolver = "2" [workspace.package] -version = "0.1.0-beta.5" +version = "0.1.0-beta.6" edition = "2024" rust-version = "1.89" license = "Apache-2.0" diff --git a/crates/binding-c/Cargo.toml b/crates/binding-c/Cargo.toml index 216f3683..df8ae605 100644 --- a/crates/binding-c/Cargo.toml +++ b/crates/binding-c/Cargo.toml @@ -9,7 +9,7 @@ license.workspace = true crate-type = ["rlib", "staticlib", "cdylib"] [dependencies] -cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.5" } +cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.6" } libc = "0.2" bitvec.workspace = true ndarray.workspace = true diff --git a/crates/binding-c/src/circuit/parameter.rs b/crates/binding-c/src/circuit/parameter.rs index 2f9dfa0e..8ae00dfd 100644 --- a/crates/binding-c/src/circuit/parameter.rs +++ b/crates/binding-c/src/circuit/parameter.rs @@ -445,7 +445,8 @@ pub extern "C" fn param_replace( box_param(unsafe { (*a).inner.replace(symbol, replacement) }) } -/// Substitute several symbols in one call and simplify the result. +/// Simultaneously substitute symbols and simplify the result. Replacement +/// expressions are not substituted again; the original parameter is unchanged. /// `names` and `values` are parallel arrays of length `len`. Returns a new /// parameter, or NULL on NULL input, non-UTF-8 names or NULL elements. #[unsafe(no_mangle)] diff --git a/crates/binding-c/src/qis/density_matrix_noise.rs b/crates/binding-c/src/qis/density_matrix_noise.rs index c84d7c08..85d7be0b 100644 --- a/crates/binding-c/src/qis/density_matrix_noise.rs +++ b/crates/binding-c/src/qis/density_matrix_noise.rs @@ -682,6 +682,10 @@ pub extern "C" fn density_matrix_noise_probabilities_with_readout_len( /// Copy the probability distribution including readout errors configured in /// the noise model for `qubits` into `buffer`. Returns 0 on success. /// +/// `qubits` contains storage positions. For a circuit-bound state, noise is +/// looked up using the circuit ID at each position; unbound states use identity +/// IDs. The full output distribution retains storage order. +/// /// `qubits` may be NULL when `qubit_len` is 0, in which case no readout /// error is applied and the result matches `density_matrix_noise_probabilities`. #[unsafe(no_mangle)] diff --git a/crates/binding-c/tests/parameter_test.rs b/crates/binding-c/tests/parameter_test.rs index 7afb1e84..87dc059d 100644 --- a/crates/binding-c/tests/parameter_test.rs +++ b/crates/binding-c/tests/parameter_test.rs @@ -22,6 +22,63 @@ fn eval(ptr: *const CParameter, bindings: &str) -> f64 { param_evaluate(ptr, cstr(bindings).as_ptr()) } +#[test] +fn param_substitution_swaps_symbols_simultaneously() { + let x = param_symbol(cstr("x").as_ptr()); + let y = param_symbol(cstr("y").as_ptr()); + let expression = param_parse(cstr("x + 2*y").as_ptr()); + let name_x = cstr("x"); + let name_y = cstr("y"); + for (names, values) in [ + ( + [name_x.as_ptr(), name_y.as_ptr()], + [y.cast_const(), x.cast_const()], + ), + ( + [name_y.as_ptr(), name_x.as_ptr()], + [x.cast_const(), y.cast_const()], + ), + ] { + for _ in 0..100 { + let result = param_substitute_many(expression, names.as_ptr(), values.as_ptr(), 2); + assert!(!result.is_null()); + let value = eval(result, "x:1,y:3"); + param_free(result); + assert_eq!(value, 5.0); + } + } + assert_eq!(eval(expression, "x:1,y:3"), 7.0); + for ptr in [x, y, expression] { + param_free(ptr); + } +} + +#[test] +fn param_substitution_reaches_polynomial_symbols() { + let x = param_symbol(cstr("x").as_ptr()); + let y = param_symbol(cstr("y").as_ptr()); + let square = param_mul(x, x); + let sum = param_add(square, x); + let one = param_from_double(1.0); + let expression = param_add(sum, one); + let name = cstr("x"); + let replaced = param_replace(expression, name.as_ptr(), y); + let substituted = param_substitute_many( + expression, + [name.as_ptr()].as_ptr(), + [y.cast_const()].as_ptr(), + 1, + ); + for result in [replaced, substituted] { + assert!(!result.is_null()); + assert_eq!(eval(result, "y:2"), 7.0); + param_free(result); + } + for ptr in [x, y, square, sum, one, expression] { + param_free(ptr); + } +} + #[test] fn param_construction_and_constants() { let theta = param_symbol(cstr("theta").as_ptr()); diff --git a/crates/binding-c/tests/qis_noise_test.rs b/crates/binding-c/tests/qis_noise_test.rs index 989e236e..21f69a3b 100644 --- a/crates/binding-c/tests/qis_noise_test.rs +++ b/crates/binding-c/tests/qis_noise_test.rs @@ -2,7 +2,9 @@ #![allow(dead_code)] -use binding_c::circuit::{circuit_cx, circuit_free, circuit_h, circuit_new, circuit_x}; +use binding_c::circuit::{ + circuit_cx, circuit_free, circuit_from_qubits, circuit_h, circuit_new, circuit_x, +}; use binding_c::cqlib_string_free; use binding_c::device::{ NOISE_BIT_FLIP, noise_model_add_readout, noise_model_add_single_qubit, noise_model_free, @@ -216,6 +218,42 @@ fn test_readout_error_difference() { noise_model_free(nm); } +#[test] +fn test_readout_uses_circuit_ids_at_storage_positions() { + for ids in [[0, 1], [1, 0], [2, 0], [9, 2]] { + let circuit = circuit_from_qubits(ids.as_ptr(), ids.len()); + assert!(!circuit.is_null()); + let model = noise_model_new(); + assert_eq!(noise_model_add_readout(model, ids[0], 0.0, 1.0), 0); + let sim = density_matrix_noise_from_circuit(circuit, model); + assert!(!sim.is_null()); + assert_eq!(density_matrix_noise_probabilities_with_readout_len(sim), 4); + for (positions, expected) in [ + (&[0][..], [0.0, 1.0, 0.0, 0.0]), + (&[1][..], [1.0, 0.0, 0.0, 0.0]), + (&[1, 0][..], [0.0, 1.0, 0.0, 0.0]), + (&[][..], [1.0, 0.0, 0.0, 0.0]), + ] { + let mut actual = [0.0; 4]; + assert_eq!( + density_matrix_noise_probabilities_with_readout( + sim, + actual.as_mut_ptr(), + actual.len(), + positions.as_ptr(), + positions.len() + ), + 0 + ); + assert_eq!(actual, expected); + } + assert_eq!(read_probs(sim, 4), vec![1.0, 0.0, 0.0, 0.0]); + density_matrix_noise_free(sim); + noise_model_free(model); + circuit_free(circuit); + } +} + #[test] fn test_multi_qubit_readout() { // Mirrors core `test_multi_qubit_readout_error`. diff --git a/crates/binding-c/tests/qis_test.rs b/crates/binding-c/tests/qis_test.rs index 05e2174d..ddc13132 100644 --- a/crates/binding-c/tests/qis_test.rs +++ b/crates/binding-c/tests/qis_test.rs @@ -165,6 +165,61 @@ fn test_density_matrix_partial_trace() { density_matrix_free(dm); } +#[test] +fn density_matrix_duplicate_gate_targets_preserve_state() { + let dm = density_matrix_new(2); + assert!(!dm.is_null()); + assert_eq!(density_matrix_apply_h(dm, 0), 0); + let mut before = vec![num_complex::Complex64::default(); 16]; + assert_eq!(density_matrix_data(dm, before.as_mut_ptr(), 16), 0); + for q in [0, 1] { + assert_eq!(density_matrix_apply_cx(dm, q, q), -8); + assert_eq!(density_matrix_apply_cz(dm, q, q), -8); + let mut after = before.clone(); + assert_eq!(density_matrix_data(dm, after.as_mut_ptr(), 16), 0); + assert_eq!(after, before); + } + assert_eq!(density_matrix_apply_cx(dm, 0, 1), 0); + density_matrix_free(dm); +} + +#[test] +fn density_matrix_duplicate_kraus_targets_preserve_state() { + let state: Vec<_> = (1..=4) + .map(|i| num_complex::Complex64::new(i as f64 / 30.0_f64.sqrt(), 0.0)) + .collect(); + let dm = density_matrix_from_state(2, state.as_ptr(), state.len()); + assert!(!dm.is_null()); + let vector = [4.0, -1.0, 2.0, 3.0]; + let operator: Vec<_> = (0..4) + .flat_map(|r| { + (0..4).map(move |c| { + num_complex::Complex64::new( + if r == c { 1.0 } else { 0.0 } - 2.0 * vector[r] * vector[c] / 30.0, + 0.0, + ) + }) + }) + .collect(); + let mut before = vec![num_complex::Complex64::default(); 16]; + assert_eq!(density_matrix_data(dm, before.as_mut_ptr(), 16), 0); + for q in [0, 1] { + let targets = [q, q]; + assert_eq!( + density_matrix_apply_kraus(dm, operator.as_ptr(), 1, 4, targets.as_ptr(), 2), + -8 + ); + let mut after = before.clone(); + assert_eq!(density_matrix_data(dm, after.as_mut_ptr(), 16), 0); + assert_eq!(after, before); + } + assert_eq!( + density_matrix_apply_kraus(dm, operator.as_ptr(), 1, 4, [0, 1].as_ptr(), 2), + 0 + ); + density_matrix_free(dm); +} + #[test] fn test_stabilizer_new_and_clifford() { let ss = stabilizer_new(2); diff --git a/crates/binding-c/tests/test_circuit.c b/crates/binding-c/tests/test_circuit.c index f1f33944..ba48098f 100644 --- a/crates/binding-c/tests/test_circuit.c +++ b/crates/binding-c/tests/test_circuit.c @@ -165,11 +165,98 @@ static void test_extended_gates(void) { circuit_free(circuit); } +static void test_circuit_gate_operation_matrix(void) { + const uint32_t orders[][2] = {{0, 1}, {1, 0}, {7, 2}}; + const uint32_t targets[] = {0, 1}; + for (size_t order = 0; order < 3; ++order) { + for (int parameterized = 0; parameterized < 2; ++parameterized) { + CCircuit* body = circuit_new(0); + assert(body != NULL); + assert(circuit_add_qubits(body, orders[order], 2) == 0); + /* Build on dense ids, then map into the explicit storage order. */ + CCircuit* definition = circuit_new(2); + assert(definition != NULL); + CParameter* theta = param_parse("theta"); + assert(theta != NULL); + if (parameterized) { + assert(circuit_rx_param(definition, 0, theta) == 0); + } else { + assert(circuit_x(definition, 0) == 0); + } + assert(circuit_compose(body, definition, orders[order], 2) == 0); + circuit_free(definition); + assert(circuit_set_global_phase(body, 0.13) == 0); + CCircuitGate* gate = circuit_to_gate(body, "first_gate"); + assert(gate != NULL); + for (size_t application = 0; application < 2; ++application) { + double angle = application == 0 ? 0.37 : -0.61; + CParameter* argument = param_from_double(angle); + const CParameter* arguments[] = {argument}; + assert(argument != NULL); + CCircuit* circuit = circuit_new(2); + assert(circuit != NULL); + assert(circuit_circuit_gate(circuit, gate, targets, 2, + parameterized ? arguments : NULL, + parameterized ? 1 : 0) == 0); + /* CParameter arguments are stored as parameter-table indices. + * Resolve them in the owning circuit before extracting the raw + * DAG Operation, which has no access to that parameter table. + * The gate definition still contains theta; matrix() must bind it. */ + CCircuit* resolved = circuit_assign_params(circuit, ""); + assert(resolved != NULL && circuit_num_parameters(resolved) == 0); + circuit_free(circuit); + circuit = resolved; + CCircuitDag* dag = circuit_dag_from_circuit(circuit); + assert(dag != NULL && circuit_dag_op_nodes_len(dag) == 1); + uint32_t node; + assert(circuit_dag_op_nodes(dag, &node, 1) == 0); + COperation* operation = circuit_dag_operation(dag, node); + assert(operation != NULL); + size_t rows = 0, cols = 0; + assert(operation_matrix_len(operation) == 16); + assert(operation_matrix_dims(operation, &rows, &cols) == 0); + assert(rows == 4 && cols == 4); + double local[32], full[32]; + assert(operation_matrix(operation, local, 32) == 16); + assert(circuit_to_matrix(circuit, NULL, 0, full, 32) == 16); + for (size_t row = 0; row < 4; ++row) { + for (size_t col = 0; col < 4; ++col) { + for (int local_basis = 0; local_basis < 2; ++local_basis) { + size_t flip = local_basis ? 2 : 1; + double re = 0, im = 0; + if (parameterized) { + if (row == col) re = cos(angle / 2); + if (row == (col ^ flip)) im = -sin(angle / 2); + } else if (row == (col ^ flip)) { + re = 1; + } + double expected_re = re * cos(0.13) - im * sin(0.13); + double expected_im = re * sin(0.13) + im * cos(0.13); + const double* matrix = local_basis ? local : full; + size_t index = 2 * (row * 4 + col); + assert(fabs(matrix[index] - expected_re) < 1e-12); + assert(fabs(matrix[index + 1] - expected_im) < 1e-12); + } + } + } + operation_free(operation); + circuit_dag_free(dag); + circuit_free(circuit); + param_free(argument); + } + circuit_gate_free(gate); + param_free(theta); + circuit_free(body); + } + } +} + int main(void) { test_basic_circuit(); test_errors(); test_symbolic_parameters(); test_extended_gates(); + test_circuit_gate_operation_matrix(); printf("binding-c circuit tests passed\n"); return 0; } diff --git a/crates/binding-c/tests/test_density_matrix_targets.c b/crates/binding-c/tests/test_density_matrix_targets.c new file mode 100644 index 00000000..342c8401 --- /dev/null +++ b/crates/binding-c/tests/test_density_matrix_targets.c @@ -0,0 +1,62 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 2026 +// +// Licensed under the Apache License, Version 2.0. See LICENSE.txt or +// http://www.apache.org/licenses/LICENSE-2.0. + +#include +#include +#include +#include + +#include "cqlib_c.h" + +int main(void) { + Complex64 amplitudes[4]; + Complex64 operator[16]; + const double vector[4] = {4.0, -1.0, 2.0, 3.0}; + for (int i = 0; i < 4; ++i) { + amplitudes[i] = (Complex64){(i + 1) / sqrt(30.0), 0.0}; + for (int j = 0; j < 4; ++j) { + operator[4 * i + j] = + (Complex64){(i == j ? 1.0 : 0.0) - 2.0 * vector[i] * vector[j] / 30.0, 0.0}; + } + } + CDensityMatrix* dm = density_matrix_from_state(2, amplitudes, 4); + assert(dm != NULL); + Complex64 before[16], after[16]; + assert(density_matrix_data(dm, before, 16) == 0); + for (uint32_t q = 0; q < 2; ++q) { + const uint32_t targets[2] = {q, q}; + /* Check Kraus first: a panic must fail the host-process test. */ + assert(density_matrix_apply_kraus(dm, operator, 1, 4, targets, 2) == -8); + assert(density_matrix_data(dm, after, 16) == 0); + assert(memcmp(before, after, sizeof(before)) == 0); + assert(density_matrix_apply_cx(dm, q, q) == -8); + assert(density_matrix_apply_cz(dm, q, q) == -8); + assert(density_matrix_data(dm, after, 16) == 0); + assert(memcmp(before, after, sizeof(before)) == 0); + } + /* A valid channel must still work after the rejected operations. */ + const uint32_t targets[2] = {0, 1}; + assert(density_matrix_apply_kraus(dm, operator, 1, 4, targets, 2) == 0); + assert(density_matrix_data(dm, after, 16) == 0); + double evolved[4] = {0}; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + evolved[i] += operator[4 * i + j].re * amplitudes[((j & 1) << 1) | (j >> 1)].re; + } + } + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + const int row = ((i & 1) << 1) | (i >> 1); + const int col = ((j & 1) << 1) | (j >> 1); + assert(fabs(after[4 * row + col].re - evolved[i] * evolved[j]) < 1e-12); + assert(fabs(after[4 * row + col].im) < 1e-12); + } + } + density_matrix_free(dm); + puts("density matrix target validation passed"); + return 0; +} diff --git a/crates/binding-python/Cargo.toml b/crates/binding-python/Cargo.toml index 942962ac..926ce4ec 100644 --- a/crates/binding-python/Cargo.toml +++ b/crates/binding-python/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "binding-python" -version = "2.0.0-beta5" +version = "2.0.0-beta6" edition.workspace = true rust-version.workspace = true license-file.workspace = true @@ -14,7 +14,7 @@ crate-type = ["cdylib", "rlib"] [dependencies] pyo3 = { version = "0.29", features = ["abi3-py311", "num-complex", "chrono"] } chrono = "0.4" -cqlib-core = { path = "../cqlib-core", version = '0.1.0-beta.5' } +cqlib-core = { path = "../cqlib-core", version = '0.1.0-beta.6' } num-complex.workspace = true numpy.workspace = true rustworkx-core.workspace = true diff --git a/crates/binding-python/cqlib/circuit/parameter.pyi b/crates/binding-python/cqlib/circuit/parameter.pyi index 3e67a27d..fd7a8be4 100644 --- a/crates/binding-python/cqlib/circuit/parameter.pyi +++ b/crates/binding-python/cqlib/circuit/parameter.pyi @@ -123,7 +123,10 @@ class Parameter: ... def substitute(self, bindings: dict[str, Parameter]) -> Parameter: - """Substitute multiple symbols and simplify the resulting expression.""" + """Substitute simultaneously and simplify, without revisiting replacements. + + The original parameter is unchanged. + """ ... def provably_equal(self, other: Parameter, tolerance: float = 1e-12) -> bool: diff --git a/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi b/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi index f0fadb76..67706b5f 100644 --- a/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi +++ b/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi @@ -121,10 +121,17 @@ class DensityMatrixNoise: """Computes measurement probabilities with readout error modeling. Args: - qubits: Indices of qubits to measure + qubits: Storage positions on which to apply readout errors. For a + circuit-bound state, noise is looked up using the circuit ID at + each position; unbound states use identity IDs. Returns: - A vector of probabilities for all 2^n computational basis states. + A vector of probabilities for all 2^num_qubits computational basis + states in storage order, independent of the selection order or size. + An empty selection applies no readout errors. + + Raises: + IndexError: A position is out of range, even without a noise model. """ ... diff --git a/crates/binding-python/cqlib/visualization/__init__.py b/crates/binding-python/cqlib/visualization/__init__.py index aef02055..f7c2bdd0 100644 --- a/crates/binding-python/cqlib/visualization/__init__.py +++ b/crates/binding-python/cqlib/visualization/__init__.py @@ -196,7 +196,12 @@ def plot_state_paulivec( figsize=None, output_path=None, ): - """Render Pauli-basis expectation values as SVG.""" + """Render Pauli-basis expectation values as SVG. + + Labels follow PauliString by default: the rightmost character acts on + qubit 0, so IX measures X on qubit 0 of a two-qubit state. Setting + reverse_bits=True reverses each label without changing its coefficient. + """ svg = _visualization_module.plot_state_paulivec( state, title=title, diff --git a/crates/binding-python/src/circuit/parameter.rs b/crates/binding-python/src/circuit/parameter.rs index dae61327..7f11c551 100644 --- a/crates/binding-python/src/circuit/parameter.rs +++ b/crates/binding-python/src/circuit/parameter.rs @@ -228,7 +228,8 @@ impl PyParameter { self.inner.as_symbol() } - /// Substitutes multiple symbols and simplifies the resulting expression. + /// Simultaneously substitutes symbols and simplifies the resulting expression. + /// Inserted expressions are not substituted again; the original is unchanged. fn substitute(&self, bindings: HashMap) -> Self { let bindings = bindings .into_iter() diff --git a/crates/binding-python/src/qis/hamiltonian.rs b/crates/binding-python/src/qis/hamiltonian.rs index afae3ee4..39d21b75 100644 --- a/crates/binding-python/src/qis/hamiltonian.rs +++ b/crates/binding-python/src/qis/hamiltonian.rs @@ -76,16 +76,6 @@ impl fmt::Display for PyHamiltonian { /// - complex -> complex number /// - tuple (real, imag) -> complex number fn extract_complex(value: &Bound<'_, PyAny>) -> PyResult { - // Try to extract as float first - if let Ok(val) = value.extract::() { - return Ok(Complex64::new(val, 0.0)); - } - // Try to extract as Python complex - if let Ok(py_complex) = value.cast::() { - return Ok(Complex64::new(py_complex.real(), py_complex.imag())); - } - - // Try to extract as tuple (real, imag) if let Ok(tuple) = value.cast::() && tuple.len() == 2 { @@ -94,6 +84,13 @@ fn extract_complex(value: &Bound<'_, PyAny>) -> PyResult { return Ok(Complex64::new(real, imag)); } + // Use the complex protocol before real conversion: NumPy complex scalars + // permit float extraction, but discard the imaginary part with a warning. + // Propagate conversion errors instead of falling back to a lossy float cast. + if let Some(coefficient) = super::fermion::scalar_operand(value)? { + return Ok(coefficient); + } + Err(PyValueError::new_err( "Coefficient must be a float, int, complex, or tuple (real, imag)", )) diff --git a/crates/binding-python/src/qis/state/density_matrix_noise.rs b/crates/binding-python/src/qis/state/density_matrix_noise.rs index 8916923b..a70358d3 100644 --- a/crates/binding-python/src/qis/state/density_matrix_noise.rs +++ b/crates/binding-python/src/qis/state/density_matrix_noise.rs @@ -170,10 +170,17 @@ impl PyDensityMatrixNoise { /// Computes measurement probabilities with readout error modeling. /// /// Args: - /// qubits: Indices of qubits to measure + /// qubits: Storage positions on which to apply readout errors. For a + /// circuit-bound state, noise is looked up using the circuit ID at + /// each position; unbound states use identity IDs. /// /// Returns: - /// A vector of probabilities for all 2^n computational basis states. + /// A vector of probabilities for all 2^num_qubits computational basis + /// states in storage order, independent of the selection order or size. + /// An empty selection applies no readout errors. + /// + /// Raises: + /// IndexError: A position is out of range, even without a noise model. fn probabilities_with_readout(&self, py: Python<'_>, qubits: Vec) -> PyResult> { py.detach(|| self.inner.probabilities_with_readout(&qubits)) .map_err(qis_error_to_py_err) diff --git a/crates/binding-python/tests/circuit/test_decompose.py b/crates/binding-python/tests/circuit/test_decompose.py index 668edb6f..7833ecd8 100644 --- a/crates/binding-python/tests/circuit/test_decompose.py +++ b/crates/binding-python/tests/circuit/test_decompose.py @@ -13,6 +13,115 @@ """Decomposition tests for the Python circuit API.""" from cqlib import Circuit, Parameter +from cqlib.qis import DensityMatrix, Statevector +import numpy as np +import pytest + + +@pytest.mark.parametrize( + "width,has_operations", [(0, False), (1, False), (1, True), (2, True)] +) +@pytest.mark.parametrize("repetitions", [0, 1, 2]) +def test_decompose_preserves_global_phase_of_each_definition( + width, has_operations, repetitions +): + inner = Circuit(width) + inner.set_global_phase(0.37) + if has_operations: + inner.x(0) + if width == 2: + inner.cx(0, 1) + gate = inner.to_gate("Phased") + outer = Circuit(2) + outer.set_global_phase(0.12) + expected = Circuit(2) + expected.set_global_phase(0.12 + repetitions * 0.37) + for _ in range(repetitions): + outer.append_circuit_gate(gate, [1, 0][:width]) + if has_operations: + expected.x(1) + if width == 2: + expected.cx(1, 0) + flat = outer.decompose() + assert flat.global_phase.evaluate() == pytest.approx(0.12 + repetitions * 0.37) + for candidate in (outer, flat, flat.decompose()): + np.testing.assert_allclose( + candidate.to_matrix(), expected.to_matrix(), rtol=0, atol=1e-10 + ) + assert len(flat) == len(expected) + assert outer.global_phase.evaluate() == 0.12 + reference_state = expected.to_matrix()[:, 0] + np.testing.assert_allclose( + Statevector.from_circuit(outer).data, reference_state, rtol=0, atol=1e-10 + ) + np.testing.assert_allclose( + DensityMatrix.from_circuit(outer).data, + np.outer(reference_state, reference_state.conj()), + rtol=0, + atol=1e-10, + ) + + +@pytest.mark.parametrize("inverse", [False, True]) +def test_decompose_preserves_nested_symbolic_global_phase(inverse): + theta = Parameter("theta") + inner = Circuit([7, 2]) + inner.set_global_phase(theta / 2) + inner.rx(7, theta) + inner.cx(7, 2) + middle = Circuit([5, 1]) + middle.set_global_phase(theta.sin()) + middle.append_circuit_gate(inner.to_gate("Inner"), [5, 1], [theta + 1]) + outer = Circuit(2) + x = Parameter("x") + outer.set_global_phase(x / 3) + outer.append_circuit_gate(middle.to_gate("Middle"), [1, 0], [2 * x]) + if inverse: + outer = outer.inverse() + flat = outer.decompose() + symbolic = flat.to_symbolic_matrix() + for value in (0.37, -0.61): + phase = value / 3 + np.sin(2 * value) + (2 * value + 1) / 2 + expected = Circuit(2) + expected.set_global_phase(phase) + expected.rx(1, 2 * value + 1) + expected.cx(1, 0) + if inverse: + expected = expected.inverse() + phase = -phase + bindings = {"x": value} + reference = expected.to_matrix() + assert flat.global_phase.evaluate(bindings) == pytest.approx(phase) + for candidate in ( + flat.assign_parameters(bindings), + outer.assign_parameters(bindings).decompose(), + outer.assign_parameters(bindings), + ): + np.testing.assert_allclose( + candidate.to_matrix(), reference, rtol=0, atol=1e-10 + ) + np.testing.assert_allclose( + symbolic.evaluate(bindings), reference, rtol=0, atol=1e-10 + ) + + +def test_decompose_keeps_symbols_used_only_by_global_phase(): + inner = Circuit(1) + inner.set_global_phase(Parameter("theta")) + outer = Circuit(1) + outer.set_global_phase(Parameter("phi")) + outer.append_circuit_gate(inner.to_gate("PhaseOnly"), [0], [2 * Parameter("x")]) + flat = outer.decompose() + assert len(flat) == 0 + assert set(flat.used_symbols) == {"phi", "x"} + bindings = {"phi": 0.12, "x": 0.37} + assert flat.global_phase.evaluate(bindings) == pytest.approx(0.86) + np.testing.assert_allclose( + flat.assign_parameters(bindings).to_matrix(), + np.exp(0.86j) * np.eye(2), + rtol=0, + atol=1e-10, + ) def test_decompose_keeps_standard_gates_unchanged(): diff --git a/crates/binding-python/tests/compile/test_commutation.py b/crates/binding-python/tests/compile/test_commutation.py index 4eb5692c..4696ad8e 100644 --- a/crates/binding-python/tests/compile/test_commutation.py +++ b/crates/binding-python/tests/compile/test_commutation.py @@ -17,6 +17,7 @@ import pytest from cqlib.circuit import ( + Circuit, Instruction, Parameter, Qubit, @@ -64,6 +65,46 @@ def test_builtin_checker_proves_exact_and_global_phase_commutation(): assert "Parameter(" in repr(global_phase) +@pytest.mark.parametrize("order", [[0, 1], [1, 0], [7, 2]]) +@pytest.mark.parametrize("parameterized", [False, True]) +def test_circuit_gate_matrix_commutation_uses_positional_targets(order, parameterized): + body = Circuit(order) + body.rx(order[0], Parameter("theta") if parameterized else 0.37) + circuit = Circuit(2) + circuit.append_circuit_gate( + body.to_gate("first_rx"), [0, 1], [0.37] if parameterized else [] + ) + matrix_checker = CommutationChecker.with_config( + CommutationConfig(enable_rule_oracle=False, enable_matrix_fallback=True) + ) + for target in (0, 1): + z = operation(StandardGate.Z, [target]) + expected = None if target == 0 else Commutation.exact() + for lhs, rhs in ((circuit[0], z), (z, circuit[0])): + assert check_commutation(lhs, rhs) == expected + assert matrix_checker.check(lhs, rhs) == expected + + +def test_circuit_gate_matrix_commutation_preserves_global_phase_proof(): + body = Circuit([7, 2]) + body.x(7) + circuit = Circuit(2) + circuit.append_circuit_gate(body.to_gate("first_x"), [0, 1]) + z = operation(StandardGate.Z, [0]) + for lhs, rhs in ((circuit[0], z), (z, circuit[0])): + proof = check_commutation(lhs, rhs) + assert proof is not None and not proof.is_exact() + assert math.isclose(abs(proof.phase.evaluate()), math.pi, abs_tol=1e-12) + + +def test_circuit_gate_unresolved_matrix_does_not_produce_a_proof(): + body = Circuit(2) + body.rx(0, Parameter("theta")) + circuit = Circuit(2) + circuit.append_circuit_gate(body.to_gate("symbolic"), [0, 1], [Parameter("x")]) + assert check_commutation(circuit[0], operation(StandardGate.Z, [0])) is None + + def test_symbolic_parameters_are_preserved_for_algebraic_proofs(): lhs = operation(StandardGate.RZ(Parameter("a")), [0]) rhs = operation(StandardGate.RZ(Parameter("b")), [0]) diff --git a/crates/binding-python/tests/gates/test_circuit_gate.py b/crates/binding-python/tests/gates/test_circuit_gate.py index a18d8416..969ef952 100644 --- a/crates/binding-python/tests/gates/test_circuit_gate.py +++ b/crates/binding-python/tests/gates/test_circuit_gate.py @@ -20,12 +20,202 @@ - Parametric CircuitGate creation and usage """ +from itertools import permutations + +import numpy as np import pytest from cqlib.circuit import Circuit, Parameter from cqlib.circuit.gates import CircuitGate +def _embed_operation_matrix(matrix, targets, width): + """Independent basis-index embedding: first gate argument is the MSB.""" + result = np.zeros((2**width, 2**width), dtype=complex) + for column in range(2**width): + local_column = 0 + for target in targets: + local_column = 2 * local_column + ((column >> target) & 1) + for local_row in range(2 ** len(targets)): + row = column + for position, target in enumerate(targets): + bit = (local_row >> (len(targets) - 1 - position)) & 1 + row = (row & ~(1 << target)) | (bit << target) + result[row, column] = matrix[local_row, local_column] + return result + + +class TestCircuitGateOperationMatrix: + @staticmethod + def assert_operations(circuit, expected): + # Cover snapshots returned by both collection and indexing APIs. + for operation in (circuit.operations[0], circuit[0]): + np.testing.assert_allclose(operation.matrix(), expected, rtol=0, atol=1e-12) + + @pytest.mark.parametrize("order", [[0, 1], [1, 0], [7, 2]]) + def test_positional_matrix(self, order): + body = Circuit(order) + body.x(order[0]) + circuit = Circuit(2) + circuit.append_circuit_gate(body.to_gate("first_x"), [0, 1]) + self.assert_operations(circuit, np.kron([[0, 1], [1, 0]], np.eye(2))) + + def test_bound_parameter_and_repeated_application(self): + body = Circuit([7]) + body.rx(7, Parameter("theta")) + gate = body.to_gate("rx") + for angle in (0.37, -0.61, 0.37): + circuit = Circuit(1) + circuit.append_circuit_gate(gate, [0], [angle]) + c, s = np.cos(angle / 2), -1j * np.sin(angle / 2) + self.assert_operations(circuit, np.array([[c, s], [s, c]])) + assert gate.used_symbols == ["theta"] + + @pytest.mark.parametrize("order", list(permutations([7, 2, 5]))) + @pytest.mark.parametrize("targets", [[3, 0, 2], [2, 3, 0]]) + def test_three_qubit_permutations_and_spectator(self, order, targets): + body = Circuit(list(order)) + body.h(order[0]) + body.rz(order[1], 0.29) + body.cx(order[0], order[2]) + body.set_global_phase(0.17) + expected = np.zeros((8, 8), dtype=complex) + for column in range(8): + first, middle, last = column >> 2, (column >> 1) & 1, column & 1 + for output_first in (0, 1): + row = (output_first << 2) | (middle << 1) | (last ^ output_first) + expected[row, column] = ( + (-1) ** (first * output_first) + / np.sqrt(2) + * np.exp(1j * (0.17 + (2 * middle - 1) * 0.29 / 2)) + ) + circuit = Circuit(4) + circuit.append_circuit_gate(body.to_gate("three"), targets) + self.assert_operations(circuit, expected) + full = _embed_operation_matrix(expected, targets, 4) + np.testing.assert_allclose(circuit.to_matrix(), full, rtol=0, atol=1e-12) + np.testing.assert_allclose( + circuit.to_symbolic_matrix().evaluate(), full, rtol=0, atol=1e-12 + ) + np.testing.assert_allclose( + circuit.decompose().to_matrix(), full, rtol=0, atol=1e-12 + ) + + @pytest.mark.parametrize("nested", [False, True]) + @pytest.mark.parametrize("inverse", [False, True]) + def test_signature_phase_and_inverse(self, nested, inverse): + theta, phi = Parameter("theta"), Parameter("phi") + body = Circuit([7, 2]) + body.rx(7, theta) + body.rz(2, phi) + body.cx(7, 2) + body.set_global_phase(theta + phi / 2) + gate = CircuitGate.with_signature( + "declared", body.to_gate("definition").circuit, ["unused", "phi", "theta"] + ) + arguments = [12.0, -0.21, 0.37] + if nested: + middle = Circuit([5, 1]) + middle.append_circuit_gate(gate, [5, 1], [0, phi, theta]) + gate = CircuitGate.with_signature( + "nested", middle.to_gate("middle").circuit, ["unused", "phi", "theta"] + ) + circuit = Circuit(2) + circuit.append_circuit_gate(gate, [1, 0], arguments) + if inverse: + circuit = circuit.inverse() + c, s = np.cos(0.37 / 2), -1j * np.sin(0.37 / 2) + rx = np.array([[c, s], [s, c]]) + rz = np.diag(np.exp(1j * np.array([0.21, -0.21]) / 2)) + cx = np.eye(4)[[0, 1, 3, 2]] + expected = np.exp(1j * (0.37 - 0.21 / 2)) * cx @ np.kron(rx, rz) + if inverse: + expected = expected.conj().T + self.assert_operations(circuit, expected) + np.testing.assert_allclose(circuit.to_matrix(), expected, rtol=0, atol=1e-12) + np.testing.assert_allclose( + circuit.to_symbolic_matrix().evaluate(), expected, rtol=0, atol=1e-12 + ) + + def test_zero_qubit_phase(self): + body = Circuit(0) + body.set_global_phase(Parameter("theta")) + circuit = Circuit(0) + circuit.append_circuit_gate(body.to_gate("phase"), [], [0.37]) + expected = np.array([[np.exp(0.37j)]]) + self.assert_operations(circuit, expected) + np.testing.assert_allclose(circuit.to_matrix(), expected, rtol=0, atol=1e-12) + + def test_unresolved_arguments_keep_existing_error_behavior(self): + body = Circuit(1) + body.rx(0, Parameter("theta")) + circuit = Circuit(1) + circuit.append_circuit_gate(body.to_gate("unresolved"), [0], [Parameter("x")]) + with pytest.raises(ValueError): + circuit.operations[0].matrix() + with pytest.raises(ValueError): + circuit[0].matrix() + + +class TestCircuitGateStorageOrder: + @pytest.mark.parametrize("order", [[7, 2], [1, 0], [2, 7], [0, 1]]) + @pytest.mark.parametrize("targets", [[0, 1], [1, 0]]) + def test_matrix_uses_internal_storage_order(self, order, targets): + inner = Circuit(order) + inner.x(order[0]) + outer = Circuit(2) + outer.append_circuit_gate(inner.to_gate("StoredOrderX"), targets) + expected = Circuit(2) + expected.x(targets[0]) + for output_order in (None, [1, 0]): + reference = expected.to_matrix(output_order) + np.testing.assert_allclose( + outer.to_matrix(output_order), reference, rtol=0, atol=1e-10 + ) + np.testing.assert_allclose( + outer.to_symbolic_matrix(output_order).evaluate(), + reference, + rtol=0, + atol=1e-10, + ) + np.testing.assert_allclose( + outer.decompose().to_matrix(output_order), reference, rtol=0, atol=1e-10 + ) + assert outer.to_matrix()[1 << targets[0], 0] == 1 + + @pytest.mark.parametrize("middle_order", [None, [0, 1], [5, 1]]) + @pytest.mark.parametrize("inverse", [False, True]) + def test_nested_parameterized_matrix_and_inverse(self, middle_order, inverse): + theta = Parameter("theta") + body = Circuit([7, 2]) + body.rx(7, theta) + body.cx(7, 2) + if middle_order is not None: + gate = body.to_gate("Inner") + body = Circuit(middle_order) + body.append_circuit_gate(gate, middle_order, [theta]) + outer = Circuit(2) + outer.append_circuit_gate(body.to_gate("Body"), [0, 1], [Parameter("x")]) + if inverse: + outer = outer.inverse() + symbolic = outer.to_symbolic_matrix() + for angle in (0.37, -0.61): + expected = Circuit(2) + expected.rx(0, angle) + expected.cx(0, 1) + if inverse: + expected = expected.inverse() + reference = expected.to_matrix() + bound = outer.assign_parameters({"x": angle}) + np.testing.assert_allclose(bound.to_matrix(), reference, rtol=0, atol=1e-10) + np.testing.assert_allclose( + symbolic.evaluate({"x": angle}), reference, rtol=0, atol=1e-10 + ) + np.testing.assert_allclose( + bound.decompose().to_matrix(), reference, rtol=0, atol=1e-10 + ) + + class TestCircuitGateCreation: """Test CircuitGate creation from circuits""" diff --git a/crates/binding-python/tests/gates/test_gphase.py b/crates/binding-python/tests/gates/test_gphase.py new file mode 100644 index 00000000..2e9e8df8 --- /dev/null +++ b/crates/binding-python/tests/gates/test_gphase.py @@ -0,0 +1,89 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 2026 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Zero-qubit phase and controlled-phase matrix contracts.""" + +import numpy as np +import pytest + +from cqlib import Circuit, Parameter +from cqlib.circuit.gates import GPhase + + +@pytest.mark.parametrize("theta", [0.0, 0.37, -0.51, np.pi]) +@pytest.mark.parametrize("controls", [None, 0, 1, 2]) +def test_gphase_matrix_metadata_and_inverse(theta, controls): + base = GPhase(theta) + gate = base if controls is None else base.control(controls) + count = controls or 0 + assert gate.num_qubits == count + assert gate.num_ctrl_qubits == count + assert gate.num_params == 1 + expected = np.eye(1 << count, dtype=complex) + expected[-1, -1] = np.exp(1j * theta) + actual = gate.matrix() + assert actual.shape == expected.shape + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-12) + inverse = gate.inverse().matrix() + assert inverse.shape == expected.shape + np.testing.assert_allclose(inverse, expected.conj().T, rtol=0, atol=1e-12) + + +@pytest.mark.parametrize("num_qubits", [0, 1, 3]) +@pytest.mark.parametrize("symbolic", [False, True]) +def test_gphase_markers_accumulate_in_circuit_matrices(num_qubits, symbolic): + circuit = Circuit(num_qubits) + if num_qubits: + circuit.h(0) + baseline = circuit.to_matrix() + circuit.append_gate(GPhase(-0.21), []) + circuit.append_gate(GPhase(Parameter("theta") if symbolic else 0.37), []) + circuit.append_gate(GPhase(0.12), []) + circuit.set_global_phase(0.17) + expected = np.exp(1j * (-0.21 + 0.37 + 0.12 + 0.17)) * baseline + if not symbolic: + actual = circuit.to_matrix() + assert actual.shape == expected.shape + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-12) + expression = circuit.to_symbolic_matrix() + assert expression.shape == expected.shape + np.testing.assert_allclose( + expression.evaluate({"theta": 0.37}), expected, rtol=0, atol=1e-12 + ) + + +@pytest.mark.parametrize("controls", [[], [2], [2, 0]]) +@pytest.mark.parametrize("order", [[0, 1, 2], [2, 0, 1]]) +@pytest.mark.parametrize("symbolic", [False, True]) +def test_gphase_controlled_circuit_matrices(controls, order, symbolic): + circuit = Circuit(3) + circuit.h(1) + theta = Parameter("theta") if symbolic else 0.37 + circuit.append_mc_gate(GPhase(theta).control(len(controls)), controls) + spectator_bit = order.index(1) + mask = sum(1 << order.index(control) for control in controls) + expected = np.zeros((8, 8), dtype=complex) + for row in range(8): + for col in range(8): + if (row ^ col) & ~(1 << spectator_bit) == 0: + sign = -1 if row & col & (1 << spectator_bit) else 1 + phase = np.exp(0.37j) if row & mask == mask else 1 + expected[row, col] = phase * sign / np.sqrt(2) + if not symbolic: + actual = circuit.to_matrix(order) + assert actual.shape == (8, 8) + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-12) + expression = circuit.to_symbolic_matrix(order) + assert expression.shape == (8, 8) + np.testing.assert_allclose( + expression.evaluate({"theta": 0.37}), expected, rtol=0, atol=1e-12 + ) diff --git a/crates/binding-python/tests/gates/test_standard.py b/crates/binding-python/tests/gates/test_standard.py index 4c7770fe..a0d64e27 100644 --- a/crates/binding-python/tests/gates/test_standard.py +++ b/crates/binding-python/tests/gates/test_standard.py @@ -396,12 +396,12 @@ def test_phase_gate(self): assert np.allclose(mat, expected), "Phase gate matrix mismatch" def test_global_phase_gate(self): - """GPhase(θ) = e^(iθ) I (2x2)""" + """The zero-qubit GPhase matrix contains one scalar phase.""" theta = np.pi / 4 mat = GPhase(theta).matrix() - # GPhase is 1-qubit gate with 1 parameter, matrix is 2x2 - expected = np.exp(1j * theta) * np.eye(2) - assert np.allclose(mat, expected), "GPhase matrix mismatch" + expected = np.array([[np.exp(1j * theta)]]) + assert mat.shape == (1, 1) + np.testing.assert_allclose(mat, expected, rtol=0, atol=1e-12) def test_u_gate_general(self): """U gate is general single-qubit gate""" diff --git a/crates/binding-python/tests/ir/test_qasm3_gphase.py b/crates/binding-python/tests/ir/test_qasm3_gphase.py new file mode 100644 index 00000000..dca06a10 --- /dev/null +++ b/crates/binding-python/tests/ir/test_qasm3_gphase.py @@ -0,0 +1,85 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 2026 +# +# Licensed under the Apache License, Version 2.0. See LICENSE.txt or +# http://www.apache.org/licenses/LICENSE-2.0. + +"""Global phase must survive serialization of reusable gate definitions.""" + +import numpy as np +import pytest + +from cqlib import Circuit, Parameter +from cqlib.circuit import UnitaryGate +from cqlib.circuit.gates import GPhase +from cqlib.ir import qasm3 +from cqlib.qis import Statevector + + +def assert_round_trip(circuit, expected, bindings=None): + output = qasm3.dumps(circuit) + loaded = qasm3.loads(output) + original = circuit.assign_parameters(bindings) if bindings else circuit + loaded = loaded.assign_parameters(bindings) if bindings else loaded + for candidate in [original, loaded]: + np.testing.assert_allclose(candidate.to_matrix(), expected, rtol=0, atol=1e-12) + np.testing.assert_allclose( + Statevector.from_circuit(candidate).data, + expected[:, 0], + rtol=0, + atol=1e-12, + ) + return output + + +@pytest.mark.parametrize("symbolic", [False, True]) +@pytest.mark.parametrize("explicit_phase", [False, True]) +@pytest.mark.parametrize("nested", [False, True]) +def test_nested_gate_global_phase_round_trip(symbolic, explicit_phase, nested): + theta = Parameter("theta") if symbolic else 0.37 + inner = Circuit([12]) + inner.h(12) + inner.set_global_phase(theta) + if explicit_phase: + inner.append_gate(GPhase(-0.21), []) + gate = inner.to_gate("fragment") + if nested: + middle = Circuit([5]) + middle.set_global_phase(0.11) + middle.append_circuit_gate(gate, [5], [theta] if symbolic else []) + gate = middle.to_gate("middle") + outer = Circuit(1) + outer.set_global_phase(-0.07) + outer.append_circuit_gate(gate, [0], [Parameter("outer_theta")] if symbolic else []) + phase = 0.37 - 0.07 + (0.11 if nested else 0) - (0.21 if explicit_phase else 0) + expected = np.exp(1j * phase) * np.array([[1, 1], [1, -1]]) / np.sqrt(2) + assert_round_trip(outer, expected, {"outer_theta": 0.37} if symbolic else None) + + +def test_phase_only_parameter_and_repeated_gate_calls(): + inner = Circuit(1) + inner.h(0) + inner.set_global_phase(Parameter("theta")) + gate = inner.to_gate("fragment") + outer = Circuit(1) + outer.append_circuit_gate(gate, [0], [0.37]) + outer.append_circuit_gate(gate, [0], [-0.51]) + output = assert_round_trip(outer, np.exp(-0.14j) * np.eye(2)) + assert output.count("gphase(theta);") == 1 + + +def test_unitary_gate_definition_global_phase_round_trip(): + inner = Circuit(1) + inner.h(0) + inner.set_global_phase(0.37) + inner.append_gate(GPhase(-0.21), []) + expected = np.exp(0.16j) * np.array([[1, 1], [1, -1]]) / np.sqrt(2) + gate = ( + UnitaryGate("phase_h", 1) + .with_matrix(expected) + .with_circuit(inner.to_gate("definition").circuit) + ) + outer = Circuit(1) + outer.append_unitary_gate(gate, [0]) + assert_round_trip(outer, expected) diff --git a/crates/binding-python/tests/parameter/test_substitution.py b/crates/binding-python/tests/parameter/test_substitution.py new file mode 100644 index 00000000..cb650b98 --- /dev/null +++ b/crates/binding-python/tests/parameter/test_substitution.py @@ -0,0 +1,107 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 2026 +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +import numpy as np +import pytest + +from cqlib import Circuit, Parameter + + +@pytest.mark.parametrize("reverse", [False, True]) +def test_substitute_swaps_symbols_using_the_same_dictionary(reverse): + x, y = Parameter("x"), Parameter("y") + original = x + 2 * y + bindings = dict(reversed([("x", y), ("y", x)])) if reverse else {"x": y, "y": x} + + for _ in range(100): + result = original.substitute(bindings) + assert set(result.symbols) == {"x", "y"} + assert result.evaluate({"x": 1.0, "y": 3.0}) == 5.0 + + assert original.evaluate({"x": 1.0, "y": 3.0}) == 7.0 + assert bindings["x"] == y + assert bindings["y"] == x + + +@pytest.mark.parametrize( + "case", ["chain", "self", "cycle", "multiple_keys", "independent"] +) +def test_substitute_does_not_process_inserted_expressions(case): + x, y, z = (Parameter(name) for name in ("x", "y", "z")) + original, bindings, expected = { + "chain": (x, {"x": y + 1, "y": Parameter(2)}, 11.0), + "self": (x, {"x": x + 1}, 3.0), + "cycle": (x + 2 * y + 3 * z, {"x": y, "y": z, "z": x}, 56.0), + "multiple_keys": (x, {"x": y + z, "y": Parameter(1), "z": Parameter(2)}, 30.0), + "independent": (x + y, {"x": z + 1, "y": Parameter(2)}, 23.0), + }[case] + assert ( + original.substitute(bindings).evaluate({"x": 2.0, "y": 10.0, "z": 20.0}) + == expected + ) + + +@pytest.mark.parametrize("method", ["replace", "substitute"]) +@pytest.mark.parametrize("kind", ["polynomial", "function_polynomial", "nested"]) +def test_substitution_reaches_polynomials_and_nested_functions(method, kind): + x, y = Parameter("x"), Parameter("y") + base = x.sin() + original = { + "polynomial": x * x + x + 1, + "function_polynomial": base * base + 2 * base + 1, + "nested": ((x * x + x + 1) / (x + 2)).cos() + x ** (y + 1), + }[kind] + result = ( + original.replace("x", y) + if method == "replace" + else original.substitute({"x": y}) + ) + assert set(result.symbols) == {"y"} + assert result.evaluate({"y": 0.5}) == pytest.approx( + original.evaluate({"x": 0.5, "y": 0.5}), abs=1e-12 + ) + + +def test_substitute_preserves_unbound_symbols_constants_and_empty_bindings(): + x, y = Parameter("x"), Parameter("y") + original = x + y + for bindings in ({}, {"missing": Parameter(9)}): + result = original.substitute(bindings) + assert set(result.symbols) == {"x", "y"} + assert result.evaluate({"x": 2.0, "y": 3.0}) == 5.0 + partial = original.substitute({"x": Parameter(2)}) + assert partial.symbols == ["y"] + assert partial.evaluate({"y": 3.0}) == 5.0 + assert Parameter(0.25).substitute({"x": y}).evaluate() == 0.25 + + +@pytest.mark.parametrize( + "coefficient", [1e-16, np.nextafter(0.0, 1.0), 1.0 + np.finfo(float).eps] +) +def test_substitution_preserves_affine_coefficients(coefficient): + x, y = Parameter("x"), Parameter("y") + original = x * float(coefficient) + for result in (original.replace("x", y), original.substitute({"x": y})): + assert result.evaluate({"y": 1.0}) == coefficient + + +def test_numeric_circuit_assignment_remains_equivalent(): + x, y = Parameter("x"), Parameter("y") + circuit = Circuit(1) + circuit.rx(0, x + 2 * y) + expected = Circuit(1) + expected.rx(0, 7.0) + for bound in ( + circuit.assign_parameters({"x": 1.0, "y": 3.0}), + circuit.assign_parameters({"x": 1.0}).assign_parameters({"y": 3.0}), + ): + np.testing.assert_allclose(bound.to_matrix(), expected.to_matrix(), atol=1e-12) diff --git a/crates/binding-python/tests/qis/test_controlled_gphase.py b/crates/binding-python/tests/qis/test_controlled_gphase.py new file mode 100644 index 00000000..2b29f446 --- /dev/null +++ b/crates/binding-python/tests/qis/test_controlled_gphase.py @@ -0,0 +1,226 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 2026 +# +# Licensed under the Apache License, Version 2.0. See LICENSE.txt or +# http://www.apache.org/licenses/LICENSE-2.0. + +"""Controlled zero-target phase evolution in the numerical simulators.""" + +import numpy as np +import pytest + +from cqlib import Circuit, Parameter +from cqlib.circuit import StandardGate +from cqlib.circuit.gates import GPhase +from cqlib.device import NoiseModel, SingleQubitNoise +from cqlib.qis import DensityMatrix, DensityMatrixNoise, Statevector + + +SIMULATORS = [Statevector, DensityMatrix, DensityMatrixNoise] + + +def assert_state(state, amplitudes): + if isinstance(state, Statevector): + actual, expected = state.data, amplitudes + else: + actual = state.state if isinstance(state, DensityMatrixNoise) else state.data + expected = np.outer(amplitudes, amplitudes.conj()) + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-12) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +@pytest.mark.parametrize("num_qubits", [0, 1, 3]) +def test_controlled_gphase_preserves_uncontrolled_phase(simulator, num_qubits): + circuit = Circuit(num_qubits) + expected = np.zeros(1 << num_qubits, dtype=complex) + expected[0] = 1 + if num_qubits: + circuit.h(0) + expected[:2] = 1 / np.sqrt(2) + circuit.append_gate(GPhase(0.37), []) + circuit.set_global_phase(-0.21) + expected *= np.exp(0.16j) + assert_state(simulator.from_circuit(circuit), expected) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +@pytest.mark.parametrize("theta", [0.0, 0.37, -0.51, np.pi]) +@pytest.mark.parametrize("initial", ["zero", "one", "plus"]) +@pytest.mark.parametrize("incremental", [False, True]) +def test_controlled_gphase_single_qubit(simulator, theta, initial, incremental): + circuit = Circuit(1) + expected = np.array([1, 0], dtype=complex) + if initial == "one": + circuit.x(0) + expected = np.array([0, 1], dtype=complex) + elif initial == "plus": + circuit.h(0) + expected = np.ones(2, dtype=complex) / np.sqrt(2) + if incremental: + state = simulator.from_circuit(circuit) + circuit = Circuit(1) + circuit.append_mc_gate(GPhase(theta).control(1), [0]) + expected[1] *= np.exp(1j * theta) + if incremental: + state.apply_circuit(circuit) + else: + state = simulator.from_circuit(circuit) + assert_state(state, expected) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +@pytest.mark.parametrize("labels", [[0, 1, 2], [12, 5, 31]]) +@pytest.mark.parametrize("positions", [[2], [2, 0], [0, 2]]) +def test_controlled_gphase_labels_spectators_and_inverse(simulator, labels, positions): + circuit = Circuit(labels) + for qubit in labels: + circuit.h(qubit) + circuit.s(labels[1]) + controls = [labels[position] for position in positions] + phase = Circuit(labels) + gate = GPhase(0.37).control(len(controls)) + phase.append_mc_gate(gate, controls) + incremental = simulator.from_circuit(circuit) + incremental.apply_circuit(phase) + circuit.append_mc_gate(gate, controls) + direct = simulator.from_circuit(circuit) + indices = np.arange(8) + baseline = np.where(indices & 2, 1j, 1) / np.sqrt(8) + mask = sum(1 << position for position in positions) + expected = baseline * np.where(indices & mask == mask, np.exp(0.37j), 1) + assert_state(direct, expected) + assert_state(incremental, expected) + inverse = Circuit(labels) + inverse.append_mc_gate(gate.inverse(), controls) + incremental.apply_circuit(inverse) + assert_state(incremental, baseline) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +def test_controlled_gphase_parameter_resolution(simulator): + circuit = Circuit(1) + circuit.h(0) + circuit.append_mc_gate(GPhase(Parameter("theta")).control(1), [0]) + with pytest.raises(ValueError, match="(?i)symbolic"): + simulator.from_circuit(circuit) + state = simulator(1) + with pytest.raises(ValueError, match="(?i)symbolic"): + state.apply_circuit(circuit) + bound = circuit.assign_parameters({"theta": -0.51}) + expected = np.array([1, np.exp(-0.51j)]) / np.sqrt(2) + assert_state(simulator.from_circuit(bound), expected) + + +@pytest.mark.parametrize("empty_model", [False, True]) +def test_controlled_gphase_noisy_simulator_without_errors(empty_model): + circuit = Circuit(1) + circuit.h(0) + circuit.append_mc_gate(GPhase(0.37).control(1), [0]) + model = NoiseModel() if empty_model else None + state = DensityMatrixNoise.from_circuit(circuit, model) + assert_state(state, np.array([1, np.exp(0.37j)]) / np.sqrt(2)) + + +@pytest.mark.parametrize("incremental", [False, True]) +def test_controlled_gphase_noise_follows_unitary_once_with_original_labels(incremental): + labels = [12, 5] + circuit = Circuit(labels) + circuit.x(12) + circuit.h(5) + model = NoiseModel() + probability = 0.23 + model.add_single_qubit_error( + StandardGate.GPhase, 5, SingleQubitNoise.bit_flip(probability) + ) + model.add_single_qubit_error( + StandardGate.GPhase, 12, SingleQubitNoise.bit_flip(1.0) + ) + if incremental: + state = DensityMatrixNoise.from_circuit(circuit, model) + circuit = Circuit(labels) + circuit.append_mc_gate(GPhase(0.37).control(1), [5]) + if incremental: + state.apply_circuit(circuit) + else: + state = DensityMatrixNoise.from_circuit(circuit, model) + amplitudes = np.array([0, 1, 0, np.exp(0.37j)]) / np.sqrt(2) + rho = np.outer(amplitudes, amplitudes.conj()) + flipped = np.arange(4) ^ 2 + expected = (1 - probability) * rho + probability * rho[np.ix_(flipped, flipped)] + np.testing.assert_allclose(state.state, expected, rtol=0, atol=1e-12) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +@pytest.mark.parametrize("num_qubits", [0, 1, 3]) +@pytest.mark.parametrize("theta", [0.0, 0.37, -0.51, np.pi]) +@pytest.mark.parametrize("incremental", [False, True]) +def test_zero_control_gphase_evolution(simulator, num_qubits, theta, incremental): + preparation = Circuit(num_qubits) + expected = np.zeros(1 << num_qubits, dtype=complex) + expected[0] = 1 + if num_qubits: + preparation.h(0) + expected[:2] = 1 / np.sqrt(2) + circuit = Circuit(num_qubits) if incremental else preparation + circuit.append_mc_gate(GPhase(theta).control(0), []) + circuit.append_mc_gate(GPhase(-0.21).control(0), []) + circuit.set_global_phase(0.11) + expected *= np.exp(1j * (theta - 0.10)) + if incremental: + state = simulator.from_circuit(preparation) + state.apply_circuit(circuit) + else: + state = simulator.from_circuit(circuit) + assert_state(state, expected) + inverse = Circuit(num_qubits) + inverse.append_mc_gate(GPhase(theta - 0.10).control(0).inverse(), []) + state.apply_circuit(inverse) + assert_state(state, expected * np.exp(-1j * (theta - 0.10))) + + +@pytest.mark.parametrize("simulator", SIMULATORS) +def test_zero_control_gphase_parameter_resolution(simulator): + circuit = Circuit([12, 5]) + circuit.h(5) + circuit.append_mc_gate(GPhase(Parameter("theta")).control(0), []) + with pytest.raises(ValueError, match="(?i)symbolic"): + simulator.from_circuit(circuit) + with pytest.raises(ValueError, match="(?i)symbolic"): + simulator(2).apply_circuit(circuit) + bound = circuit.assign_parameters({"theta": -0.51}) + expected = np.array([1, 0, 1, 0], dtype=complex) * np.exp(-0.51j) / np.sqrt(2) + assert_state(simulator.from_circuit(bound), expected) + + +@pytest.mark.parametrize("empty_model", [False, True]) +@pytest.mark.parametrize("incremental", [False, True]) +def test_zero_control_gphase_noise_matches_plain_gate(empty_model, incremental): + model = NoiseModel() + if not empty_model: + model.add_single_qubit_error( + StandardGate.GPhase, 12, SingleQubitNoise.bit_flip(1.0) + ) + model.add_single_qubit_error( + StandardGate.H, 5, SingleQubitNoise.depolarizing(0.23) + ) + preparation = Circuit([12, 5]) + preparation.x(12) + preparation.h(5) + results = [] + for wrapped in [False, True]: + circuit = Circuit([12, 5]) + if not incremental: + circuit.x(12) + circuit.h(5) + if wrapped: + circuit.append_mc_gate(GPhase(0.37).control(0), []) + else: + circuit.append_gate(GPhase(0.37), []) + if incremental: + state = DensityMatrixNoise.from_circuit(preparation, model) + state.apply_circuit(circuit) + else: + state = DensityMatrixNoise.from_circuit(circuit, model) + results.append(state.state) + np.testing.assert_allclose(results[0], results[1], rtol=0, atol=1e-12) diff --git a/crates/binding-python/tests/qis/test_density_matrix_noise.py b/crates/binding-python/tests/qis/test_density_matrix_noise.py index 22184776..21687db0 100644 --- a/crates/binding-python/tests/qis/test_density_matrix_noise.py +++ b/crates/binding-python/tests/qis/test_density_matrix_noise.py @@ -101,6 +101,35 @@ def test_probabilities_with_readout(self): assert len(probs) == 4 assert abs(sum(probs) - 1.0) < 1e-10 + @pytest.mark.parametrize("ids", [[0, 1], [1, 0], [2, 0], [9, 2]]) + @pytest.mark.parametrize("bind_with_apply", [False, True]) + def test_readout_uses_circuit_ids_at_storage_positions(self, ids, bind_with_apply): + circuit = Circuit(ids) + receipt = circuit.measure(ids[0]) + model = NoiseModel() + model.add_readout_error(ids[0], ReadoutError(0.0, 1.0)) + if bind_with_apply: + sim = DensityMatrixNoise(2, model) + sim.apply_circuit(circuit) + else: + sim = DensityMatrixNoise.from_circuit(circuit, model) + + assert sim.probabilities_with_readout([0]) == [0.0, 1.0, 0.0, 0.0] + assert sim.probabilities_with_readout([1]) == [1.0, 0.0, 0.0, 0.0] + assert sim.probabilities_with_readout([1, 0]) == [0.0, 1.0, 0.0, 0.0] + assert sim.probabilities_with_readout([]) == sim.probabilities() + assert sim.probabilities() == [1.0, 0.0, 0.0, 0.0] + marginal = { + o.to_bitstring(1): p for o, p in sim.probs_with_readout(receipt).items() + } + assert marginal == {"1": 1.0} + + @pytest.mark.parametrize("model", [None, NoiseModel()]) + def test_readout_rejects_out_of_range_positions_without_errors(self, model): + sim = DensityMatrixNoise.from_circuit(Circuit([9, 2]), model) + with pytest.raises(IndexError): + sim.probabilities_with_readout([9]) + class TestDensityMatrixNoiseSingleQubitGates: """Test single-qubit gate operations with noise.""" diff --git a/crates/binding-python/tests/qis/test_density_matrix_targets.py b/crates/binding-python/tests/qis/test_density_matrix_targets.py new file mode 100644 index 00000000..061f03fc --- /dev/null +++ b/crates/binding-python/tests/qis/test_density_matrix_targets.py @@ -0,0 +1,141 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 2026 +# +# Licensed under the Apache License, Version 2.0. See LICENSE.txt or +# http://www.apache.org/licenses/LICENSE-2.0. + +"""Target validation must precede any density matrix mutation.""" + +import numpy as np +import pytest + +from cqlib.qis import DensityMatrix, DensityMatrixNoise + + +def initial_state(): + amplitudes = np.arange(1, 9) + 1j * np.arange(8, 0, -1) + return amplitudes / np.linalg.norm(amplitudes) + + +def householder(dim): + vector = np.arange(1, dim + 1, dtype=float) + return (np.eye(dim) - 2 * np.outer(vector, vector) / (vector @ vector)).astype( + complex + ) + + +def state_data(state): + return state.state if isinstance(state, DensityMatrixNoise) else state.data + + +@pytest.mark.parametrize("simulator", [DensityMatrix, DensityMatrixNoise]) +@pytest.mark.parametrize("target", [0, 1, 2]) +@pytest.mark.parametrize( + "method,params", + [ + ("apply_cx", []), + ("apply_cz", []), + ("apply_cy", []), + ("apply_swap", []), + ("apply_rxx", [0.37]), + ("apply_ryy", [0.37]), + ("apply_rzz", [0.37]), + ("apply_rzx", [0.37]), + ("apply_crx", [0.37]), + ("apply_cry", [0.37]), + ("apply_crz", [0.37]), + ("apply_fsim", [0.37, -0.21]), + ], +) +def test_duplicate_two_qubit_targets(simulator, target, method, params): + state = simulator(3) + state.apply_h(0) + state.apply_s(0) + before = state_data(state).copy() + with pytest.raises(ValueError, match=rf"(?i)duplicate qubit index {target}\b"): + getattr(state, method)(target, target, *params) + np.testing.assert_array_equal(state_data(state), before) + + +@pytest.mark.parametrize("simulator", [DensityMatrix, DensityMatrixNoise]) +@pytest.mark.parametrize("targets", [(0, 0, 2), (0, 2, 0), (2, 0, 0)]) +def test_duplicate_ccx_targets(simulator, targets): + state = simulator(3) + state.apply_h(0) + before = state_data(state).copy() + with pytest.raises(ValueError, match=r"(?i)duplicate qubit index 0\b"): + state.apply_ccx(*targets) + np.testing.assert_array_equal(state_data(state), before) + + +@pytest.mark.parametrize("target", [0, 1, 2]) +def test_duplicate_custom_gate_targets(target): + state = DensityMatrix.from_state(3, initial_state()) + before = state.data.copy() + with pytest.raises(ValueError, match=rf"(?i)duplicate qubit index {target}\b"): + state.apply_double_qubits_gate(target, target, householder(4)) + np.testing.assert_array_equal(state.data, before) + + +@pytest.mark.parametrize( + "targets,duplicate", + [ + ([0, 0], 0), + ([1, 1], 1), + ([2, 2], 2), + ([0, 2, 0], 0), + ([0, 0, 2], 0), + ([2, 0, 0], 0), + ], +) +def test_duplicate_kraus_targets(targets, duplicate): + state = DensityMatrix.from_state(3, initial_state()) + before = state.data.copy() + with pytest.raises(ValueError, match=rf"(?i)duplicate qubit index {duplicate}\b"): + state.apply_kraus(targets, [householder(1 << len(targets)).ravel()]) + np.testing.assert_array_equal(state.data, before) + + +def embed_operator(operator, targets, num_qubits): + """Build the full operator independently, preserving target argument order.""" + dim = 1 << num_qubits + result = np.zeros((dim, dim), dtype=complex) + mask = sum(1 << target for target in targets) + for column in range(dim): + local_column = 0 + for target in targets: + local_column = (local_column << 1) | ((column >> target) & 1) + for local_row in range(1 << len(targets)): + row = column & ~mask + for offset, target in enumerate(reversed(targets)): + row |= ((local_row >> offset) & 1) << target + result[row, column] = operator[local_row, local_column] + return result + + +@pytest.mark.parametrize("targets", [[0, 1], [1, 0], [0, 2], [2, 0]]) +@pytest.mark.parametrize("method", ["custom", "kraus", "unitary", "mixed_channel"]) +def test_distinct_targets_preserve_matrix_order(targets, method): + amplitudes = initial_state() + state = DensityMatrix.from_state(3, amplitudes) + operator = householder(4) + full_operator = embed_operator(operator, targets, 3) + evolved = full_operator @ amplitudes + expected = np.outer(evolved, evolved.conj()) + if method == "custom": + state.apply_double_qubits_gate(*targets, operator) + elif method == "unitary": + state.apply_unitary_gate(targets, operator) + elif method == "kraus": + state.apply_kraus(targets, [operator.ravel()]) + else: + state.apply_kraus( + targets, + [ + np.sqrt(0.7) * operator.ravel(), + np.sqrt(0.3) * np.eye(4, dtype=complex).ravel(), + ], + ) + expected = 0.7 * expected + 0.3 * np.outer(amplitudes, amplitudes.conj()) + np.testing.assert_allclose(state.data, expected, rtol=0, atol=1e-12) diff --git a/crates/binding-python/tests/qis/test_hamiltonian.py b/crates/binding-python/tests/qis/test_hamiltonian.py index 1abcf710..e57894ea 100644 --- a/crates/binding-python/tests/qis/test_hamiltonian.py +++ b/crates/binding-python/tests/qis/test_hamiltonian.py @@ -12,10 +12,113 @@ import pytest import math +import warnings import numpy as np from cqlib.qis import Hamiltonian, PauliString, Phase, TrotterMode +@pytest.fixture(params=["from_list", "add_term", "scale"]) +def hamiltonian_with_coefficient(request): + def make(pauli, coefficient): + if request.param == "from_list": + return Hamiltonian.from_list([(pauli, coefficient)]) + if request.param == "add_term": + h = Hamiltonian(1) + assert h.add_term(pauli, coefficient) is None + else: + h = Hamiltonian.from_pauli(pauli) + assert h.scale(coefficient) is None + return h + + return make + + +@pytest.mark.parametrize("scalar_type", [complex, np.complex128, np.complex64]) +@pytest.mark.parametrize( + "label, coefficient", [("iX", -1j), ("X", 1.25 - 2.5j), ("X", 2 + 0j)] +) +def test_hamiltonian_complex_scalar_coefficients( + hamiltonian_with_coefficient, scalar_type, label, coefficient +): + # In particular, (-i) * (iX) must be X rather than the zero operator. + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + h = hamiltonian_with_coefficient( + PauliString.from_str(label), scalar_type(coefficient) + ) + assert h.terms[0][1] == coefficient + assert not caught + phase = 1j if label == "iX" else 1 + expected = np.array([[0, coefficient * phase], [coefficient * phase, 0]]) + np.testing.assert_array_equal(h.to_matrix(), expected) + + +@pytest.mark.parametrize( + "coefficient, expected", + [ + (2, 2), + (2.5, 2.5), + (np.int64(2), 2), + (np.float32(2.5), 2.5), + (np.float64(2.5), 2.5), + ((1.25, -2.5), 1.25 - 2.5j), + ((np.float32(1.25), np.float64(-2.5)), 1.25 - 2.5j), + ], +) +def test_hamiltonian_real_and_tuple_coefficients( + hamiltonian_with_coefficient, coefficient, expected +): + with warnings.catch_warnings(): + warnings.simplefilter("error") + h = hamiltonian_with_coefficient(PauliString.from_str("X"), coefficient) + assert h.terms[0][1] == expected + np.testing.assert_array_equal(h.to_matrix(), [[0, expected], [expected, 0]]) + + +@pytest.mark.parametrize("coefficient", ["1+2j", object(), (), (1,), (1, 2, 3)]) +def test_hamiltonian_rejects_invalid_coefficient( + hamiltonian_with_coefficient, coefficient +): + with pytest.raises(ValueError, match="Coefficient must be"): + hamiltonian_with_coefficient(PauliString.from_str("X"), coefficient) + + +class _BrokenComplexCoefficient: + def __complex__(self): + raise RuntimeError("complex conversion failed") + + def __float__(self): + return 0.0 # A failed complex conversion must not fall back to this. + + +def test_hamiltonian_propagates_complex_conversion_error(hamiltonian_with_coefficient): + with pytest.raises(RuntimeError, match="complex conversion failed"): + hamiltonian_with_coefficient( + PauliString.from_str("X"), _BrokenComplexCoefficient() + ) + + +@pytest.mark.parametrize("method", ["add_term", "scale"]) +@pytest.mark.parametrize( + "coefficient, error", + [ + (object(), ValueError), + (("invalid", 1), TypeError), + (_BrokenComplexCoefficient(), RuntimeError), + ], +) +def test_hamiltonian_coefficient_error_preserves_state(method, coefficient, error): + h = Hamiltonian.from_list([(PauliString.from_str("Z"), 2)]) + before = h.copy() + with pytest.raises(error): + if method == "add_term": + h.add_term(PauliString.from_str("X"), coefficient) + else: + h.scale(coefficient) + assert h == before + np.testing.assert_array_equal(h.to_matrix(), [[2, 0], [0, -2]]) + + def test_checked_sparse_algebra_and_legacy_compatibility(): x = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1)]) z = Hamiltonian.from_pauli(PauliString.from_str("Z")) diff --git a/crates/binding-python/tests/qis/test_stabilizer.py b/crates/binding-python/tests/qis/test_stabilizer.py index 8114759d..7aa02ba5 100644 --- a/crates/binding-python/tests/qis/test_stabilizer.py +++ b/crates/binding-python/tests/qis/test_stabilizer.py @@ -272,3 +272,19 @@ def test_two_qubit_gates_reject_duplicate_qubits(self): with pytest.raises(ValueError): state.apply_swap(1, 1) + + +@pytest.mark.parametrize("theta", [0.0, 0.37, math.pi / 2]) +@pytest.mark.parametrize("entrypoint", ["from_circuit", "apply_circuit", "run_circuit"]) +def test_controlled_gphase_returns_normal_unsupported_error(theta, entrypoint): + from cqlib.circuit.gates import GPhase + + circuit = Circuit([12]) + circuit.append_mc_gate(GPhase(theta).control(1), [12]) + execute = ( + StabilizerState(1).apply_circuit + if entrypoint == "apply_circuit" + else getattr(StabilizerState, entrypoint) + ) + with pytest.raises(ValueError, match="Controlled-GPhase.*not supported"): + execute(circuit) diff --git a/crates/binding-python/tests/visualization/test_visualization.py b/crates/binding-python/tests/visualization/test_visualization.py index 42b1098d..331cf9d5 100644 --- a/crates/binding-python/tests/visualization/test_visualization.py +++ b/crates/binding-python/tests/visualization/test_visualization.py @@ -23,6 +23,7 @@ import subprocess import tempfile import zlib +from xml.etree import ElementTree import pytest import resvg_py @@ -641,6 +642,34 @@ def test_state_visualization_accepts_density_matrix(): assert vis.plot_state_paulivec(sv) == vis.plot_state_paulivec(dm) +@pytest.mark.parametrize("state_type", [Statevector, DensityMatrix]) +@pytest.mark.parametrize("reverse_bits", [False, True]) +def test_paulivec_labels_identify_the_correct_qubit(state_type, reverse_bits): + """The displayed IX bar measures q0 unless display bits are reversed.""" + state = state_type(2) + state.apply_h(0) + svg = vis.plot_state_paulivec(state, reverse_bits=reverse_bits) + elements = list(ElementTree.fromstring(svg)) + bars = { + element.text: elements[index - 1].attrib + for index, element in enumerate(elements) + if element.tag.endswith("text") + and element.text in {"II", "IX", "XI", "IZ", "ZI"} + } + one_y = float(bars["II"]["y"]) + one_height = float(bars["II"]["height"]) + assert one_height > 1.0 + for label, expected in [("IX", 1), ("XI", 0), ("IZ", 0), ("ZI", 1)]: + displayed = label[::-1] if reverse_bits else label + y = float(bars[displayed]["y"]) + height = float(bars[displayed]["height"]) + if expected: + assert (y, height) == (one_y, one_height) + else: + assert height == 1.0 + assert y == pytest.approx(one_y + one_height) + + def test_state_visualization_validates_python_inputs(): """Invalid Python inputs should fail with clear ValueError messages.""" with pytest.raises(ValueError, match="exactly 3"): diff --git a/crates/cqlib-core/README.md b/crates/cqlib-core/README.md index 696a7225..5549ab8e 100644 --- a/crates/cqlib-core/README.md +++ b/crates/cqlib-core/README.md @@ -14,13 +14,13 @@ Add the crate to your project: ```toml [dependencies] -cqlib-core = "0.1.0-beta.5" +cqlib-core = "0.1.0-beta.6" ``` Or use Cargo: ```shell -cargo add cqlib-core@0.1.0-beta.5 +cargo add cqlib-core@0.1.0-beta.6 ``` ## Quick start diff --git a/crates/cqlib-core/src/circuit/circuit_impl.rs b/crates/cqlib-core/src/circuit/circuit_impl.rs index 3cfffaf6..d670328f 100644 --- a/crates/cqlib-core/src/circuit/circuit_impl.rs +++ b/crates/cqlib-core/src/circuit/circuit_impl.rs @@ -2255,6 +2255,8 @@ impl Circuit { /// the result is `Rx(beta+1)`. /// 2. **Qubit Mapping**: Virtual qubits in the sub-circuit definition are mapped to the /// physical qubits in the parent circuit. + /// 3. **Global Phase**: Each sub-circuit's phase is substituted using its call arguments + /// and accumulated into the flattened circuit, including empty definitions. /// /// # Returns /// @@ -2263,8 +2265,9 @@ impl Circuit { /// /// # Errors /// - /// Returns [`CircuitError::UnresolvedParameter`] if a symbolic parameter in a sub-circuit - /// or control flow body cannot be evaluated to a concrete value. + /// Returns an error if a fully resolved parameter or accumulated global phase is invalid, + /// or if the circuit contains unsupported classical control. Unbound symbolic parameters + /// and phases are preserved. pub fn decompose(&self) -> Result { let mut new_circuit = Circuit::from_qubits(self.qubits()).unwrap(); new_circuit.classical_vars = self.classical_vars.clone(); @@ -2367,6 +2370,21 @@ impl Circuit { } } + // Phase is circuit metadata, so it must be accumulated even when + // the definition has no operations. Each gate application adds it once. + if cg.circuit.circuit.global_phase_param() != &CircuitParam::Fixed(0.0) { + let phase = + Self::apply_param_map(cg.circuit.circuit.global_phase(), &next_param_map); + if phase.get_symbols().is_empty() { + phase.evaluate(&None)?; + } + let phase = target_circuit.global_phase() + phase; + if phase.get_symbols().is_empty() { + phase.evaluate(&None)?; + } + target_circuit.set_global_phase(phase); + } + // Qubit Map: Inner Qubit -> Outer Qubit // op.qubits are the qubits in 'context_circuit' that the gate acts on. // We need to map them through 'qubit_map' to get 'target_circuit' qubits. diff --git a/crates/cqlib-core/src/circuit/circuit_test.rs b/crates/cqlib-core/src/circuit/circuit_test.rs index d8f476c8..373ef8dc 100644 --- a/crates/cqlib-core/src/circuit/circuit_test.rs +++ b/crates/cqlib-core/src/circuit/circuit_test.rs @@ -1447,6 +1447,259 @@ fn test_decompose() { } } +#[test] +fn test_decompose_global_phase_constant_and_empty_definitions() { + use crate::circuit::circuit_to_matrix; + use crate::circuit::test_utils::assert_matrix_approx_eq; + + for (width, has_operations) in [(0, false), (1, false), (1, true), (2, true)] { + let mut inner = Circuit::new(width); + inner.set_global_phase(Parameter::from(0.37)); + if has_operations { + inner.x(Qubit::new(0)).unwrap(); + if width == 2 { + inner.cx(Qubit::new(0), Qubit::new(1)).unwrap(); + } + } + let gate = inner.to_gate("Phased").unwrap(); + for repetitions in 0..=2 { + let mut outer = Circuit::new(2); + outer.set_global_phase(Parameter::from(0.12)); + let mut expected = Circuit::new(2); + expected.set_global_phase(Parameter::from(0.12 + repetitions as f64 * 0.37)); + for _ in 0..repetitions { + let targets = [Qubit::new(1), Qubit::new(0)]; + outer + .append(gate.clone(), targets[..width].to_vec(), [], None) + .unwrap(); + if has_operations { + expected.x(Qubit::new(1)).unwrap(); + if width == 2 { + expected.cx(Qubit::new(1), Qubit::new(0)).unwrap(); + } + } + } + let flat = outer.decompose().unwrap(); + assert!( + (flat.global_phase().evaluate(&None).unwrap() - (0.12 + repetitions as f64 * 0.37)) + .abs() + < 1e-12 + ); + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq( + &circuit_to_matrix(&outer, None).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&flat, None).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&flat.decompose().unwrap(), None).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_eq!(flat.operations().len(), expected.operations().len()); + assert_eq!(outer.global_phase().evaluate(&None).unwrap(), 0.12); + } + } +} + +#[test] +fn test_decompose_global_phase_only_symbols_and_explicit_signature() { + use crate::circuit::circuit_to_matrix; + use crate::circuit::test_utils::assert_matrix_approx_eq; + use std::collections::HashMap; + + let mut inner = Circuit::new(1); + inner.set_global_phase(Parameter::symbol("theta") / 2.0); + let gate = CircuitGate::with_signature( + "PhaseOnly", + FrozenCircuit::new(inner), + ["unused".to_string(), "theta".to_string()], + ) + .unwrap(); + let mut outer = Circuit::new(1); + outer.set_global_phase(Parameter::symbol("phi") / 3.0); + outer + .circuit_gate( + gate, + vec![Qubit::new(0)], + [ + ParameterValue::Fixed(0.7), + ParameterValue::from(Parameter::symbol("x") * 2.0), + ], + ) + .unwrap(); + let flat = outer.decompose().unwrap(); + assert!(flat.operations().is_empty()); + assert_eq!( + flat.used_symbols().iter().cloned().collect::>(), + HashSet::from(["phi".to_string(), "x".to_string()]) + ); + for (x, phi) in [(0.37, -0.23), (-0.61, 0.44)] { + let bindings = Some(HashMap::from([("x", x), ("phi", phi)])); + let mut expected = Circuit::new(1); + expected.set_global_phase(Parameter::from(phi / 3.0 + x)); + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq( + &circuit_to_matrix(&flat.assign_parameters(&bindings).unwrap(), None).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix( + &outer + .assign_parameters(&bindings) + .unwrap() + .decompose() + .unwrap(), + None, + ) + .unwrap(), + &expected_matrix, + 1e-10, + ); + } +} + +#[test] +fn test_decompose_global_phase_nested_symbolic_and_inverse() { + use crate::circuit::circuit_to_matrix; + use crate::circuit::symbolic_matrix::{circuit_to_symbolic_matrix, evaluate_symbolic_matrix}; + use crate::circuit::test_utils::assert_matrix_approx_eq; + use std::collections::HashMap; + + let a = Parameter::symbol("a"); + let b = Parameter::symbol("b"); + let mut inner = Circuit::from_qubits(vec![Qubit::new(7), Qubit::new(2)]).unwrap(); + inner.set_global_phase(a.clone() + b.clone() * 2.0); + inner.rx(Qubit::new(7), a.clone() - b.clone()).unwrap(); + inner.cx(Qubit::new(7), Qubit::new(2)).unwrap(); + let inner_gate = CircuitGate::with_signature( + "Inner", + FrozenCircuit::new(inner), + ["a".to_string(), "b".to_string()], + ) + .unwrap(); + let mut middle = Circuit::from_qubits(vec![Qubit::new(5), Qubit::new(1)]).unwrap(); + middle.set_global_phase(a.clone() / 2.0); + middle + .circuit_gate( + inner_gate, + vec![Qubit::new(5), Qubit::new(1)], + [ParameterValue::from(b), ParameterValue::from(a)], + ) + .unwrap(); + let middle_gate = CircuitGate::with_signature( + "Middle", + FrozenCircuit::new(middle), + ["a".to_string(), "b".to_string()], + ) + .unwrap(); + let x = Parameter::symbol("x"); + let mut outer = Circuit::new(2); + outer.set_global_phase(x.clone() / 3.0); + outer + .circuit_gate( + middle_gate, + vec![Qubit::new(1), Qubit::new(0)], + [ + ParameterValue::from(x + 1.0), + ParameterValue::from(Parameter::symbol("y")), + ], + ) + .unwrap(); + + for inverse in [false, true] { + let circuit = if inverse { + outer.inverse().unwrap() + } else { + outer.clone() + }; + let flat = circuit.decompose().unwrap(); + assert_eq!( + flat.used_symbols().iter().cloned().collect::>(), + HashSet::from(["x".to_string(), "y".to_string()]) + ); + let symbolic = circuit_to_symbolic_matrix(&flat, None).unwrap(); + for (x, y) in [(0.37, -0.23), (-0.61, 0.44)] { + let bindings = Some(HashMap::from([("x", x), ("y", y)])); + let mut expected = Circuit::new(2); + let expected_phase = x / 3.0 + (x + 1.0) / 2.0 + y + 2.0 * (x + 1.0); + expected.set_global_phase(Parameter::from(expected_phase)); + expected.rx(Qubit::new(1), y - (x + 1.0)).unwrap(); + expected.cx(Qubit::new(1), Qubit::new(0)).unwrap(); + if inverse { + expected = expected.inverse().unwrap(); + } + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + assert!( + (flat.global_phase().evaluate(&bindings).unwrap() + - if inverse { + -expected_phase + } else { + expected_phase + }) + .abs() + < 1e-12 + ); + for candidate in [ + flat.assign_parameters(&bindings).unwrap(), + circuit + .assign_parameters(&bindings) + .unwrap() + .decompose() + .unwrap(), + circuit.assign_parameters(&bindings).unwrap(), + ] { + assert_matrix_approx_eq( + &circuit_to_matrix(&candidate, None).unwrap(), + &expected_matrix, + 1e-10, + ); + } + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&symbolic, &bindings).unwrap(), + &expected_matrix, + 1e-10, + ); + } + } +} + +#[test] +fn test_decompose_global_phase_rejects_invalid_constants() { + let mut inner = Circuit::new(1); + inner.set_global_phase(Parameter::symbol("theta").ln()); + let mut outer = Circuit::new(1); + outer + .append( + inner.to_gate("LogPhase").unwrap(), + [Qubit::new(0)], + [ParameterValue::Fixed(-1.0)], + None, + ) + .unwrap(); + assert!(outer.decompose().is_err()); + + let mut inner = Circuit::new(1); + inner.set_global_phase(Parameter::from(1e308)); + let mut outer = Circuit::new(1); + outer.set_global_phase(Parameter::from(1e308)); + outer + .append( + inner.to_gate("LargePhase").unwrap(), + [Qubit::new(0)], + [], + None, + ) + .unwrap(); + assert!(outer.decompose().is_err()); +} + #[test] fn test_decompose_nested() { let mut l1 = Circuit::new(1); diff --git a/crates/cqlib-core/src/circuit/circuit_to_matrix.rs b/crates/cqlib-core/src/circuit/circuit_to_matrix.rs index 9dc5a5d4..f72835fb 100644 --- a/crates/cqlib-core/src/circuit/circuit_to_matrix.rs +++ b/crates/cqlib-core/src/circuit/circuit_to_matrix.rs @@ -198,23 +198,11 @@ pub fn circuit_to_matrix( } } Instruction::CircuitGate(circuit_gate) => { - let symbols = circuit_gate.signature_params(); - let expected = symbols.len(); - let actual = params.len(); - if actual != expected { - return Err(CircuitError::ParameterCountMismatch { expected, actual }); - } - let mut bindings = HashMap::new(); - for (sym, val) in symbols.iter().zip(params.iter()) { - bindings.insert(sym.as_str(), *val); - } - let sub_circuit = circuit_gate - .circuit - .circuit - .assign_parameters(&Some(bindings)) - .map_err(|_| CircuitError::SymbolicParameterError)?; - let sub_matrix = circuit_to_matrix(&sub_circuit, None)?; - apply_gate_to_matrix(&mut matrix, &sub_matrix, &bits)?; + let gate_matrix = circuit_gate.matrix_for_params(params.as_slice())?; + // Gate-local matrices put the first positional argument at MSB; + // the full circuit uses an LSB-first basis, as for standard gates. + let reversed_bits: SmallVec<[usize; 3]> = bits.iter().copied().rev().collect(); + apply_gate_to_matrix(&mut matrix, &gate_matrix, &reversed_bits)?; } Instruction::ClassicalControl(_) => { return Err(CircuitError::InvalidOperation( diff --git a/crates/cqlib-core/src/circuit/circuit_to_matrix_test.rs b/crates/cqlib-core/src/circuit/circuit_to_matrix_test.rs index 86754c4d..5f223d51 100644 --- a/crates/cqlib-core/src/circuit/circuit_to_matrix_test.rs +++ b/crates/cqlib-core/src/circuit/circuit_to_matrix_test.rs @@ -419,6 +419,76 @@ fn test_unitary_gate_two_qubit() { assert!((matrix[[1, 2]] - c(1.0, 0.0)).norm() < 1e-10); } +#[test] +fn test_circuit_gate_storage_order_mapping() { + for order in [[7, 2], [1, 0], [2, 7], [0, 1]] { + let mut inner = Circuit::from_qubits(order.map(Qubit::new).to_vec()).unwrap(); + inner.x(Qubit::new(order[0])).unwrap(); + let gate = inner.to_gate("StoredOrderX").unwrap(); + for targets in [[0, 1], [1, 0]] { + let mut circuit = Circuit::new(2); + circuit + .append(gate.clone(), targets.map(Qubit::new), [], None) + .unwrap(); + let mut expected = Circuit::new(2); + expected.x(Qubit::new(targets[0])).unwrap(); + let matrix = circuit_to_matrix(&circuit, None).unwrap(); + let expected = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq(&matrix, &expected, 1e-10); + assert_eq!(matrix[[1 << targets[0], 0]], c(1.0, 0.0)); + assert_matrix_approx_eq( + &matrix, + &circuit_to_matrix(&circuit.decompose().unwrap(), None).unwrap(), + 1e-10, + ); + } + } +} + +#[test] +fn test_circuit_gate_storage_order_three_qubit_permutations() { + let permutations = [ + [0, 1, 2], + [0, 2, 1], + [1, 0, 2], + [1, 2, 0], + [2, 0, 1], + [2, 1, 0], + ]; + let inner_ids = [2, 7, 11]; + let outer_ids = [4, 9, 15]; + for permutation in permutations { + let order = permutation.map(|i| Qubit::new(inner_ids[i])); + let mut inner = Circuit::from_qubits(order.to_vec()).unwrap(); + inner.x(order[0]).unwrap(); + inner.ry(order[1], 0.37).unwrap(); + inner.cx(order[0], order[2]).unwrap(); + let gate = inner.to_gate("Permuted").unwrap(); + for target_permutation in permutations { + let targets = target_permutation.map(|i| Qubit::new(outer_ids[i])); + let mut circuit = Circuit::from_qubits(outer_ids.map(Qubit::new).to_vec()).unwrap(); + circuit.append(gate.clone(), targets, [], None).unwrap(); + let mut expected = Circuit::from_qubits(outer_ids.map(Qubit::new).to_vec()).unwrap(); + expected.x(targets[0]).unwrap(); + expected.ry(targets[1], 0.37).unwrap(); + expected.cx(targets[0], targets[2]).unwrap(); + for output_order in [None, Some([15, 4, 9].as_slice())] { + let matrix = circuit_to_matrix(&circuit, output_order).unwrap(); + assert_matrix_approx_eq( + &matrix, + &circuit_to_matrix(&expected, output_order).unwrap(), + 1e-10, + ); + assert_matrix_approx_eq( + &matrix, + &circuit_to_matrix(&circuit.decompose().unwrap(), output_order).unwrap(), + 1e-10, + ); + } + } + } +} + #[test] fn test_circuit_gate_bell_state() { // Create a Bell state preparation circuit as a gate @@ -866,6 +936,74 @@ fn test_global_phase_is_included_for_empty_circuit() { assert_matrix_approx_eq(&matrix, &expected, 1e-10); } +#[test] +fn test_gphase_markers_accumulate_with_circuit_global_phase() { + for num_qubits in [0, 1, 3] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut circuit = Circuit::new(num_qubits); + if num_qubits > 0 { + circuit.h(Qubit::new(0)).unwrap(); + } + let baseline = circuit_to_matrix(&circuit, None).unwrap(); + for phase in [theta, -0.21] { + circuit + .append( + Instruction::Standard(StandardGate::GPhase), + [] as [Qubit; 0], + [phase.into()], + None, + ) + .unwrap(); + } + circuit.set_global_phase(Parameter::from(0.17)); + let expected = baseline * Complex64::from_polar(1.0, theta - 0.21 + 0.17); + let actual = circuit_to_matrix(&circuit, None).unwrap(); + assert_matrix_approx_eq(&actual, &expected, 1e-12); + } + } +} + +#[test] +fn test_gphase_controlled_embedding_with_spectator_and_custom_order() { + for controls in [vec![], vec![2usize], vec![2, 0]] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut circuit = Circuit::new(3); + circuit.h(Qubit::new(1)).unwrap(); + circuit + .multi_control( + StandardGate::GPhase, + controls.iter().map(|&q| Qubit::new(q as u32)), + [] as [Qubit; 0], + [theta.into()], + ) + .unwrap(); + for order in [[0, 1, 2], [2, 0, 1]] { + let spectator_bit = order.iter().position(|&q| q == 1).unwrap(); + let mask = controls.iter().fold(0, |mask, control| { + mask | (1 << order.iter().position(|q| q == control).unwrap()) + }); + let mut expected = Array2::zeros((8, 8)); + for row in 0..8 { + for col in 0..8 { + if (row ^ col) & !(1 << spectator_bit) == 0 { + let sign = if row & col & (1 << spectator_bit) != 0 { + -1.0 + } else { + 1.0 + }; + let phase = if row & mask == mask { theta } else { 0.0 }; + expected[[row, col]] = + Complex64::from_polar(1.0, phase) * sign / SQRT_2; + } + } + } + let actual = circuit_to_matrix(&circuit, Some(&order)).unwrap(); + assert_matrix_approx_eq(&actual, &expected, 1e-12); + } + } + } +} + #[test] fn test_global_phase_is_included_with_gate() { let mut circuit = Circuit::new(1); diff --git a/crates/cqlib-core/src/circuit/gate/circuit_gate.rs b/crates/cqlib-core/src/circuit/gate/circuit_gate.rs index 7cbeca3d..5eaf87a9 100644 --- a/crates/cqlib-core/src/circuit/gate/circuit_gate.rs +++ b/crates/cqlib-core/src/circuit/gate/circuit_gate.rs @@ -16,10 +16,15 @@ //! quantum circuits to be used as reusable gate components within other circuits. //! This enables hierarchical circuit construction and custom composite gates. +use crate::circuit::Qubit; use crate::circuit::circuit_impl::Circuit; +use crate::circuit::circuit_to_matrix; use crate::circuit::error::CircuitError; use crate::circuit::symbolic_matrix::{SymbolicMatrix, circuit_to_symbolic_matrix}; use indexmap::IndexSet; +use ndarray::Array2; +use num_complex::Complex64; +use std::collections::HashMap; use std::sync::{Arc, OnceLock}; /// An immutable, frozen circuit for use in gate definitions. @@ -46,6 +51,7 @@ pub struct FrozenCircuit { pub(crate) circuit: Circuit, used_symbols: Arc>, symbolic_matrix_cache: Arc>>, + storage_order_symbolic_matrix_cache: Arc>>, } impl FrozenCircuit { @@ -67,6 +73,7 @@ impl FrozenCircuit { circuit, used_symbols, symbolic_matrix_cache: Arc::new(OnceLock::new()), + storage_order_symbolic_matrix_cache: Arc::new(OnceLock::new()), } } @@ -95,6 +102,16 @@ impl FrozenCircuit { &self.used_symbols } + /// The default matrix order, from least to most significant bit. + /// + /// Circuit-backed UnitaryGate definitions map this order positionally to + /// the caller's qubits. CircuitGate instead uses circuit storage order. + pub(crate) fn qubits_in_numeric_order(&self) -> Vec { + let mut qubits = self.circuit.qubits(); + qubits.sort_unstable(); + qubits + } + /// Returns the cached symbolic matrix for this frozen circuit. /// /// The matrix is computed in the circuit's default qubit order on the @@ -105,10 +122,32 @@ impl FrozenCircuit { return Ok(matrix.clone()); } - let matrix = Arc::new(circuit_to_symbolic_matrix(&self.circuit, None)?); + let order: Vec = self + .qubits_in_numeric_order() + .iter() + .map(|q| q.index()) + .collect(); + let matrix = Arc::new(circuit_to_symbolic_matrix(&self.circuit, Some(&order))?); let _ = self.symbolic_matrix_cache.set(matrix.clone()); Ok(self.symbolic_matrix_cache.get().cloned().unwrap_or(matrix)) } + + /// CircuitGate's positional basis is distinct from the numeric-sorted + /// default used by the public FrozenCircuit matrix and UnitaryGate. + fn symbolic_matrix_in_storage_order(&self) -> Result, CircuitError> { + if let Some(matrix) = self.storage_order_symbolic_matrix_cache.get() { + return Ok(matrix.clone()); + } + + let order: Vec = self.circuit.qubits().iter().map(|q| q.index()).collect(); + let matrix = Arc::new(circuit_to_symbolic_matrix(&self.circuit, Some(&order))?); + let _ = self.storage_order_symbolic_matrix_cache.set(matrix.clone()); + Ok(self + .storage_order_symbolic_matrix_cache + .get() + .cloned() + .unwrap_or(matrix)) + } } impl PartialEq for FrozenCircuit { @@ -331,6 +370,40 @@ impl CircuitGate { self.num_params } + /// Resolves a numeric gate-local matrix, with the first positional + /// argument as the most-significant bit, like a standard gate matrix. + /// This differs from `symbolic_matrix`'s storage-order, little-endian basis. + pub(crate) fn matrix_for_params( + &self, + params: &[f64], + ) -> Result, CircuitError> { + if params.len() != self.signature_params.len() { + return Err(CircuitError::ParameterCountMismatch { + expected: self.signature_params.len(), + actual: params.len(), + }); + } + for (index, &value) in params.iter().enumerate() { + if !value.is_finite() { + return Err(CircuitError::InvalidParameterValue(index, value)); + } + } + let bindings: HashMap<&str, f64> = self + .signature_params + .iter() + .map(String::as_str) + .zip(params.iter().copied()) + .collect(); + let assigned = self + .circuit + .circuit + .assign_parameters(&Some(bindings)) + .map_err(|_| CircuitError::SymbolicParameterError)?; + // circuit_to_matrix takes LSB-first ordering; gate arguments are MSB-first. + let order: Vec = assigned.qubits().iter().rev().map(|q| q.index()).collect(); + circuit_to_matrix(&assigned, Some(&order)) + } + /// Returns a clone of the internal frozen circuit. /// /// # Examples @@ -350,9 +423,12 @@ impl CircuitGate { self.circuit.clone() } - /// Returns the cached symbolic matrix for this gate's frozen circuit. + /// Returns the cached symbolic matrix in this gate's positional basis. + /// + /// The first qubit in the inner circuit's storage order is the + /// least-significant matrix bit, matching the first application argument. pub fn symbolic_matrix(&self) -> Result, CircuitError> { - self.circuit.symbolic_matrix() + self.circuit.symbolic_matrix_in_storage_order() } /// Returns the name of this circuit gate. @@ -412,6 +488,82 @@ mod tests { use super::*; use crate::circuit::{Circuit, Parameter, Qubit}; + #[test] + fn numeric_matrix_binds_declared_signature_and_global_phase_without_mutation() { + let mut circuit = Circuit::from_qubits(vec![Qubit::new(7)]).unwrap(); + let theta = Parameter::symbol("theta"); + let phi = Parameter::symbol("phi"); + circuit.rx(Qubit::new(7), theta.clone()).unwrap(); + circuit.rz(Qubit::new(7), phi.clone()).unwrap(); + circuit.set_global_phase(theta + phi / Parameter::from(2.0)); + let gate = CircuitGate::with_signature( + "declared", + FrozenCircuit::new(circuit), + ["unused".to_string(), "phi".to_string(), "theta".to_string()], + ) + .unwrap(); + let original = gate.clone(); + for (theta, phi) in [(0.37_f64, -0.21_f64), (-0.61, 0.42), (0.37, -0.21)] { + let matrix = gate.matrix_for_params(&[12.0, phi, theta]).unwrap(); + for row in 0..2 { + let phase = theta + phi / 2.0 + if row == 0 { -phi / 2.0 } else { phi / 2.0 }; + for col in 0..2 { + let rx = if row == col { + Complex64::new((theta / 2.0).cos(), 0.0) + } else { + Complex64::new(0.0, -(theta / 2.0).sin()) + }; + assert!( + (matrix[[row, col]] - Complex64::from_polar(1.0, phase) * rx).norm() + < 1e-12 + ); + } + } + assert_eq!(gate, original); + assert_eq!(gate.used_symbols().len(), 2); + } + } + + #[test] + fn numeric_matrix_rejects_invalid_arguments_and_nonunitary_definitions() { + let gate = CircuitGate::with_signature( + "unused", + FrozenCircuit::new(Circuit::new(1)), + ["unused".to_string()], + ) + .unwrap(); + for params in [vec![], vec![0.1, 0.2]] { + assert!(matches!( + gate.matrix_for_params(¶ms), + Err(CircuitError::ParameterCountMismatch { expected: 1, .. }) + )); + } + for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert!(matches!( + gate.matrix_for_params(&[value]), + Err(CircuitError::InvalidParameterValue(0, _)) + )); + } + for directive in [ + crate::circuit::Directive::Measure, + crate::circuit::Directive::Reset, + ] { + let mut body = Circuit::new(1); + body.append( + crate::circuit::Instruction::Directive(directive), + [Qubit::new(0)], + [], + None, + ) + .unwrap(); + let gate = CircuitGate::new("nonunitary", FrozenCircuit::new(body)).unwrap(); + assert!(matches!( + gate.matrix_for_params(&[]), + Err(CircuitError::NoMatrixRepresentation) + )); + } + } + #[test] fn signature_params_are_distinct_from_used_symbols() { let mut circuit = Circuit::new(1); diff --git a/crates/cqlib-core/src/circuit/gate/gate_matrix.rs b/crates/cqlib-core/src/circuit/gate/gate_matrix.rs index f4792fc8..5477911f 100644 --- a/crates/cqlib-core/src/circuit/gate/gate_matrix.rs +++ b/crates/cqlib-core/src/circuit/gate/gate_matrix.rs @@ -581,17 +581,14 @@ pub fn phase_gate(lambda: f64) -> Array2> { array![[ONE, ZERO], [ZERO, exp_i_lambda]] } -/// Returns a global phase gate. +/// Returns the 1×1 matrix of a zero-qubit global phase gate. /// -/// This applies a global phase factor to the quantum state, physically equivalent to identity -/// but useful for mathematical consistency in simulations. +/// When embedded in a circuit, this scalar multiplies every state amplitude +/// by the same phase factor without requiring a target qubit. /// /// Defined as: /// ```text -/// [ -/// [exp(i*theta), 0 ], -/// [ 0, exp(i*theta)] -/// ] +/// [[exp(i*theta)]] /// ``` /// /// # Arguments @@ -601,7 +598,7 @@ pub fn phase_gate(lambda: f64) -> Array2> { pub fn global_phase_gate(theta: f64) -> Array2> { let (sin_val, cos_val) = theta.sin_cos(); let exp_i_theta = Complex::new(cos_val, sin_val); - array![[exp_i_theta, ZERO], [ZERO, exp_i_theta]] + array![[exp_i_theta]] } /// Returns the general single-qubit unitary gate (U3). diff --git a/crates/cqlib-core/src/circuit/gate/gate_matrix_test.rs b/crates/cqlib-core/src/circuit/gate/gate_matrix_test.rs index a6958d55..9444296e 100644 --- a/crates/cqlib-core/src/circuit/gate/gate_matrix_test.rs +++ b/crates/cqlib-core/src/circuit/gate/gate_matrix_test.rs @@ -11,6 +11,8 @@ // that they have been altered from the originals. use super::*; +use crate::circuit::Parameter; +use crate::circuit::gate::{MCGate, StandardGate}; use crate::circuit::test_utils::assert_is_unitary; use std::f64::consts::PI; @@ -27,6 +29,45 @@ fn assert_matrix_approx_eq(a: &Array2>, b: &Array2>) { } } +#[test] +fn test_gphase_scalar_matrix() { + for theta in [0.0, 0.37, -0.51, PI] { + let matrix = global_phase_gate(theta); + assert_eq!(matrix.dim(), (1, 1)); + assert!((matrix[[0, 0]] - Complex::from_polar(1.0, theta)).norm() < 1e-12); + assert_is_unitary(&matrix, 1e-12); + } +} + +#[test] +fn test_gphase_controlled_matrix_matches_metadata_and_inverse() { + assert_eq!(StandardGate::GPhase.num_qubits(), 0); + assert_eq!(StandardGate::GPhase.num_ctrl_qubits(), 0); + assert_eq!(StandardGate::GPhase.num_params(), 1); + for theta in [0.0, 0.37, -0.51, PI] { + for controls in 0..=2 { + let gate = MCGate::new(controls, StandardGate::GPhase); + assert_eq!(gate.num_qubits(), controls as usize); + let dim = 1 << gate.num_qubits(); + let matrix = gate.matrix(&[theta]).unwrap(); + let mut expected = Array2::eye(dim); + expected[[dim - 1, dim - 1]] = Complex::from_polar(1.0, theta); + crate::circuit::test_utils::assert_matrix_approx_eq(&matrix, &expected, 1e-12); + assert_is_unitary(&matrix, 1e-12); + + let (inverse, params) = gate.inverse(&[Parameter::from(theta)]).unwrap(); + let inverse_theta = params[0].evaluate(&None).unwrap(); + let inverse_matrix = inverse.matrix(&[inverse_theta]).unwrap(); + let expected_inverse = expected.t().mapv(|value| value.conj()); + crate::circuit::test_utils::assert_matrix_approx_eq( + &inverse_matrix, + &expected_inverse, + 1e-12, + ); + } + } +} + #[test] fn test_fixed_gates_unitary() { // Batch verify unitarity for all fixed gates diff --git a/crates/cqlib-core/src/circuit/gate/instruction.rs b/crates/cqlib-core/src/circuit/gate/instruction.rs index 31b9513c..7047f713 100644 --- a/crates/cqlib-core/src/circuit/gate/instruction.rs +++ b/crates/cqlib-core/src/circuit/gate/instruction.rs @@ -26,7 +26,7 @@ use crate::circuit::gate::circuit_gate::{CircuitGate, FrozenCircuit}; use crate::circuit::gate::directive::Directive; use crate::circuit::gate::standard_gate::StandardGate; use crate::circuit::gate::{ClassicalDataOp, MCGate, UnitaryGate, gate_matrix}; -use crate::circuit::{ClassicalControlOp, ClassicalValue, Parameter, circuit_to_matrix}; +use crate::circuit::{ClassicalControlOp, ClassicalValue, Parameter}; use alloc::borrow::Cow; use ndarray::Array2; use num_complex::Complex64; @@ -318,23 +318,26 @@ impl Instruction { /// Computes the unitary matrix representation of the instruction. /// + /// The first qubit argument is the most-significant bit of the local + /// matrix for every supported gate, including circuit-backed UnitaryGate + /// definitions. Full circuit matrices instead use an LSB-first basis. + /// /// # Arguments /// /// * `params` - A slice of floating-point parameters (e.g., rotation angles). + /// CircuitGate values bind positionally to its declared signature. /// /// # Returns /// - /// - `Some(Cow)`: The unitary matrix for `Standard` and `Extended` gates. - /// - `None`: If the instruction is non-unitary (e.g., `Measure`, `Barrier`, `Reset`, - /// or classical control). + /// - `Some(Cow)`: A numeric matrix for a supported gate instruction. + /// - `None`: If parameters are invalid or unresolved, or the instruction has + /// no matrix (e.g., `Measure`, `Barrier`, `Reset`, or classical control). pub fn matrix(&self, params: &[f64]) -> Option>> { match self { Instruction::Standard(g) => g.matrix(params).ok(), Instruction::McGate(g) => g.matrix(params).ok(), - Instruction::UnitaryGate(g) => g.matrix_for_params(params).ok(), - Instruction::CircuitGate(g) => circuit_to_matrix(&g.circuit.circuit, None) - .ok() - .map(Cow::Owned), + Instruction::UnitaryGate(g) => g.matrix_for_params_in_gate_order(params).ok(), + Instruction::CircuitGate(g) => g.matrix_for_params(params).ok().map(Cow::Owned), Instruction::Directive(_) => None, Instruction::ClassicalData(_) => None, Instruction::Delay => None, diff --git a/crates/cqlib-core/src/circuit/gate/instruction_test.rs b/crates/cqlib-core/src/circuit/gate/instruction_test.rs index 33f33779..426cf7b1 100644 --- a/crates/cqlib-core/src/circuit/gate/instruction_test.rs +++ b/crates/cqlib-core/src/circuit/gate/instruction_test.rs @@ -13,9 +13,179 @@ use super::*; use crate::circuit::{ - Circuit, CircuitId, ClassicalExpr, ClassicalType, ClassicalValue, ClassicalVar, + Circuit, CircuitId, ClassicalExpr, ClassicalType, ClassicalValue, ClassicalVar, Qubit, }; +#[test] +fn circuit_backed_unitary_instruction_matrix_uses_msb_first_qubits() { + use crate::circuit::circuit_to_matrix; + use crate::circuit::test_utils::assert_matrix_approx_eq; + + for order in [vec![], vec![7], vec![0, 1], vec![1, 0], vec![9, 2, 7]] { + let mut numeric_order = order.clone(); + numeric_order.sort_unstable(); + for parameterized in [false, true] { + let mut body = + Circuit::from_qubits(order.iter().copied().map(Qubit::new).collect()).unwrap(); + if let Some(&first) = numeric_order.first() { + body.rx( + Qubit::new(first), + if parameterized { + Parameter::symbol("theta") + } else { + Parameter::from(0.7) + }, + ) + .unwrap(); + if numeric_order.len() > 1 { + body.cx( + Qubit::new(first), + Qubit::new(*numeric_order.last().unwrap()), + ) + .unwrap(); + } + } + body.set_global_phase(if parameterized { + Parameter::symbol("theta") + } else { + Parameter::from(0.7) + }); + let gate = UnitaryGate::new("positional", order.len() as u16, u16::from(parameterized)) + .with_circuit(Arc::new(FrozenCircuit::new(body.clone()))) + .unwrap(); + let params = if parameterized { vec![0.7] } else { vec![] }; + let bindings = Some(std::collections::HashMap::from([("theta", 0.7)])); + let bound = body.assign_parameters(&bindings).unwrap(); + let lsb_order: Vec = numeric_order.iter().map(|&q| q as usize).collect(); + // The definition API keeps its existing numeric-ID, LSB-first basis. + assert_matrix_approx_eq( + gate.matrix_for_params(¶ms).unwrap().as_ref(), + &circuit_to_matrix(&bound, Some(&lsb_order)).unwrap(), + 1e-12, + ); + + let instruction = Instruction::UnitaryGate(Box::new(gate)); + let matrix = instruction.matrix(¶ms).unwrap(); + let qubits: Vec = [11, 3, 8][..order.len()] + .iter() + .copied() + .map(Qubit::new) + .collect(); + let mut expected = Circuit::from_qubits(qubits.clone()).unwrap(); + if let Some(&first) = qubits.first() { + expected.rx(first, Parameter::from(0.7)).unwrap(); + if qubits.len() > 1 { + expected.cx(first, *qubits.last().unwrap()).unwrap(); + } + } + expected.set_global_phase(Parameter::from(0.7)); + let msb_order: Vec = qubits.iter().rev().map(|q| q.index()).collect(); + assert_matrix_approx_eq( + matrix.as_ref(), + &circuit_to_matrix(&expected, Some(&msb_order)).unwrap(), + 1e-12, + ); + let operation = crate::circuit::Operation { + instruction: instruction.clone(), + qubits: qubits.into_iter().collect(), + params: params + .into_iter() + .map(crate::circuit::CircuitParam::Fixed) + .collect(), + label: None, + }; + assert_matrix_approx_eq(operation.matrix().unwrap().as_ref(), matrix.as_ref(), 1e-12); + if parameterized { + for invalid in [vec![], vec![f64::NAN], vec![f64::INFINITY], vec![0.7, 0.2]] { + assert!(instruction.matrix(&invalid).is_none()); + } + } + } + } +} + +#[test] +fn unitary_instruction_matrix_keeps_explicit_matrix_order_and_precedence() { + use crate::circuit::symbolic_matrix::standard_gate_symbolic_matrix; + + let mut body = Circuit::from_qubits(vec![Qubit::new(7), Qubit::new(2)]).unwrap(); + body.cx(Qubit::new(2), Qubit::new(7)).unwrap(); + let frozen = Arc::new(FrozenCircuit::new(body)); + let numeric = StandardGate::CX.matrix(&[]).unwrap(); + for with_circuit in [false, true] { + for symbolic in [false, true] { + let gate = UnitaryGate::new("explicit_cx", 2, 0); + let gate = if symbolic { + gate.with_symbolic_matrix( + std::iter::empty::(), + standard_gate_symbolic_matrix(StandardGate::CX, &[]).unwrap(), + ) + } else { + gate.with_matrix(numeric.clone().into_owned()) + } + .unwrap(); + let gate = if with_circuit { + gate.with_circuit(frozen.clone()).unwrap() + } else { + gate + }; + let instruction = Instruction::UnitaryGate(Box::new(gate)); + let matrix = instruction.matrix(&[]).unwrap(); + assert_eq!(matrix.as_ref(), numeric.as_ref()); + assert_eq!(matches!(matrix, Cow::Borrowed(_)), !symbolic); + } + } +} + +#[test] +fn circuit_gate_operation_matrix_respects_positional_qubits() { + for order in [[0, 1], [1, 0], [7, 2]] { + let mut body = Circuit::from_qubits(order.map(Qubit::new).to_vec()).unwrap(); + body.x(Qubit::new(order[0])).unwrap(); + let instruction = body.to_gate("first_x").unwrap(); + let matrix = instruction.matrix(&[]).unwrap(); + // First argument is the MSB of a gate-local matrix: X tensor I. + for row in 0..4 { + for col in 0..4 { + let expected = if row == (col ^ 2) { 1.0 } else { 0.0 }; + assert_eq!(matrix[[row, col]], Complex64::new(expected, 0.0)); + } + } + } +} + +#[test] +fn circuit_gate_operation_matrix_binds_numeric_arguments() { + let mut body = Circuit::from_qubits(vec![Qubit::new(7)]).unwrap(); + body.rx(Qubit::new(7), Parameter::symbol("theta")).unwrap(); + let instruction = body.to_gate("rx").unwrap(); + for angle in [0.37_f64, -0.61] { + let matrix = instruction.matrix(&[angle]).unwrap(); + let mut operation = crate::circuit::Operation { + instruction: instruction.clone(), + qubits: smallvec::smallvec![Qubit::new(7)], + params: smallvec::smallvec![crate::circuit::CircuitParam::Fixed(angle)], + label: None, + }; + assert_eq!(operation.matrix().unwrap().as_ref(), matrix.as_ref()); + operation.params[0] = crate::circuit::CircuitParam::Index(0); + assert!(matches!( + operation.matrix(), + Err(crate::circuit::CircuitError::SymbolicParameterError) + )); + for row in 0..2 { + for col in 0..2 { + let expected = if row == col { + Complex64::new((angle / 2.0).cos(), 0.0) + } else { + Complex64::new(0.0, -(angle / 2.0).sin()) + }; + assert!((matrix[[row, col]] - expected).norm() < 1e-12); + } + } + } +} + #[test] fn standard_gate_query_distinguishes_instruction_kinds() { assert_eq!( diff --git a/crates/cqlib-core/src/circuit/gate/unitary_gate.rs b/crates/cqlib-core/src/circuit/gate/unitary_gate.rs index 9e574e15..35e104c6 100644 --- a/crates/cqlib-core/src/circuit/gate/unitary_gate.rs +++ b/crates/cqlib-core/src/circuit/gate/unitary_gate.rs @@ -286,9 +286,31 @@ impl UnitaryGate { } /// Resolves this gate into a concrete matrix for the supplied parameter values. + /// + /// Explicit matrices use MSB-first argument order. A circuit-only definition + /// retains its numeric-ID, LSB-first basis for compatibility with circuit + /// matrix construction. Use [`Instruction::matrix`](super::Instruction::matrix) + /// for a uniform MSB-first gate-local matrix. pub fn matrix_for_params( &self, params: &[f64], + ) -> Result>>, CircuitError> { + self.resolve_matrix_for_params(params, false) + } + + /// Resolves a gate-local matrix with the first argument at MSB, regardless + /// of whether the definition supplies a matrix or a circuit. + pub(crate) fn matrix_for_params_in_gate_order( + &self, + params: &[f64], + ) -> Result>>, CircuitError> { + self.resolve_matrix_for_params(params, true) + } + + fn resolve_matrix_for_params( + &self, + params: &[f64], + circuit_msb_first: bool, ) -> Result>>, CircuitError> { self.validate_param_values(params)?; @@ -323,7 +345,17 @@ impl UnitaryGate { let assigned = inner .assign_parameters(&Some(bindings)) .map_err(|_| CircuitError::SymbolicParameterError)?; - return circuit_to_matrix(&assigned, None).map(Cow::Owned); + let mut order: Vec = circuit + .qubits_in_numeric_order() + .iter() + .map(|q| q.index()) + .collect(); + if circuit_msb_first { + // circuit_to_matrix takes an LSB-first order. Generate the + // instruction basis directly instead of permuting a dense matrix. + order.reverse(); + } + return circuit_to_matrix(&assigned, Some(&order)).map(Cow::Owned); } Err(CircuitError::NoMatrixRepresentation) @@ -407,6 +439,12 @@ impl UnitaryGate { /// /// This allows the gate to be defined by its circuit decomposition, /// which is useful for inverse operations and optimization. + /// The internal qubits are mapped in ascending numeric-ID order to the + /// caller's qubits, independently of their storage order in the circuit. + /// When no explicit matrix is attached, the first such qubit is the least + /// significant bit of the derived matrix. Explicit matrices retain their + /// MSB-first argument convention; if both representations are supplied, + /// they must describe the same operation under their respective conventions. /// /// # Arguments /// diff --git a/crates/cqlib-core/src/circuit/gate/unitary_gate_test.rs b/crates/cqlib-core/src/circuit/gate/unitary_gate_test.rs index 87e83d6a..8c197ab2 100644 --- a/crates/cqlib-core/src/circuit/gate/unitary_gate_test.rs +++ b/crates/cqlib-core/src/circuit/gate/unitary_gate_test.rs @@ -16,7 +16,9 @@ use crate::circuit::gate::StandardGate; use crate::circuit::symbolic_matrix::{ SymbolicComplex, standard_gate_symbolic_matrix, symbolic_eye, }; +use crate::circuit::test_utils::assert_matrix_approx_eq; use crate::circuit::{Circuit, Parameter, Qubit}; +use crate::compile::transform::decompose::expand_definitions; use ndarray::array; use num_complex::Complex64; @@ -646,3 +648,252 @@ fn test_custom_unitary_in_composite_circuit() { } } } + +/// Check complete operators, including spectator qubits and nested CircuitGate +/// definitions, against an independently constructed positional reference. +fn assert_circuit_definition_order(transform: impl Fn(&Circuit) -> Circuit) { + let mut inner = Circuit::from_qubits(vec![Qubit::new(2), Qubit::new(0)]).unwrap(); + inner.x(Qubit::new(2)).unwrap(); + let gate = UnitaryGate::new("asymmetric_x", 2, 0) + .with_circuit(Arc::new(FrozenCircuit::new(inner))) + .unwrap(); + let mut circuit = Circuit::new(2); + circuit + .unitary(gate, vec![Qubit::new(0), Qubit::new(1)]) + .unwrap(); + let mut expected = Circuit::new(2); + expected.x(Qubit::new(1)).unwrap(); + let expected = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq( + &circuit_to_matrix(&circuit, None).unwrap(), + &expected, + 1e-12, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&transform(&circuit), None).unwrap(), + &expected, + 1e-12, + ); + + let mut leaf = Circuit::from_qubits(vec![Qubit::new(4), Qubit::new(1)]).unwrap(); + leaf.h(Qubit::new(4)).unwrap(); + leaf.cx(Qubit::new(4), Qubit::new(1)).unwrap(); + let leaf = leaf.to_gate("storage_order_leaf").unwrap(); + + for order in [[9, 2, 5], [5, 9, 2], [2, 5, 9]] { + let mut inner = Circuit::from_qubits(order.map(Qubit::new).to_vec()).unwrap(); + inner.rx(Qubit::new(9), 0.37).unwrap(); + inner.ry(Qubit::new(2), -0.51).unwrap(); + inner.cx(Qubit::new(5), Qubit::new(9)).unwrap(); + inner.cx(Qubit::new(9), Qubit::new(2)).unwrap(); + inner.rz(Qubit::new(5), 0.29).unwrap(); + inner + .append(leaf.clone(), [Qubit::new(9), Qubit::new(5)], [], None) + .unwrap(); + // Explicit matrices use MSB-first qargs. Reverse the definition's + // numeric order when also supplying an equivalent matrix representation. + let explicit_matrix = circuit_to_matrix(&inner, Some(&[9, 5, 2])).unwrap(); + let gate = UnitaryGate::new("numeric_order_unitary", 3, 0) + .with_circuit(Arc::new(FrozenCircuit::new(inner))) + .unwrap(); + for with_matrix in [false, true] { + let gate = if with_matrix { + gate.clone().with_matrix(explicit_matrix.clone()).unwrap() + } else { + gate.clone() + }; + for args in [ + [0, 2, 3], + [0, 3, 2], + [2, 0, 3], + [2, 3, 0], + [3, 0, 2], + [3, 2, 0], + ] { + let [a, b, c] = args.map(Qubit::new); + let mut circuit = Circuit::new(4); + circuit.unitary(gate.clone(), vec![a, b, c]).unwrap(); + let mut expected = Circuit::new(4); + expected.rx(c, 0.37).unwrap(); + expected.ry(a, -0.51).unwrap(); + expected.cx(b, c).unwrap(); + expected.cx(c, a).unwrap(); + expected.rz(b, 0.29).unwrap(); + expected.h(c).unwrap(); + expected.cx(c, b).unwrap(); + let expected = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq( + &circuit_to_matrix(&circuit, None).unwrap(), + &expected, + 1e-12, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&transform(&circuit), None).unwrap(), + &expected, + 1e-12, + ); + } + } + } +} + +#[test] +fn circuit_definition_order_survives_expansion() { + assert_circuit_definition_order(|circuit| expand_definitions(circuit).unwrap()); +} + +#[test] +fn circuit_definition_order_survives_qasm2() { + assert_circuit_definition_order(|circuit| { + crate::ir::qasm2::load::loads(&crate::ir::qasm2::dump::dumps(circuit).unwrap()).unwrap() + }); +} + +#[test] +fn circuit_definition_order_survives_qasm3() { + assert_circuit_definition_order(|circuit| { + crate::ir::qasm3::load::loads(&crate::ir::qasm3::dump::dumps(circuit).unwrap()).unwrap() + }); +} + +#[test] +fn circuit_definition_order_survives_compilation() { + use crate::compile::resource::ResourcePolicy; + use crate::compile::{CompileConfig, CompileMode, CompileTarget, compile}; + for mode in [CompileMode::Normal, CompileMode::Enhanced] { + assert_circuit_definition_order(|circuit| { + compile( + circuit, + CompileConfig { + mode, + target: CompileTarget::Logical, + resource_policy: ResourcePolicy::default(), + }, + ) + .unwrap() + .circuit + }); + } +} + +#[test] +fn circuit_definition_order_preserves_device_layout_and_basis() { + use crate::compile::resource::ResourcePolicy; + use crate::compile::{CompileConfig, CompileMode, CompileTarget, DeviceCompileTarget, compile}; + use crate::device::{Device, Layout}; + + let mut inner = Circuit::from_qubits(vec![Qubit::new(2), Qubit::new(0)]).unwrap(); + inner.x(Qubit::new(2)).unwrap(); + inner.cx(Qubit::new(2), Qubit::new(0)).unwrap(); + let gate = UnitaryGate::new("device_unitary", 2, 0) + .with_circuit(Arc::new(FrozenCircuit::new(inner))) + .unwrap(); + let mut circuit = Circuit::new(2); + circuit + .unitary(gate, vec![Qubit::new(0), Qubit::new(1)]) + .unwrap(); + let device = Device::line("unitary-ports", 2) + .unwrap() + .with_native_gates(vec![ + StandardGate::X.into(), + StandardGate::H.into(), + StandardGate::CX.into(), + ]) + .unwrap(); + let layout = Layout::from_pairs(&[(0, 1), (1, 0)], 2).unwrap(); + // Internal 2 -> argument 1 -> physical 0. The CX direction is native on + // the device even though its logical arguments are in reverse order. + let mut expected = Circuit::new(2); + expected.x(Qubit::new(0)).unwrap(); + expected.cx(Qubit::new(0), Qubit::new(1)).unwrap(); + for mode in [CompileMode::Normal, CompileMode::Enhanced] { + let result = compile( + &circuit, + CompileConfig { + mode, + target: CompileTarget::Device(DeviceCompileTarget { + device: device.clone(), + initial_layout: Some(layout.clone()), + seed: Some(47), + }), + resource_policy: ResourcePolicy::default(), + }, + ) + .unwrap(); + device.validate_circuit(&result.circuit).unwrap(); + let metadata = result.device_metadata.unwrap(); + assert_eq!(metadata.initial_layout, layout); + assert_eq!(metadata.final_layout, layout); + assert_matrix_approx_eq( + &circuit_to_matrix(&result.circuit, None).unwrap(), + &circuit_to_matrix(&expected, None).unwrap(), + 1e-9, + ); + } +} + +#[test] +fn circuit_definition_order_preserves_nested_symbolic_parameters_and_phase() { + use crate::circuit::symbolic_matrix::{circuit_to_symbolic_matrix, evaluate_symbolic_matrix}; + let mut inner = Circuit::from_qubits(vec![Qubit::new(2), Qubit::new(0)]).unwrap(); + inner.rx(Qubit::new(2), Parameter::symbol("theta")).unwrap(); + inner.cx(Qubit::new(2), Qubit::new(0)).unwrap(); + inner.set_global_phase(Parameter::symbol("theta") * 0.25); + let gate = UnitaryGate::new("symbolic_unitary", 2, 1) + .with_circuit(Arc::new(FrozenCircuit::new(inner))) + .unwrap(); + let mut middle = Circuit::from_qubits(vec![Qubit::new(7), Qubit::new(3)]).unwrap(); + middle + .unitary_with_params( + gate, + vec![Qubit::new(3), Qubit::new(7)], + [Parameter::symbol("alpha").into()], + ) + .unwrap(); + let wrapper = middle.to_gate("storage_order_wrapper").unwrap(); + let mut circuit = Circuit::new(3); + circuit + .append( + wrapper, + [Qubit::new(2), Qubit::new(0)], + [crate::circuit::ParameterValue::Param( + Parameter::symbol("t") + 0.1_f64, + )], + None, + ) + .unwrap(); + let mut expected = Circuit::new(3); + expected + .rx(Qubit::new(2), Parameter::symbol("t") + 0.1) + .unwrap(); + expected.cx(Qubit::new(2), Qubit::new(0)).unwrap(); + expected.set_global_phase((Parameter::symbol("t") + 0.1) * 0.25); + + // Build the nested symbolic matrix before expansion, exercising numeric-order + // UnitaryGate and storage-order CircuitGate caches in the same expression. + let symbolic = circuit_to_symbolic_matrix(&circuit, None).unwrap(); + let expanded = expand_definitions(&circuit).unwrap(); + let expanded_symbolic = circuit_to_symbolic_matrix(&expanded, None).unwrap(); + for value in [0.0, 0.37, -0.51] { + let bindings = Some(HashMap::from([("t", value)])); + let expected = + circuit_to_matrix(&expected.assign_parameters(&bindings).unwrap(), None).unwrap(); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&symbolic, &bindings).unwrap(), + &expected, + 1e-12, + ); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&expanded_symbolic, &bindings).unwrap(), + &expected, + 1e-12, + ); + let bound = circuit.assign_parameters(&bindings).unwrap(); + assert_matrix_approx_eq(&circuit_to_matrix(&bound, None).unwrap(), &expected, 1e-12); + assert_matrix_approx_eq( + &circuit_to_matrix(&expand_definitions(&bound).unwrap(), None).unwrap(), + &expected, + 1e-12, + ); + } +} diff --git a/crates/cqlib-core/src/circuit/operation.rs b/crates/cqlib-core/src/circuit/operation.rs index 4bfafecd..1d4b9ad7 100644 --- a/crates/cqlib-core/src/circuit/operation.rs +++ b/crates/cqlib-core/src/circuit/operation.rs @@ -102,6 +102,9 @@ impl Operation { /// # Returns /// /// * `Ok(Cow)` - The unitary matrix. It may be borrowed (static) or owned (computed). + /// All gate matrices put the first qubit argument at the most-significant + /// local basis bit, including circuit-backed UnitaryGate and CircuitGate + /// definitions. Full circuit matrices instead use an LSB-first basis. /// * `Err(CircuitError)` - If the operation is non-unitary or contains an /// unresolved symbolic parameter. /// diff --git a/crates/cqlib-core/src/circuit/parameter.rs b/crates/cqlib-core/src/circuit/parameter.rs index cc497cc9..18bfe78e 100644 --- a/crates/cqlib-core/src/circuit/parameter.rs +++ b/crates/cqlib-core/src/circuit/parameter.rs @@ -42,6 +42,92 @@ use std::ops::{Add, Div, Mul, Neg, Sub}; use std::str::FromStr; use symb_anafis::{Expr, Simplify, parse, symb, visitor::ExprView}; +// The engine's map() treats polynomials as opaque and visits rebuilt nodes. +// Walk the public view instead: it exposes polynomial terms, and replacement +// expressions are returned immediately without being traversed themselves. +fn substitute_parameter_expr<'a>(expr: &Expr, lookup: &impl Fn(&str) -> Option<&'a Expr>) -> Expr { + match expr.view() { + ExprView::Number(_) => expr.clone(), + ExprView::Symbol(name) => lookup(&name).cloned().unwrap_or_else(|| expr.clone()), + ExprView::Sum(children) | ExprView::Product(children) => { + let mapped: Vec<_> = children + .iter() + .map(|child| substitute_parameter_expr(child, lookup)) + .collect(); + if children.iter().zip(&mapped).all(|(a, b)| a.as_ref() == b) { + return expr.clone(); + } + + // Smart sum/product constructors can prune tiny coefficients. + // Construct a template of distinct function calls first, then use + // map() only on that template, whose children have no polynomials. + // Inserted expressions never pass through the smart constructors. + let slots = (0..mapped.len()) + .map(|i| Expr::func("__cqlib_substitution_slot", Expr::number(i as f64))) + .collect(); + let template = if matches!(expr.view(), ExprView::Sum(_)) { + Expr::sum(slots) + } else { + Expr::product(slots) + }; + let slots = match template.view() { + ExprView::Sum(slots) | ExprView::Product(slots) => slots.into_owned(), + _ => return mapped.into_iter().next().unwrap_or(template), + }; + // Match in template order so reconstruction preserves the original + // child order even if construction sorted the placeholder nodes. + let replacements: HashMap<_, _> = slots + .into_iter() + .map(|slot| (*slot).clone()) + .zip(mapped) + .collect(); + template.map(|node| { + replacements + .get(node) + .cloned() + .unwrap_or_else(|| node.clone()) + }) + } + ExprView::Div(lhs, rhs) => { + let a = substitute_parameter_expr(lhs, lookup); + let b = substitute_parameter_expr(rhs, lookup); + if &a == lhs && &b == rhs { + expr.clone() + } else { + div_parameter_expr(a, b) + } + } + ExprView::Pow(base, exponent) => { + let a = substitute_parameter_expr(base, lookup); + let b = substitute_parameter_expr(exponent, lookup); + if &a == base && &b == exponent { + expr.clone() + } else { + Expr::pow_static(a, b) + } + } + ExprView::Function { name, args } => { + let mapped: Vec<_> = args + .iter() + .map(|arg| substitute_parameter_expr(arg, lookup)) + .collect(); + if args.iter().zip(&mapped).all(|(a, b)| a.as_ref() == b) { + expr.clone() + } else { + Expr::func_multi(name, mapped) + } + } + ExprView::Derivative { inner, var, order } => { + let mapped = substitute_parameter_expr(inner, lookup); + if &mapped == inner { + expr.clone() + } else { + Expr::derivative(mapped, var, order) + } + } + } +} + // Integral coefficients in the exactly representable integer range cannot // trigger the engine's near-zero/near-one pruning. Ordinary sums of symbols // therefore avoid coefficient collection and expression rebuilding entirely. @@ -774,6 +860,7 @@ impl Parameter { /// Substitutes every occurrence of `symbol` in this expression with /// `param`, returning the resulting expression. + /// The inserted expression is not visited again; the original is unchanged. /// /// # Examples /// @@ -788,7 +875,9 @@ impl Parameter { pub fn replace(&self, symbol: &str, param: impl Into) -> Self { let p: Self = param.into(); Self { - expr: self.expr.substitute(symbol, &p.expr), + expr: substitute_parameter_expr(&self.expr, &|name| { + (name == symbol).then_some(&p.expr) + }), } } @@ -815,12 +904,13 @@ impl Parameter { .map(|_| symbol) } - /// Substitutes every bound symbol in this expression and simplifies when possible. + /// Simultaneously substitutes bound symbols and simplifies when possible. + /// + /// Replacement expressions are inserted as-is, without applying other + /// bindings to them. The result is independent of binding iteration order, + /// and the original expression is unchanged. pub fn substitute_many(&self, bindings: &HashMap) -> Self { - let mut substituted = self.clone(); - for (symbol, value) in bindings { - substituted = substituted.replace(symbol, value.clone()); - } + let substituted = self.substitute_many_simultaneous(bindings); substituted.simplify().unwrap_or(substituted) } @@ -835,15 +925,8 @@ impl Parameter { if bindings.is_empty() { return self.clone(); } - let replacements = bindings - .iter() - .map(|(symbol, replacement)| (Self::symbol(symbol).expr, &replacement.expr)) - .collect::>(); - Self::new(self.expr.map(|node| { - replacements - .iter() - .find(|(symbol, _)| node == symbol) - .map_or_else(|| node.clone(), |(_, replacement)| (*replacement).clone()) + Self::new(substitute_parameter_expr(&self.expr, &|name| { + bindings.get(name).map(|replacement| &replacement.expr) })) } diff --git a/crates/cqlib-core/src/circuit/parameter_test.rs b/crates/cqlib-core/src/circuit/parameter_test.rs index e26af669..b427a713 100644 --- a/crates/cqlib-core/src/circuit/parameter_test.rs +++ b/crates/cqlib-core/src/circuit/parameter_test.rs @@ -1319,3 +1319,265 @@ fn simultaneous_substitution_does_not_revisit_replacement_expressions() { assert_eq!(substituted.get_symbols(), HashSet::from(["b".to_string()])); } + +#[test] +fn substitution_swaps_symbols_in_either_binding_order() { + let x = Parameter::symbol("x"); + let y = Parameter::symbol("y"); + let expression: Parameter = x.clone() + y.clone() * 2.0; + let values = Some(HashMap::from([("x", 1.0), ("y", 3.0)])); + for pairs in [ + [("x", y.clone()), ("y", x.clone())], + [("y", x.clone()), ("x", y.clone())], + ] { + let bindings = pairs.into_iter().map(|(k, v)| (k.to_string(), v)).collect(); + for result in [ + expression.substitute_many(&bindings), + expression.substitute_many_simultaneous(&bindings), + ] { + assert_eq!(result.evaluate(&values).unwrap(), 5.0); + assert_eq!( + result.get_symbols(), + HashSet::from(["x".into(), "y".into()]) + ); + } + } + assert_eq!(expression.evaluate(&values).unwrap(), 7.0); +} + +#[test] +fn substitution_inserts_replacements_without_cascading() { + let x = Parameter::symbol("x"); + let y = Parameter::symbol("y"); + let z = Parameter::symbol("z"); + let values = Some(HashMap::from([("x", 2.0), ("y", 10.0), ("z", 20.0)])); + let cases = [ + ( + x.clone(), + vec![("x", y.clone() + 1.0), ("y", Parameter::from(2.0))], + 11.0, + ), + (x.clone(), vec![("x", x.clone() + 1.0)], 3.0), + ( + x.clone(), + vec![ + ("x", y.clone() + z.clone()), + ("y", Parameter::from(1.0)), + ("z", Parameter::from(2.0)), + ], + 30.0, + ), + ( + x.clone() + y.clone() * 2.0 + z.clone() * 3.0, + vec![("x", y), ("y", z), ("z", x)], + 56.0, + ), + ]; + for (expression, pairs, expected) in cases { + let bindings = pairs.into_iter().map(|(k, v)| (k.to_string(), v)).collect(); + for result in [ + expression.substitute_many(&bindings), + expression.substitute_many_simultaneous(&bindings), + ] { + assert_eq!(result.evaluate(&values).unwrap(), expected); + } + } +} + +#[test] +fn substitution_preserves_unbound_symbols_constants_and_noops() { + let expression = Parameter::symbol("x") + Parameter::symbol("y"); + let values = Some(HashMap::from([("x", 2.0), ("y", 3.0), ("z", 4.0)])); + for bindings in [ + HashMap::new(), + HashMap::from([("missing".into(), Parameter::symbol("z"))]), + ] { + assert_eq!( + expression + .substitute_many(&bindings) + .evaluate(&values) + .unwrap(), + 5.0 + ); + assert_eq!( + expression.substitute_many_simultaneous(&bindings), + expression + ); + } + let bindings = HashMap::from([ + ("x".into(), Parameter::symbol("z") + 1.0), + ("y".into(), Parameter::from(2.0)), + ]); + assert_eq!( + expression + .substitute_many(&bindings) + .evaluate(&values) + .unwrap(), + 7.0 + ); + assert_eq!( + Parameter::from(0.25).substitute_many(&bindings), + Parameter::from(0.25) + ); + assert_eq!( + expression.replace("missing", Parameter::from(0.0)), + expression + ); +} + +#[test] +fn substitution_reaches_polynomial_bases_and_nested_operations() { + let x = Parameter::symbol("x"); + let y = Parameter::symbol("y"); + let values = Some(HashMap::from([("y", 0.5)])); + let natural_polynomial: Parameter = x.clone() * x.clone() + x.clone() + 1.0; + let base = x.sin(); + let function_polynomial: Parameter = base.clone() * base.clone() + base * 2.0 + 1.0; + let nested = (natural_polynomial.clone() / (x.clone() + Parameter::from(2.0))).cos() + + x.clone().pow(y.clone() + 1.0); + let bindings = HashMap::from([("x".into(), y.clone())]); + for expression in [natural_polynomial, function_polynomial, nested] { + let expected = expression + .evaluate(&Some(HashMap::from([("x", 0.5), ("y", 0.5)]))) + .unwrap(); + for result in [ + expression.replace("x", y.clone()), + expression.substitute_many(&bindings), + expression.substitute_many_simultaneous(&bindings), + ] { + assert_eq!(result.get_symbols(), HashSet::from(["y".into()])); + assert!((result.evaluate(&values).unwrap() - expected).abs() < 1e-12); + } + } +} + +#[test] +fn substitution_does_not_revisit_rebuilt_division_or_power() { + let x = Parameter::symbol("x"); + let y = Parameter::symbol("y"); + let bindings = HashMap::from([ + ("x".into(), Parameter::symbol("z")), + ("y".into(), Parameter::from(1.0)), + ("z".into(), Parameter::from(2.0)), + ]); + for expression in [x.clone() / y.clone(), x.pow(y)] { + for result in [ + expression.substitute_many(&bindings), + expression.substitute_many_simultaneous(&bindings), + ] { + assert_eq!(result.get_symbols(), HashSet::from(["z".into()])); + assert_eq!( + result + .evaluate(&Some(HashMap::from([("z", 10.0)]))) + .unwrap(), + 10.0 + ); + } + } +} + +#[test] +fn substitution_preserves_formal_derivative_metadata() { + use symb_anafis::{Expr, visitor::ExprView}; + let expression = Parameter::new(Expr::derivative( + Parameter::symbol("x").sin().into_expr(), + "x", + 2, + )); + let bindings = HashMap::from([("x".into(), Parameter::symbol("y"))]); + for result in [ + expression.replace("x", Parameter::symbol("y")), + expression.substitute_many_simultaneous(&bindings), + ] { + let ExprView::Derivative { inner, var, order } = result.as_expr().view() else { + panic!("expected a formal derivative") + }; + assert_eq!(var, "x"); + assert_eq!(order, 2); + assert_eq!(inner.variables(), HashSet::from(["y".into()])); + } +} + +#[test] +fn substitution_preserves_small_and_near_unit_affine_coefficients() { + let x = Parameter::symbol("x"); + let y = Parameter::symbol("y"); + let bindings = HashMap::from([("x".into(), y.clone())]); + for coefficient in [ + 1e-16, + f64::from_bits(1), + 1.0 + f64::EPSILON, + 1.0 - f64::EPSILON, + ] { + let expression = x.clone() * coefficient; + for result in [ + expression.replace("x", y.clone()), + expression.substitute_many(&bindings), + expression.substitute_many_simultaneous(&bindings), + ] { + assert_eq!( + result.evaluate(&Some(HashMap::from([("y", 1.0)]))).unwrap(), + coefficient + ); + } + let cancelled = expression.clone() - expression; + assert_eq!( + cancelled + .substitute_many(&bindings) + .evaluate(&None) + .unwrap(), + 0.0 + ); + } +} + +#[test] +fn substitution_agrees_with_independent_numeric_composition() { + use rand::{Rng, SeedableRng, rngs::StdRng}; + fn expression(rng: &mut StdRng, depth: u32) -> Parameter { + if depth == 0 { + return Parameter::symbol(["x", "y", "z"][rng.random_range(0..3)]); + } + let lhs = expression(rng, depth - 1); + match rng.random_range(0..5) { + 0 => lhs + expression(rng, depth - 1), + 1 => lhs * expression(rng, depth - 1), + 2 => lhs.sin(), + 3 => lhs.pow(2.0), + _ => lhs / (expression(rng, depth - 1).pow(2.0) + 1.0), + } + } + let mut rng = StdRng::seed_from_u64(0xB03); + for _ in 0..64 { + let original = expression(&mut rng, 3); + let bindings = ["x", "y", "z"] + .into_iter() + .map(|name| (name.to_string(), expression(&mut rng, 1))) + .collect::>(); + let values = HashMap::from([ + ("x", rng.random_range(-0.8..0.8)), + ("y", rng.random_range(-0.8..0.8)), + ("z", rng.random_range(-0.8..0.8)), + ]); + let composed = bindings + .iter() + .map(|(name, replacement)| { + ( + name.as_str(), + replacement.evaluate(&Some(values.clone())).unwrap(), + ) + }) + .collect(); + let expected = original.evaluate(&Some(composed)).unwrap(); + for result in [ + original.substitute_many(&bindings), + original.substitute_many_simultaneous(&bindings), + ] { + let actual = result.evaluate(&Some(values.clone())).unwrap(); + assert!( + (actual - expected).abs() <= 1e-10 * expected.abs().max(1.0), + "{original}: actual={actual}, expected={expected}" + ); + } + } +} diff --git a/crates/cqlib-core/src/circuit/symbolic_matrix/gate.rs b/crates/cqlib-core/src/circuit/symbolic_matrix/gate.rs index 3cfb3983..27530d0e 100644 --- a/crates/cqlib-core/src/circuit/symbolic_matrix/gate.rs +++ b/crates/cqlib-core/src/circuit/symbolic_matrix/gate.rs @@ -155,10 +155,7 @@ pub fn standard_gate_symbolic_matrix( [o.clone(), z.clone()], [z, SymbolicComplex::exp_i(params[0].clone())] ], - StandardGate::GPhase => { - let phase = SymbolicComplex::exp_i(params[0].clone()); - ndarray::array![[phase.clone(), z.clone()], [z, phase]] - } + StandardGate::GPhase => ndarray::array![[SymbolicComplex::exp_i(params[0].clone())]], StandardGate::RXX => { let c = cos_half(¶ms[0]); let s = neg_i_times(sin_half(¶ms[0])); @@ -430,7 +427,9 @@ pub fn circuit_to_symbolic_matrix( flush_numeric_run(&mut matrix, &mut numeric_run, &full_bits)?; let matrix = matrix.get_or_insert_with(|| symbolic_eye(dim)); let base = standard_gate_symbolic_matrix(*mc_gate.base_gate(), ¶ms)?; - let gate_matrix = control_matrix(&base, mc_gate.num_ctrl_qubits()); + // The base matrix already includes any inherent controls. + let added_controls = mc_gate.num_qubits() - mc_gate.base_gate().num_qubits(); + let gate_matrix = control_matrix(&base, added_controls); apply_gate_to_matrix(matrix, &gate_matrix, &reversed_bits)?; } } @@ -546,7 +545,10 @@ pub fn circuit_to_symbolic_matrix( .circuit() .assign_parameters(&Some(bindings)) .map_err(|_| CircuitError::SymbolicParameterError)?; - let sub_matrix = circuit_to_matrix(&sub_circuit, None)?; + // Match CircuitGate argument positions to the inner storage order. + let order: Vec = + sub_circuit.qubits().iter().map(|q| q.index()).collect(); + let sub_matrix = circuit_to_matrix(&sub_circuit, Some(&order))?; push_numeric_gate( &mut matrix, &mut numeric_run, diff --git a/crates/cqlib-core/src/circuit/symbolic_matrix/gate_test.rs b/crates/cqlib-core/src/circuit/symbolic_matrix/gate_test.rs index 7d5d88ca..82fb2379 100644 --- a/crates/cqlib-core/src/circuit/symbolic_matrix/gate_test.rs +++ b/crates/cqlib-core/src/circuit/symbolic_matrix/gate_test.rs @@ -11,7 +11,7 @@ // that they have been altered from the originals. use super::*; -use crate::circuit::gate::{FrozenCircuit, UnitaryGate}; +use crate::circuit::gate::{FrozenCircuit, MCGate, UnitaryGate}; use crate::circuit::symbolic_matrix::gate::{ apply_gate_to_matrix, apply_standard_gate_to_matrix, circuit_to_symbolic_matrix, standard_gate_symbolic_matrix, @@ -36,6 +36,117 @@ fn assert_standard_gate_matches(gate: StandardGate, params: &[f64]) { assert_matrix_approx_eq(&evaluated, expected.as_ref(), 1e-10); } +#[test] +fn test_gphase_symbolic_gate_and_controls_match_scalar_formula() { + let theta = Parameter::symbol("theta"); + for controls in 0..=2 { + let base = + standard_gate_symbolic_matrix(StandardGate::GPhase, std::slice::from_ref(&theta)) + .unwrap(); + assert_eq!(base.dim(), (1, 1)); + let matrix = control_matrix(&base, controls); + let dim = 1 << controls; + assert_eq!(matrix.dim(), (dim, dim)); + for angle in [0.0, 0.37, -0.51, PI] { + let bindings = HashMap::from([("theta", angle)]); + let actual = evaluate_symbolic_matrix(&matrix, &Some(bindings)).unwrap(); + let mut expected = Array2::eye(dim); + expected[[dim - 1, dim - 1]] = Complex64::from_polar(1.0, angle); + assert_matrix_approx_eq(&actual, &expected, 1e-12); + } + } +} + +#[test] +fn test_gphase_symbolic_markers_accumulate_across_numeric_runs() { + for num_qubits in [0, 1, 3] { + let mut circuit = Circuit::new(num_qubits); + if num_qubits > 0 { + circuit.h(Qubit::new(0)).unwrap(); + } + for phase in [Parameter::from(-0.21), Parameter::symbol("theta")] { + circuit + .append( + Instruction::Standard(StandardGate::GPhase), + [] as [Qubit; 0], + [phase.into()], + None, + ) + .unwrap(); + } + if num_qubits > 0 { + circuit + .rz(Qubit::new(0), Parameter::symbol("alpha")) + .unwrap(); + } + circuit + .append( + Instruction::Standard(StandardGate::GPhase), + [] as [Qubit; 0], + [0.12.into()], + None, + ) + .unwrap(); + circuit.set_global_phase(Parameter::from(0.17)); + let symbolic = circuit_to_symbolic_matrix(&circuit, None).unwrap(); + let mut baseline = Circuit::new(num_qubits); + if num_qubits > 0 { + baseline.h(Qubit::new(0)).unwrap(); + baseline.rz(Qubit::new(0), Parameter::from(0.23)).unwrap(); + } + let baseline = circuit_to_matrix(&baseline, None).unwrap(); + for theta in [0.0, 0.37, -0.51, PI] { + let bindings = HashMap::from([("theta", theta), ("alpha", 0.23)]); + let actual = evaluate_symbolic_matrix(&symbolic, &Some(bindings)).unwrap(); + let expected = &baseline * Complex64::from_polar(1.0, theta - 0.21 + 0.12 + 0.17); + assert_matrix_approx_eq(&actual, &expected, 1e-12); + } + } +} + +#[test] +fn test_gphase_symbolic_controlled_circuit_with_constant_and_symbolic_params() { + for controls in [vec![], vec![2usize], vec![2, 0]] { + for phase in [Parameter::from(0.37), Parameter::symbol("theta")] { + let mut circuit = Circuit::new(3); + circuit.h(Qubit::new(1)).unwrap(); + circuit + .multi_control( + StandardGate::GPhase, + controls.iter().map(|&q| Qubit::new(q as u32)), + [] as [Qubit; 0], + [phase.into()], + ) + .unwrap(); + for order in [[0, 1, 2], [2, 0, 1]] { + let spectator_bit = order.iter().position(|&q| q == 1).unwrap(); + let mask = controls.iter().fold(0, |mask, control| { + mask | (1 << order.iter().position(|q| q == control).unwrap()) + }); + let symbolic = circuit_to_symbolic_matrix(&circuit, Some(&order)).unwrap(); + let bindings = HashMap::from([("theta", 0.37)]); + let actual = evaluate_symbolic_matrix(&symbolic, &Some(bindings)).unwrap(); + let mut expected = Array2::zeros((8, 8)); + for row in 0..8 { + for col in 0..8 { + if (row ^ col) & !(1 << spectator_bit) == 0 { + let sign = if row & col & (1 << spectator_bit) != 0 { + -1.0 + } else { + 1.0 + }; + let angle = if row & mask == mask { 0.37 } else { 0.0 }; + expected[[row, col]] = + Complex64::from_polar(1.0, angle) * sign / std::f64::consts::SQRT_2; + } + } + } + assert_matrix_approx_eq(&actual, &expected, 1e-12); + } + } + } +} + #[test] fn test_fixed_standard_gate_matches_numeric_matrix() { let symbolic = standard_gate_symbolic_matrix(StandardGate::H, &[]).unwrap(); @@ -795,6 +906,202 @@ fn test_substitution_cross_dependency() { assert_matrix_approx_eq(&evaluated, expected.as_ref(), 1e-10); } +#[test] +fn test_circuit_gate_storage_order_constant_symbolic_matrix() { + for order in [[7, 2], [1, 0], [2, 7], [0, 1]] { + let order = order.map(Qubit::new); + let mut inner = Circuit::from_qubits(order.to_vec()).unwrap(); + inner.x(order[0]).unwrap(); + let gate = inner.to_gate("StoredOrderX").unwrap(); + for targets in [[4, 9], [9, 4]] { + let targets = targets.map(Qubit::new); + let mut circuit = Circuit::from_qubits(vec![Qubit::new(9), Qubit::new(4)]).unwrap(); + circuit.append(gate.clone(), targets, [], None).unwrap(); + let mut expected = Circuit::from_qubits(vec![Qubit::new(9), Qubit::new(4)]).unwrap(); + expected.x(targets[0]).unwrap(); + for output_order in [None, Some([9, 4].as_slice())] { + let symbolic = circuit_to_symbolic_matrix(&circuit, output_order).unwrap(); + let evaluated = evaluate_symbolic_matrix(&symbolic, &None).unwrap(); + assert_matrix_approx_eq( + &evaluated, + &circuit_to_matrix(&expected, output_order).unwrap(), + 1e-10, + ); + } + } + } +} + +#[test] +fn test_circuit_gate_storage_order_nested_parameterized_inverse() { + for middle_order in [None, Some([0, 1]), Some([5, 1])] { + let mut body = Circuit::from_qubits(vec![Qubit::new(7), Qubit::new(2)]).unwrap(); + body.rx(Qubit::new(7), Parameter::symbol("theta")).unwrap(); + body.cx(Qubit::new(7), Qubit::new(2)).unwrap(); + if let Some(order) = middle_order { + let gate = body.to_gate("Inner").unwrap(); + body = Circuit::from_qubits(order.map(Qubit::new).to_vec()).unwrap(); + body.append( + gate, + order.map(Qubit::new), + [ParameterValue::from(Parameter::symbol("theta"))], + None, + ) + .unwrap(); + } + let mut outer = Circuit::new(2); + outer + .append( + body.to_gate("Body").unwrap(), + [Qubit::new(0), Qubit::new(1)], + [ParameterValue::from(Parameter::symbol("x"))], + None, + ) + .unwrap(); + for inverse in [false, true] { + let circuit = if inverse { + outer.inverse().unwrap() + } else { + outer.clone() + }; + let symbolic = circuit_to_symbolic_matrix(&circuit, None).unwrap(); + for angle in [0.37, -0.61] { + let bindings = Some(HashMap::from([("x", angle)])); + let bound = circuit.assign_parameters(&bindings).unwrap(); + let mut expected = Circuit::new(2); + expected.rx(Qubit::new(0), angle).unwrap(); + expected.cx(Qubit::new(0), Qubit::new(1)).unwrap(); + if inverse { + expected = expected.inverse().unwrap(); + } + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&symbolic, &bindings).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&bound, None).unwrap(), + &expected_matrix, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&bound.decompose().unwrap(), None).unwrap(), + &expected_matrix, + 1e-10, + ); + } + } + } +} + +#[test] +fn test_circuit_gate_storage_order_cache_isolation() { + use crate::circuit::gate::CircuitGate; + + let mut inner = Circuit::from_qubits(vec![Qubit::new(7), Qubit::new(2)]).unwrap(); + inner.rx(Qubit::new(7), Parameter::symbol("theta")).unwrap(); + inner.cx(Qubit::new(7), Qubit::new(2)).unwrap(); + let frozen = FrozenCircuit::new(inner.clone()); + let gate = CircuitGate::new("Cached", frozen.clone()).unwrap(); + // Materializing either cache first must not determine the other's basis. + for gate_first in [false, true] { + let frozen = if gate_first { + FrozenCircuit::new(inner.clone()) + } else { + frozen.clone() + }; + let gate = if gate_first { + CircuitGate::new("Cached", frozen.clone()).unwrap() + } else { + gate.clone() + }; + if gate_first { + gate.symbolic_matrix().unwrap(); + } else { + frozen.symbolic_matrix().unwrap(); + } + let gate_matrix = gate.symbolic_matrix().unwrap(); + assert!(Arc::ptr_eq( + &gate_matrix, + &gate.clone().symbolic_matrix().unwrap() + )); + assert_eq!(frozen, FrozenCircuit::new(inner.clone())); + let bindings = Some(HashMap::from([("theta", 0.37)])); + let assigned = inner.assign_parameters(&bindings).unwrap(); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&gate_matrix, &bindings).unwrap(), + &circuit_to_matrix(&assigned, Some(&[7, 2])).unwrap(), + 1e-10, + ); + let default_matrix = frozen.symbolic_matrix().unwrap(); + assert!(Arc::ptr_eq( + &default_matrix, + &frozen.symbolic_matrix().unwrap() + )); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&default_matrix, &bindings).unwrap(), + &circuit_to_matrix(&assigned, None).unwrap(), + 1e-10, + ); + // Circuit-backed UnitaryGate keeps its existing numeric-sorted basis. + let unitary = UnitaryGate::new("Frozen", 2, 1) + .with_circuit(Arc::new(frozen)) + .unwrap(); + let mut unitary_circuit = Circuit::new(2); + unitary_circuit + .append( + Instruction::UnitaryGate(Box::new(unitary)), + [Qubit::new(0), Qubit::new(1)], + [ParameterValue::from(Parameter::symbol("x"))], + None, + ) + .unwrap(); + let mut unitary_expected = Circuit::new(2); + unitary_expected.rx(Qubit::new(1), 0.37).unwrap(); + unitary_expected.cx(Qubit::new(1), Qubit::new(0)).unwrap(); + let bindings = Some(HashMap::from([("x", 0.37)])); + let expected = circuit_to_matrix(&unitary_expected, None).unwrap(); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix( + &circuit_to_symbolic_matrix(&unitary_circuit, None).unwrap(), + &bindings, + ) + .unwrap(), + &expected, + 1e-10, + ); + assert_matrix_approx_eq( + &circuit_to_matrix(&unitary_circuit.assign_parameters(&bindings).unwrap(), None) + .unwrap(), + &expected, + 1e-10, + ); + let mut circuit = Circuit::new(2); + let mut expected = Circuit::new(2); + for (targets, scale) in [([0, 1], 1.0), ([1, 0], 2.0)] { + circuit + .append( + Instruction::CircuitGate(Box::new(gate.clone())), + targets.map(Qubit::new), + [ParameterValue::from(Parameter::symbol("x") * scale)], + None, + ) + .unwrap(); + expected.rx(Qubit::new(targets[0]), 0.37 * scale).unwrap(); + expected + .cx(Qubit::new(targets[0]), Qubit::new(targets[1])) + .unwrap(); + } + let symbolic = circuit_to_symbolic_matrix(&circuit, None).unwrap(); + assert_matrix_approx_eq( + &evaluate_symbolic_matrix(&symbolic, &bindings).unwrap(), + &circuit_to_matrix(&expected, None).unwrap(), + 1e-10, + ); + } +} + #[test] fn test_frozen_circuit_symbolic_matrix_cache_reuses_arc() { let mut inner = Circuit::new(1); @@ -1175,6 +1482,82 @@ fn test_mcgate_parametric_symbolic_matrix_matches_numeric() { assert_matrix_approx_eq(&evaluated, &expected, 1e-10); } +#[test] +fn test_mcgate_zero_added_controls_preserves_symbolic_controlled_rotations() { + for gate in [StandardGate::CRX, StandardGate::CRY, StandardGate::CRZ] { + let mut circuit = Circuit::new(2); + circuit + .append( + Instruction::McGate(Box::new(MCGate::new(0, gate))), + [Qubit::new(0), Qubit::new(1)], + [ParameterValue::from(Parameter::symbol("theta"))], + None, + ) + .unwrap(); + let symbolic = circuit_to_symbolic_matrix(&circuit, Some(&[1, 0])).unwrap(); + assert_eq!(symbolic.dim(), (4, 4)); + for theta in [0.0, 0.37, -0.51, PI] { + let bindings = HashMap::from([("theta", theta)]); + let actual = evaluate_symbolic_matrix(&symbolic, &Some(bindings)).unwrap(); + let expected = gate.matrix(&[theta]).unwrap(); + assert_matrix_approx_eq(&actual, expected.as_ref(), 1e-12); + } + } +} + +#[test] +fn test_mcgate_inherent_controls_symbolic_matrix_matches_bound_numeric() { + for gate in [StandardGate::CRX, StandardGate::CRY, StandardGate::CRZ] { + for added_controls in 0..=2 { + let mc_gate = MCGate::new(added_controls, gate); + assert_eq!(mc_gate.num_ctrl_qubits(), added_controls as usize + 1); + let num_qubits = mc_gate.num_qubits() + 1; + let qubit_ids = [9usize, 2, 17, 5, 30][..num_qubits].to_vec(); + let qubits: Vec = qubit_ids.iter().map(|&id| Qubit::new(id as u32)).collect(); + // Leave a spectator between gate operands in the circuit's stored order. + let operands: Vec = qubits + .iter() + .enumerate() + .filter_map(|(index, &qubit)| (index != 1).then_some(qubit)) + .collect(); + for reverse in [false, true] { + let mut operands = operands.clone(); + if reverse { + operands.reverse(); + } + let mut circuit = Circuit::from_qubits(qubits.clone()).unwrap(); + circuit.h(operands[0]).unwrap(); + circuit.x(qubits[1]).unwrap(); + circuit + .append( + Instruction::McGate(Box::new(mc_gate.clone())), + operands.clone(), + [ParameterValue::from(Parameter::symbol("theta"))], + None, + ) + .unwrap(); + circuit.rz(*operands.last().unwrap(), 0.23).unwrap(); + for order in [None, Some(qubit_ids.as_slice())] { + let symbolic = circuit_to_symbolic_matrix(&circuit, order).unwrap(); + let dim = 1 << num_qubits; + assert_eq!(symbolic.dim(), (dim, dim)); + for theta in [0.0, 0.37, -0.51, PI] { + let bindings = HashMap::from([("theta", theta)]); + let actual = + evaluate_symbolic_matrix(&symbolic, &Some(bindings.clone())).unwrap(); + let bound = circuit.assign_parameters(&Some(bindings)).unwrap(); + let expected = circuit_to_matrix(&bound, order).unwrap(); + assert_matrix_approx_eq(&actual, &expected, 1e-12); + let constant_symbolic = circuit_to_symbolic_matrix(&bound, order).unwrap(); + let constant = evaluate_symbolic_matrix(&constant_symbolic, &None).unwrap(); + assert_matrix_approx_eq(&constant, &expected, 1e-12); + } + } + } + } + } +} + #[test] fn test_numeric_path_matches_symbolic_path_for_no_param_gates() { // H, CX, X2P — all non-parametric gates should take the numeric fast path diff --git a/crates/cqlib-core/src/compile/commutation/commutation_test.rs b/crates/cqlib-core/src/compile/commutation/commutation_test.rs index 15deebd2..4cea713b 100644 --- a/crates/cqlib-core/src/compile/commutation/commutation_test.rs +++ b/crates/cqlib-core/src/compile/commutation/commutation_test.rs @@ -80,6 +80,215 @@ fn assert_pi_phase(result: CommutationResult) { assert!((phase.evaluate(&None).unwrap() - PI).abs() < 1e-10); } +#[test] +fn circuit_backed_unitary_matrix_fallback_respects_targets_and_parameters() { + for order in [vec![0, 1], vec![1, 0], vec![7, 2], vec![9, 2, 7]] { + let mut numeric_order = order.clone(); + numeric_order.sort_unstable(); + let qubits: Vec = [5, 1, 8][..order.len()] + .iter() + .copied() + .map(Qubit::new) + .collect(); + for active in 0..order.len() { + for parameterized in [false, true] { + let mut body = + Circuit::from_qubits(order.iter().copied().map(Qubit::new).collect()).unwrap(); + body.rx( + Qubit::new(numeric_order[active]), + if parameterized { + Parameter::symbol("theta") + } else { + Parameter::from(0.7) + }, + ) + .unwrap(); + body.set_global_phase(Parameter::from(0.19)); + let instruction = Instruction::UnitaryGate(Box::new( + UnitaryGate::new("rx_one_port", order.len() as u16, u16::from(parameterized)) + .with_circuit(std::sync::Arc::new(FrozenCircuit::new(body))) + .unwrap(), + )); + let params = if parameterized { + vec![Parameter::from(0.7)] + } else { + vec![] + }; + let z = Instruction::Standard(StandardGate::Z); + for checker in [ + CommutationChecker::builtin(), + CommutationChecker::with_config(CommutationConfig { + enable_rule_oracle: false, + ..CommutationConfig::default() + }), + ] { + for (position, &qubit) in qubits.iter().enumerate() { + for result in [ + checker.check(&instruction, &qubits, ¶ms, &z, &[qubit], &[]), + checker.check(&z, &[qubit], &[], &instruction, &qubits, ¶ms), + ] { + if position == active { + assert!(result.is_none(), "order={order:?}, active={active}"); + } else { + assert_exact(result); + } + } + } + } + if parameterized { + assert!( + super::matrix::matrix_commutation( + &instruction, + &qubits, + &[Parameter::symbol("unresolved")], + &z, + &[qubits[active]], + &[], + 4, + ) + .is_none() + ); + } + assert!( + super::matrix::matrix_commutation( + &instruction, + &qubits, + ¶ms, + &z, + &[qubits[active]], + &[], + order.len() - 1, + ) + .is_none() + ); + } + } + } +} + +#[test] +fn circuit_backed_unitary_matrix_fallback_preserves_global_phase_proof() { + let mut body = Circuit::new(2); + body.x(Qubit::new(0)).unwrap(); + body.set_global_phase(Parameter::from(0.19)); + let instruction = Instruction::UnitaryGate(Box::new( + UnitaryGate::new("x_first", 2, 0) + .with_circuit(std::sync::Arc::new(FrozenCircuit::new(body))) + .unwrap(), + )); + let qubits = [Qubit::new(5), Qubit::new(1)]; + let z = Instruction::Standard(StandardGate::Z); + let checker = CommutationChecker::builtin(); + assert_pi_phase(checker.check(&instruction, &qubits, &[], &z, &[qubits[0]], &[])); + assert_pi_phase(checker.check(&z, &[qubits[0]], &[], &instruction, &qubits, &[])); + assert_exact(checker.check(&instruction, &qubits, &[], &z, &[qubits[1]], &[])); +} + +#[test] +fn circuit_gate_operation_matrix_does_not_certify_false_commutation() { + let mut body = Circuit::new(2); + body.rx(Qubit::new(0), Parameter::from(0.37)).unwrap(); + let instruction = body.to_gate("rx_first").unwrap(); + let qubits = [Qubit::new(0), Qubit::new(1)]; + let checker = CommutationChecker::with_config(CommutationConfig { + enable_rule_oracle: false, + ..CommutationConfig::default() + }); + assert!( + checker + .check( + &instruction, + &qubits, + &[], + &Instruction::Standard(StandardGate::Z), + &[qubits[0]], + &[], + ) + .is_none() + ); +} + +#[test] +fn circuit_gate_matrix_fallback_respects_targets_parameters_and_phase() { + for order in [[0, 1], [1, 0], [7, 2]] { + for parameterized in [false, true] { + let mut body = Circuit::from_qubits(order.map(Qubit::new).to_vec()).unwrap(); + body.rx( + Qubit::new(order[0]), + if parameterized { + Parameter::symbol("theta") + } else { + Parameter::from(0.37) + }, + ) + .unwrap(); + let instruction = body.to_gate("rx_first").unwrap(); + let params = if parameterized { + vec![Parameter::from(0.37)] + } else { + vec![] + }; + let qubits = [Qubit::new(5), Qubit::new(1)]; + for checker in [ + CommutationChecker::builtin(), + CommutationChecker::with_config(CommutationConfig { + enable_rule_oracle: false, + ..CommutationConfig::default() + }), + ] { + for (position, &qubit) in qubits.iter().enumerate() { + let z = Instruction::Standard(StandardGate::Z); + for result in [ + checker.check(&instruction, &qubits, ¶ms, &z, &[qubit], &[]), + checker.check(&z, &[qubit], &[], &instruction, &qubits, ¶ms), + ] { + if position == 0 { + assert!(result.is_none()); + } else { + assert_exact(result); + } + } + } + } + if parameterized { + assert!( + super::matrix::matrix_commutation( + &instruction, + &qubits, + &[Parameter::symbol("unresolved")], + &Instruction::Standard(StandardGate::Z), + &[qubits[0]], + &[], + 4, + ) + .is_none() + ); + } + } + } + let mut body = Circuit::new(2); + body.x(Qubit::new(0)).unwrap(); + let x = body.to_gate("x_first").unwrap(); + let z = Instruction::Standard(StandardGate::Z); + let checker = CommutationChecker::builtin(); + assert_pi_phase(checker.check( + &x, + &[Qubit::new(0), Qubit::new(1)], + &[], + &z, + &[Qubit::new(0)], + &[], + )); + assert_pi_phase(checker.check( + &z, + &[Qubit::new(0)], + &[], + &x, + &[Qubit::new(0), Qubit::new(1)], + &[], + )); +} + #[test] fn identity_commutes_exactly() { let checker = CommutationChecker::builtin(); diff --git a/crates/cqlib-core/src/compile/commutation/matrix.rs b/crates/cqlib-core/src/compile/commutation/matrix.rs index 448d896e..c50bd939 100644 --- a/crates/cqlib-core/src/compile/commutation/matrix.rs +++ b/crates/cqlib-core/src/compile/commutation/matrix.rs @@ -101,7 +101,7 @@ pub fn matrix_commutation( /// /// `gate_positions` maps each gate-local qubit to its position in the sorted /// combined support used by [`matrix_commutation`]. Bit extraction follows the -/// existing circuit matrix convention: earlier qubits correspond to more +/// gate-local matrix convention: earlier qubits correspond to more /// significant bits in the local basis index. fn expand_unitary( matrix: Cow<'_, Array2>, diff --git a/crates/cqlib-core/src/compile/knowledge/matcher_test.rs b/crates/cqlib-core/src/compile/knowledge/matcher_test.rs index ead1817e..ee37c6da 100644 --- a/crates/cqlib-core/src/compile/knowledge/matcher_test.rs +++ b/crates/cqlib-core/src/compile/knowledge/matcher_test.rs @@ -319,6 +319,56 @@ fn instantiate_target_rejects_unbound_rewrite_references() { ); } +#[test] +fn instantiate_target_substitutes_overlapping_symbol_names_simultaneously() { + let a = Parameter::symbol("a"); + let b = Parameter::symbol("b"); + let rule = Rule::new( + "overlapping_parameters", + vec![ + RuleItem::standard( + StandardGate::RZ, + &[0], + vec![ParameterValue::Param(a.clone())], + ), + RuleItem::standard( + StandardGate::RZ, + &[0], + vec![ParameterValue::Param(b.clone())], + ), + ], + vec![RuleItem::standard( + StandardGate::RZ, + &[0], + vec![ParameterValue::Param(a.clone() + b.clone() * 2.0)], + )], + ); + let instruction = Instruction::Standard(StandardGate::RZ); + let qubits = [Qubit::new(0)]; + let first = [b]; + let second = [a]; + let operations = [ + ConcreteOperationView::new(&instruction, &qubits, &first), + ConcreteOperationView::new(&instruction, &qubits, &second), + ]; + let bindings = rule_matches_operations(&rule, &operations) + .unwrap() + .unwrap(); + let replacements = instantiate_target(&rule.target, &bindings).unwrap(); + let ParameterValue::Param(result) = &replacements[0].params[0] else { + panic!("expected a symbolic parameter"); + }; + assert_eq!( + result + .evaluate(&Some(std::collections::HashMap::from([ + ("a", 1.0), + ("b", 3.0) + ]))) + .unwrap(), + 5.0 + ); +} + #[test] fn rule_matches_operations_returns_none_for_non_adjacent_policy_cases() { let rule = Rule::new( diff --git a/crates/cqlib-core/src/compile/transform/decompose/definition.rs b/crates/cqlib-core/src/compile/transform/decompose/definition.rs index 429e5f34..97fd9859 100644 --- a/crates/cqlib-core/src/compile/transform/decompose/definition.rs +++ b/crates/cqlib-core/src/compile/transform/decompose/definition.rs @@ -26,6 +26,7 @@ use crate::compile::transform::analysis::WorkflowCircuitAnalysis; use crate::compile::transform::rebuild::{CircuitRebuildContext, ClassicalRemap}; use crate::compile::transform::transformer::{PassApplicability, WorkflowPass}; use crate::compile::transform::{CircuitAnalysis, TransformOutcome, Transformer}; +use indexmap::IndexSet; use smallvec::{SmallVec, smallvec}; use std::collections::{HashMap, HashSet}; @@ -172,6 +173,8 @@ impl<'a> DefinitionExpander<'a> { self.expand_definition( gate.name(), definition.circuit(), + &definition.circuit().qubits(), + gate.signature_params(), gate.num_qubits(), gate.num_params(), operation, @@ -186,9 +189,13 @@ impl<'a> DefinitionExpander<'a> { Instruction::UnitaryGate(gate) => { if let Some(definition) = gate.circuit().as_ref() { self.changed = true; + // Circuit-backed unitaries use the cached live-symbol order, + // matching UnitaryGate::matrix_for_params. self.expand_definition( gate.label(), definition.circuit(), + &definition.qubits_in_numeric_order(), + definition.used_symbols(), gate.num_qubits() as usize, gate.num_params() as usize, operation, @@ -240,6 +247,8 @@ impl<'a> DefinitionExpander<'a> { &mut self, name: &str, definition: &Circuit, + definition_qubits: &[Qubit], + signature_params: &IndexSet, expected_qubits: usize, expected_params: usize, operation: &Operation, @@ -269,11 +278,10 @@ impl<'a> DefinitionExpander<'a> { ))); } - let definition_symbols = definition.symbols(); - if definition_symbols.len() != expected_params { + if signature_params.len() != expected_params { return Err(CompilerError::InvariantViolation(format!( - "definition '{name}' signature has {expected_params} parameters but backing circuit has {} symbols", - definition_symbols.len() + "definition '{name}' expects {expected_params} parameters but its signature has {} parameters", + signature_params.len() ))); } @@ -284,12 +292,13 @@ impl<'a> DefinitionExpander<'a> { // nested definition must receive a fresh binding table instead of a // merged parent table. Merging would make equal symbol names in nested // definitions accidentally capture outer bindings. - let mut next_symbol_bindings = HashMap::with_capacity(definition_symbols.len()); - for (symbol, value) in definition_symbols.iter().zip(resolved_params) { + // Only the call signature defines positional bindings. Live dependencies + // may omit unused formals, and the circuit registry may retain stale symbols. + let mut next_symbol_bindings = HashMap::with_capacity(signature_params.len()); + for (symbol, value) in signature_params.iter().zip(resolved_params) { next_symbol_bindings.insert(symbol.clone(), value); } - let definition_qubits = definition.qubits(); if definition_qubits.len() != expected_qubits { return Err(CompilerError::InvariantViolation(format!( "definition '{name}' signature has {expected_qubits} qubits but backing circuit has {} qubits", @@ -599,8 +608,9 @@ fn fresh_temp_symbol(index: usize, symbol: &str, occupied: &mut HashSet) #[cfg(test)] mod tests { use super::*; - use crate::circuit::ParameterValue; use crate::circuit::gate::{FrozenCircuit, UnitaryGate}; + use crate::circuit::test_utils::assert_matrix_approx_eq; + use crate::circuit::{CircuitGate, ParameterValue, circuit_to_matrix}; use crate::circuit::{ClassicalDataOp, ClassicalExpr}; use crate::compile::transform::TransformerTestExt; use ndarray::array; @@ -706,6 +716,192 @@ mod tests { assert_eq!(rz_param.evaluate(&Some(bindings)).unwrap(), 6.0); } + #[test] + fn expands_circuit_gate_using_declared_parameter_signature() { + let q0 = Qubit::new(0); + let mut expected = Circuit::new(1); + expected.rx(q0, -0.2).unwrap(); + expected.rz(q0, 0.7).unwrap(); + expected.set_global_phase(Parameter::from(-1.6)); + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + + for stale_registry in [false, true] { + for unused_formal in [false, true] { + let mut definition = Circuit::new(1); + if stale_registry { + definition.add_parameter(Parameter::symbol("stale")); + } + definition.rx(q0, Parameter::symbol("theta")).unwrap(); + definition.rz(q0, Parameter::symbol("phi")).unwrap(); + definition + .set_global_phase(Parameter::symbol("theta") - Parameter::symbol("phi") * 2.0); + + let mut signature = vec!["phi".to_string(), "theta".to_string()]; + let mut args = vec![ParameterValue::from(0.7), ParameterValue::from(-0.2)]; + if unused_formal { + signature.insert(0, "unused".to_string()); + args.insert(0, ParameterValue::from(11.0)); + } + let gate = CircuitGate::with_signature( + "declared", + FrozenCircuit::new(definition), + signature, + ) + .unwrap(); + let mut circuit = Circuit::new(1); + circuit.circuit_gate(gate, vec![q0], args).unwrap(); + + let expanded = expand_definitions(&circuit).unwrap(); + assert_eq!(expanded.operations().len(), 2); + for actual in [&circuit, &circuit.decompose().unwrap(), &expanded] { + assert_matrix_approx_eq( + &circuit_to_matrix(actual, None).unwrap(), + &expected_matrix, + 1e-12, + ); + } + } + } + } + + #[test] + fn expands_inferred_signatures_without_removed_operation_symbols() { + let q0 = Qubit::new(0); + let mut definition = Circuit::new(1); + definition.rx(q0, Parameter::symbol("removed")).unwrap(); + definition.remove_operation(0).unwrap(); + definition.rx(q0, Parameter::symbol("theta")).unwrap(); + definition.rz(q0, Parameter::symbol("phi")).unwrap(); + let circuit_gate = definition.clone().to_gate("inferred").unwrap(); + let unitary_gate = UnitaryGate::new("inferred_unitary", 1, 2) + .with_circuit(Arc::new(FrozenCircuit::new(definition))) + .unwrap(); + + let mut expected = Circuit::new(1); + expected.rx(q0, 0.7).unwrap(); + expected.rz(q0, -0.2).unwrap(); + let expected_matrix = circuit_to_matrix(&expected, None).unwrap(); + for gate in [ + circuit_gate, + Instruction::UnitaryGate(Box::new(unitary_gate)), + ] { + let mut circuit = Circuit::new(1); + circuit + .append( + gate, + [q0], + [ParameterValue::from(0.7), ParameterValue::from(-0.2)], + None, + ) + .unwrap(); + let expanded = expand_definitions(&circuit).unwrap(); + for actual in [&circuit, &expanded] { + assert_matrix_approx_eq( + &circuit_to_matrix(actual, None).unwrap(), + &expected_matrix, + 1e-12, + ); + } + } + } + + #[test] + fn expands_phase_only_definition_with_unused_formal() { + let mut definition = Circuit::new(1); + definition.set_global_phase(Parameter::symbol("phase")); + let gate = CircuitGate::with_signature( + "phase_only", + FrozenCircuit::new(definition), + ["unused".to_string(), "phase".to_string()], + ) + .unwrap(); + let mut circuit = Circuit::new(1); + circuit + .append( + Instruction::CircuitGate(Box::new(gate)), + [Qubit::new(0)], + [ParameterValue::from(11.0), ParameterValue::from(0.7)], + None, + ) + .unwrap(); + let expanded = expand_definitions(&circuit).unwrap(); + assert!(expanded.operations().is_empty()); + assert!((expanded.global_phase().evaluate(&None).unwrap() - 0.7).abs() < 1e-12); + assert_matrix_approx_eq( + &circuit_to_matrix(&expanded, None).unwrap(), + &circuit_to_matrix(&circuit, None).unwrap(), + 1e-12, + ); + } + + #[test] + fn expands_nested_signatures_with_simultaneous_symbolic_arguments() { + let q0 = Qubit::new(0); + let theta = Parameter::symbol("theta"); + let phi = Parameter::symbol("phi"); + let mut leaf = Circuit::new(1); + leaf.rx(q0, theta.clone()).unwrap(); + leaf.rz(q0, phi.clone()).unwrap(); + leaf.set_global_phase(theta.clone() + phi.clone() * 2.0); + let leaf_gate = CircuitGate::with_signature( + "leaf", + FrozenCircuit::new(leaf), + ["unused".to_string(), "phi".to_string(), "theta".to_string()], + ) + .unwrap(); + + let mut middle = Circuit::new(1); + middle + .append( + Instruction::CircuitGate(Box::new(leaf_gate)), + [q0], + [ + ParameterValue::from(11.0), + ParameterValue::from(theta.clone() + 1.0), + ParameterValue::from(phi.clone() - 0.5), + ], + None, + ) + .unwrap(); + let middle_gate = CircuitGate::with_signature( + "middle", + FrozenCircuit::new(middle), + ["phi".to_string(), "theta".to_string()], + ) + .unwrap(); + let mut circuit = Circuit::new(1); + circuit + .append( + Instruction::CircuitGate(Box::new(middle_gate)), + [q0], + [ + ParameterValue::from(theta.clone()), + ParameterValue::from(phi.clone()), + ], + None, + ) + .unwrap(); + + let mut expected = Circuit::new(1); + expected.rx(q0, theta.clone() - 0.5).unwrap(); + expected.rz(q0, phi.clone() + 1.0).unwrap(); + expected.set_global_phase(theta - 0.5 + (phi + 1.0) * 2.0); + let expanded = expand_definitions(&circuit).unwrap(); + for (theta, phi) in [(0.2, 0.8), (-0.4, 0.1)] { + let bindings = Some(HashMap::from([("theta", theta), ("phi", phi)])); + let expected_matrix = + circuit_to_matrix(&expected.assign_parameters(&bindings).unwrap(), None).unwrap(); + for actual in [&circuit, &expanded] { + let bound = actual.assign_parameters(&bindings).unwrap(); + assert_matrix_approx_eq( + &circuit_to_matrix(&bound, None).unwrap(), + &expected_matrix, + 1e-12, + ); + } + } + } + #[test] fn uses_simultaneous_parameter_substitution() { let mut inner = Circuit::new(1); @@ -894,15 +1090,27 @@ mod tests { #[test] fn expands_definitions_in_control_flow_body_with_body_local_phase() { let mut inner = Circuit::new(1); - inner.set_global_phase(Parameter::from(0.5)); + inner.set_global_phase(Parameter::symbol("phase")); inner.x(Qubit::new(0)).unwrap(); - let gate = inner.to_gate("phase_x").unwrap(); + let gate = CircuitGate::with_signature( + "phase_x", + FrozenCircuit::new(inner), + ["unused".to_string(), "phase".to_string()], + ) + .unwrap(); let mut circuit = Circuit::new(2); circuit .if_else( ClassicalExpr::bool_literal(true), - |body| body.append(gate, [Qubit::new(1)], [], None), + |body| { + body.append( + Instruction::CircuitGate(Box::new(gate)), + [Qubit::new(1)], + [ParameterValue::from(11.0), ParameterValue::from(0.5)], + None, + ) + }, |_body| Ok(()), ) .unwrap(); @@ -920,6 +1128,14 @@ mod tests { then_body[0].instruction, Instruction::Standard(StandardGate::GPhase) )); + assert!( + (op_param(&expanded, &then_body[0].params[0]) + .evaluate(&None) + .unwrap() + - 0.5) + .abs() + < 1e-12 + ); assert!(matches!( then_body[1].instruction, Instruction::Standard(StandardGate::X) diff --git a/crates/cqlib-core/src/compile/workflow_test.rs b/crates/cqlib-core/src/compile/workflow_test.rs index 689dd5b6..7a87dcd7 100644 --- a/crates/cqlib-core/src/compile/workflow_test.rs +++ b/crates/cqlib-core/src/compile/workflow_test.rs @@ -1365,6 +1365,81 @@ fn workflow_expands_circuit_gate_definitions_before_optimization() { assert!(!contains_high_level_gate(&result.circuit)); } +#[test] +fn workflow_preserves_qasm_gate_signatures_in_both_modes() { + let declarations = r#" + gate reordered(phi, theta) a, b { + rx(theta) a; + rz(phi) b; + cx a, b; + } + gate declared(unused, phi, theta) a, b { + reordered(phi, theta) a, b; + } + "#; + let mut expected = Circuit::new(2); + expected.rx(Qubit::new(1), -0.2).unwrap(); + expected.rz(Qubit::new(0), 0.7).unwrap(); + expected.cx(Qubit::new(1), Qubit::new(0)).unwrap(); + let expected_matrix = crate::circuit::circuit_to_matrix(&expected, None).unwrap(); + let basis = vec![ + Instruction::Standard(StandardGate::RX), + Instruction::Standard(StandardGate::RY), + Instruction::Standard(StandardGate::CZ), + ]; + let device = Device::bidirectional_line("declared-signatures", 2) + .unwrap() + .with_native_gates(basis.clone()) + .unwrap(); + + for call in ["reordered(0.7, -0.2)", "declared(11.0, 0.7, -0.2)"] { + let qasm2 = format!( + "OPENQASM 2.0; include \"qelib1.inc\"; {declarations} qreg q[2]; {call} q[1], q[0];" + ); + let qasm3 = format!( + "OPENQASM 3.0; include \"stdgates.inc\"; {declarations} qubit[2] q; {call} q[1], q[0];" + ); + for circuit in [ + crate::ir::qasm2_loads(&qasm2).unwrap(), + crate::ir::qasm3_loads(&qasm3).unwrap(), + ] { + for mode in [CompileMode::Normal, CompileMode::Enhanced] { + for target in [ + CompileTarget::Logical, + CompileTarget::Basis(basis.clone()), + CompileTarget::Device(DeviceCompileTarget { + device: device.clone(), + initial_layout: Some(Layout::from_pairs(&[(0, 0), (1, 1)], 2).unwrap()), + seed: Some(47), + }), + ] { + let device_target = matches!(target, CompileTarget::Device(_)); + let result = compile( + &circuit, + CompileConfig { + mode, + target, + resource_policy: ResourcePolicy::default(), + }, + ) + .unwrap(); + + assert!(result.step_changed("decompose.definitions")); + assert!(!contains_high_level_gate(&result.circuit)); + crate::circuit::test_utils::assert_matrix_approx_eq( + &crate::circuit::circuit_to_matrix(&result.circuit, None).unwrap(), + &expected_matrix, + 1e-9, + ); + if device_target { + device.validate_circuit(&result.circuit).unwrap(); + } + } + } + } + } +} + #[test] fn definition_expansion_invalidates_analysis_before_unitary_decomposition() { let q0 = Qubit::new(0); diff --git a/crates/cqlib-core/src/ir/qasm2/dump.rs b/crates/cqlib-core/src/ir/qasm2/dump.rs index 03af805a..a80aa55d 100644 --- a/crates/cqlib-core/src/ir/qasm2/dump.rs +++ b/crates/cqlib-core/src/ir/qasm2/dump.rs @@ -596,7 +596,7 @@ fn dump_unitary_gate_definition( // Body context let mut body_qubit_map = HashMap::new(); let inner_circuit = &frozen_circuit.circuit; - for (i, q) in inner_circuit.qubits().iter().enumerate() { + for (i, q) in frozen_circuit.qubits_in_numeric_order().iter().enumerate() { body_qubit_map.insert(*q, format!("q{}", i)); } diff --git a/crates/cqlib-core/src/ir/qasm3/dump.rs b/crates/cqlib-core/src/ir/qasm3/dump.rs index 914669cf..a6bd1313 100644 --- a/crates/cqlib-core/src/ir/qasm3/dump.rs +++ b/crates/cqlib-core/src/ir/qasm3/dump.rs @@ -516,6 +516,7 @@ fn dump_gate_definition(gate: &CircuitGate, output: &mut String) -> Result<(), Q for (index, qubit) in gate.circuit.circuit.qubits().iter().copied().enumerate() { qubit_map.insert(qubit, format!("q{index}")); } + dump_global_phase(&gate.circuit.circuit, output, &HashMap::new())?; dump_operations( &gate.circuit.circuit, gate.circuit.circuit.operations(), @@ -549,9 +550,10 @@ fn dump_unitary_gate_definition( .collect(); writeln!(output, "gate {label} {} {{", qubits.join(","))?; let mut qubit_map = HashMap::new(); - for (index, qubit) in frozen.circuit.qubits().iter().copied().enumerate() { + for (index, qubit) in frozen.qubits_in_numeric_order().into_iter().enumerate() { qubit_map.insert(qubit, format!("q{index}")); } + dump_global_phase(&frozen.circuit, output, &HashMap::new())?; dump_operations( &frozen.circuit, frozen.circuit.operations(), diff --git a/crates/cqlib-core/src/ir/qasm3/dump_test.rs b/crates/cqlib-core/src/ir/qasm3/dump_test.rs index 9b93c967..ddd548d5 100644 --- a/crates/cqlib-core/src/ir/qasm3/dump_test.rs +++ b/crates/cqlib-core/src/ir/qasm3/dump_test.rs @@ -862,3 +862,53 @@ fn export_preserves_indexed_constants_and_symbol_substitution() { "(angle0)" ); } + +#[test] +fn nested_gate_global_phase_round_trip() { + let circuit = qasm3_loads( + r#" + OPENQASM 3; + include "stdgates.inc"; + gate fragment(theta) q { gphase(theta); h q; } + gate middle(alpha) q { gphase(0.11); fragment(alpha) q; } + qubit q; + middle(0.37) q; + middle(-0.51) q; + "#, + ) + .unwrap(); + let output = dumps(&circuit).unwrap(); + let loaded = qasm3_loads(&output).unwrap(); + let expected = num_complex::Complex64::from_polar(1.0, 0.08); + for ((row, col), actual) in loaded.to_matrix(None).unwrap().indexed_iter() { + let reference = if row == col { + expected + } else { + num_complex::Complex64::new(0.0, 0.0) + }; + assert!((*actual - reference).norm() < 1e-12, "{output}"); + } + assert!(output.contains("gphase(theta);"), "{output}"); + assert!(output.contains("gphase(0.11);"), "{output}"); +} + +#[test] +fn unitary_definition_global_phase_round_trip() { + use crate::circuit::gate::UnitaryGate; + let mut inner = Circuit::new(1); + inner.h(Qubit::new(0)).unwrap(); + inner.set_global_phase(0.37.into()); + let expected = inner.to_matrix(None).unwrap(); + let gate = UnitaryGate::new("phase_h", 1, 0) + .with_matrix(expected.clone()) + .unwrap() + .with_circuit(Arc::new(FrozenCircuit::new(inner))) + .unwrap(); + let mut circuit = Circuit::new(1); + circuit.unitary(gate, vec![Qubit::new(0)]).unwrap(); + let output = dumps(&circuit).unwrap(); + let loaded = qasm3_loads(&output).unwrap(); + for (actual, expected) in loaded.to_matrix(None).unwrap().iter().zip(expected.iter()) { + assert!((actual - expected).norm() < 1e-12, "{output}"); + } +} diff --git a/crates/cqlib-core/src/qis/state/density_matrix.rs b/crates/cqlib-core/src/qis/state/density_matrix.rs index fbec87a1..904e9f9a 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix.rs @@ -151,7 +151,7 @@ impl DensityMatrix { Ok(()) } - /// Validates that two qubit indices are within bounds. + /// Validates that two qubit indices are within bounds and distinct. #[inline] fn validate_two_qubits(&self, q0: usize, q1: usize) -> Result<(), QisError> { if self.num_qubits == 0 { @@ -171,10 +171,15 @@ impl DensityMatrix { max: self.num_qubits.saturating_sub(1), }); } + if q0 == q1 { + return Err(QisError::InvalidParameterValue(format!( + "Duplicate qubit index {q0} in gate operation" + ))); + } Ok(()) } - /// Validates that all qubit indices in a slice are within bounds. + /// Validates that all qubit indices in a slice are within bounds and distinct. #[inline] fn validate_qubits(&self, qubits: &[usize]) -> Result<(), QisError> { if self.num_qubits == 0 { @@ -190,6 +195,13 @@ impl DensityMatrix { }); } } + for (i, &q) in qubits.iter().enumerate() { + if qubits[..i].contains(&q) { + return Err(QisError::InvalidParameterValue(format!( + "Duplicate qubit index {q} in gate operation" + ))); + } + } Ok(()) } @@ -464,30 +476,31 @@ impl DensityMatrix { let num_controls = mc_gate.num_ctrl_qubits(); let base_gate = mc_gate.base_gate(); - if num_controls == 1 { - let control = qubit_indices[0]; - let target = qubit_indices[1]; - match base_gate { - StandardGate::X => dm.apply_cx(control, target)?, - StandardGate::Y => dm.apply_cy(control, target)?, - StandardGate::Z => dm.apply_cz(control, target)?, - StandardGate::RX => dm.apply_crx(control, target, params[0])?, - StandardGate::RY => dm.apply_cry(control, target, params[0])?, - StandardGate::RZ => dm.apply_crz(control, target, params[0])?, - _ => { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError(CircuitError::NoMatrixRepresentation) - })?; - dm.apply_unitary_gate(&qubit_indices, &matrix)?; - } + match (num_controls, *base_gate) { + (0, StandardGate::GPhase) => { + dm.apply_standard_gate(*base_gate, &qubit_indices, ¶ms)? + } + (1, StandardGate::X) => dm.apply_cx(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::Y) => dm.apply_cy(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::Z) => dm.apply_cz(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::RX) => { + dm.apply_crx(qubit_indices[0], qubit_indices[1], params[0])? + } + (1, StandardGate::RY) => { + dm.apply_cry(qubit_indices[0], qubit_indices[1], params[0])? + } + (1, StandardGate::RZ) => { + dm.apply_crz(qubit_indices[0], qubit_indices[1], params[0])? + } + (2, StandardGate::X) => { + dm.apply_ccx(qubit_indices[0], qubit_indices[1], qubit_indices[2])? + } + _ => { + let matrix = mc_gate.matrix(¶ms).map_err(|_| { + QisError::CircuitError(CircuitError::NoMatrixRepresentation) + })?; + dm.apply_unitary_gate(&qubit_indices, &matrix)?; } - } else if num_controls == 2 && *base_gate == StandardGate::X { - dm.apply_ccx(qubit_indices[0], qubit_indices[1], qubit_indices[2])?; - } else { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError(CircuitError::NoMatrixRepresentation) - })?; - dm.apply_unitary_gate(&qubit_indices, &matrix)?; } } Instruction::UnitaryGate(u_gate) => { diff --git a/crates/cqlib-core/src/qis/state/density_matrix_noise.rs b/crates/cqlib-core/src/qis/state/density_matrix_noise.rs index d76fe5c0..6aff31b0 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_noise.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_noise.rs @@ -261,36 +261,26 @@ impl DensityMatrixNoise { let num_controls = mc_gate.num_ctrl_qubits(); let base_gate = mc_gate.base_gate(); - if num_controls == 1 { - let c = qs[0]; - let t = qs[1]; - match base_gate { - StandardGate::X => sim.apply_cx(c, t)?, - StandardGate::Y => sim.apply_cy(c, t)?, - StandardGate::Z => sim.apply_cz(c, t)?, - StandardGate::RX => sim.apply_crx(c, t, params[0])?, - StandardGate::RY => sim.apply_cry(c, t, params[0])?, - StandardGate::RZ => sim.apply_crz(c, t, params[0])?, - _ => { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError( - crate::circuit::CircuitError::NoMatrixRepresentation, - ) - })?; - sim.apply_unitary_gate(&qs, &matrix)?; - sim.apply_noise(*base_gate, &qs)?; - } + match (num_controls, *base_gate) { + (0, StandardGate::GPhase) => { + sim.apply_standard_gate_noise(*base_gate, &qs, ¶ms)? + } + (1, StandardGate::X) => sim.apply_cx(qs[0], qs[1])?, + (1, StandardGate::Y) => sim.apply_cy(qs[0], qs[1])?, + (1, StandardGate::Z) => sim.apply_cz(qs[0], qs[1])?, + (1, StandardGate::RX) => sim.apply_crx(qs[0], qs[1], params[0])?, + (1, StandardGate::RY) => sim.apply_cry(qs[0], qs[1], params[0])?, + (1, StandardGate::RZ) => sim.apply_crz(qs[0], qs[1], params[0])?, + (2, StandardGate::X) => sim.apply_ccx(qs[0], qs[1], qs[2])?, + _ => { + let matrix = mc_gate.matrix(¶ms).map_err(|_| { + QisError::CircuitError( + crate::circuit::CircuitError::NoMatrixRepresentation, + ) + })?; + sim.apply_unitary_gate(&qs, &matrix)?; + sim.apply_noise(*base_gate, &qs)?; } - } else if num_controls == 2 && *base_gate == StandardGate::X { - sim.apply_ccx(qs[0], qs[1], qs[2])?; - } else { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError( - crate::circuit::CircuitError::NoMatrixRepresentation, - ) - })?; - sim.apply_unitary_gate(&qs, &matrix)?; - sim.apply_noise(*base_gate, &qs)?; } } Instruction::UnitaryGate(u_gate) => { @@ -432,25 +422,28 @@ impl DensityMatrixNoise { .collect() } + /// Resolves a storage position to its noise-model qubit ID. + /// Callers validate positions before looking up noise. + fn noise_qubit(&self, position: usize) -> Qubit { + self.state + .qubit_map + .as_ref() + .and_then(|map| { + map.iter() + .find_map(|(&qubit, &index)| (index == position).then_some(qubit)) + }) + .unwrap_or_else(|| Qubit::new(position as u32)) + } + /// Applies noise channels associated with a gate operation. /// /// Looks up the noise model for errors associated with `gate` on `qubits`, /// converts them to Kraus operators, and applies them to the state. /// Supports single-qubit, two-qubit, and three-qubit gates. fn apply_noise(&mut self, gate: StandardGate, qubits: &[usize]) -> Result<(), QisError> { - let labelled_qubit = |position: usize| { - self.state - .qubit_map - .as_ref() - .and_then(|map| { - map.iter() - .find_map(|(&qubit, &index)| (index == position).then_some(qubit)) - }) - .unwrap_or_else(|| Qubit::new(position as u32)) - }; if let Some(noise_model) = &self.noise_model { if qubits.len() == 1 { - let q0 = labelled_qubit(qubits[0]); + let q0 = self.noise_qubit(qubits[0]); let key = OperationKey::new_single(gate, q0); if let Some(errors) = noise_model.get_single_qubit_errors(&key) { for error in errors { @@ -460,8 +453,8 @@ impl DensityMatrixNoise { } } } else if qubits.len() == 2 { - let q0 = labelled_qubit(qubits[0]); - let q1 = labelled_qubit(qubits[1]); + let q0 = self.noise_qubit(qubits[0]); + let q1 = self.noise_qubit(qubits[1]); if let Ok(key) = OperationKey::new_double(gate, q0, q1) && let Some(errors) = noise_model.get_two_qubit_errors(&key) { @@ -472,9 +465,9 @@ impl DensityMatrixNoise { } } } else if qubits.len() == 3 { - let q0 = labelled_qubit(qubits[0]); - let q1 = labelled_qubit(qubits[1]); - let q2 = labelled_qubit(qubits[2]); + let q0 = self.noise_qubit(qubits[0]); + let q1 = self.noise_qubit(qubits[1]); + let q2 = self.noise_qubit(qubits[2]); if let Ok(_key) = OperationKey::new_triple(gate, q0, q1, q2) { // Current noise model struct doesn't have `get_three_qubit_errors`. // Ready for future extension without silently failing or panicking. @@ -1042,12 +1035,21 @@ impl DensityMatrixNoise { /// /// # Arguments /// - /// * `qubits` - Indices of qubits to measure. + /// * `qubits` - Storage positions on which to apply readout errors, not circuit + /// qubit IDs. For a circuit-bound state, each position selects the noise + /// configured for `circuit.qubits()[position]`. Unbound states use identity IDs. /// /// # Returns /// /// A vector of probabilities for all 2^n computational basis states, where n - /// is the total number of qubits in the simulator. + /// is the total number of qubits in the simulator. The array retains storage + /// order regardless of the order or number of selected positions. An empty + /// selection returns the probabilities before readout noise. + /// + /// # Errors + /// + /// Returns [`QisError::IndexOutOfBounds`] for a position outside `0..num_qubits`, + /// including when no noise model is configured. /// /// # Examples /// @@ -1063,10 +1065,20 @@ impl DensityMatrixNoise { /// // For Bell state |Φ+⟩, P(|00⟩) ≈ 0.5 and P(|11⟩) ≈ 0.5 /// ``` pub fn probabilities_with_readout(&self, qubits: &[usize]) -> Result, QisError> { + for &position in qubits { + if position >= self.state.num_qubits { + return Err(QisError::IndexOutOfBounds { + index: position, + max: self.state.num_qubits.saturating_sub(1), + }); + } + } let mut probs = self.probabilities(); self.apply_readout_noise( &mut probs, - qubits.iter().map(|&q| (Qubit::new(q as u32), q)), + qubits + .iter() + .map(|&position| (self.noise_qubit(position), position)), )?; Ok(probs) } diff --git a/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs b/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs index 82775269..3c8a055d 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs @@ -19,6 +19,37 @@ use ndarray::array; use num_complex::Complex64; use std::f64::consts::PI; +#[test] +fn duplicate_targets_leave_noisy_density_matrix_unchanged() { + let mut model = NoiseModel::new(); + model + .add_two_qubit_error( + StandardGate::CX, + Qubit::new(0), + Qubit::new(2), + TwoQubitNoise::CorrelatedPauli { + op_q0: Pauli::X, + op_q1: Pauli::X, + p: 0.2, + }, + ) + .unwrap(); + let mut initial = DensityMatrixNoise::new(3, Some(model)); + initial.apply_h(0).unwrap(); + initial.apply_s(0).unwrap(); + for gate in [StandardGate::CX, StandardGate::CZ] { + for q in 0..3 { + let mut sim = initial.clone(); + assert!(matches!( + sim.apply_standard_gate_noise(gate, &[q, q], &[]), + Err(QisError::InvalidParameterValue(msg)) + if msg.contains(&format!("Duplicate qubit index {q} ")) + )); + assert_eq!(sim.state.data(), initial.state.data()); + } + } +} + #[test] fn test_bit_flip_noise() { let mut noise_model = NoiseModel::new(); @@ -1176,6 +1207,146 @@ fn test_native_gates_amplitude_damping() { assert!(probs[0] > 0.0); } +#[test] +fn test_readout_uses_circuit_ids_at_selected_storage_positions() { + use crate::device::Outcome; + + for ids in [[0, 1, 2], [2, 0, 1], [2, 5, 9], [9, 2, 5]] { + for basis in 0..8 { + for (noise_position, &noise_id) in ids.iter().enumerate() { + let mut circuit = + Circuit::from_qubits(ids.iter().copied().map(Qubit::new).collect()).unwrap(); + for (position, &id) in ids.iter().enumerate() { + if basis & (1 << position) != 0 { + circuit.x(Qubit::new(id)).unwrap(); + } + } + let receipt = circuit.measure(Qubit::new(noise_id)).unwrap(); + let mut model = NoiseModel::new(); + model + .add_readout_error( + Qubit::new(noise_id), + ReadoutError { + p_0_given_1: 0.2, + p_1_given_0: 0.7, + }, + ) + .unwrap(); + let sim = DensityMatrixNoise::from_circuit(&circuit, Some(model)).unwrap(); + let original = sim.state.data().to_vec(); + let bit = (basis >> noise_position) & 1; + let flip_probability = if bit == 0 { 0.7 } else { 0.2 }; + + for selection in 0..8 { + let positions: Vec<_> = (0..3) + .rev() + .filter(|&position| selection & (1 << position) != 0) + .collect(); + let mut expected = [0.0; 8]; + expected[basis] = 1.0; + if selection & (1 << noise_position) != 0 { + expected[basis] = 1.0 - flip_probability; + expected[basis ^ (1 << noise_position)] = flip_probability; + } + let actual = sim.probabilities_with_readout(&positions).unwrap(); + assert_eq!(actual.len(), expected.len()); + for (actual, expected) in actual.iter().zip(&expected) { + assert!( + (actual - expected).abs() < 1e-12, + "IDs={ids:?}, basis={basis}, noise_position={noise_position}, positions={positions:?}: {actual} != {expected}" + ); + } + } + + let marginal = sim.probs_with_readout(&receipt).unwrap(); + let p_one = if bit == 0 { 0.7 } else { 0.8 }; + assert!((marginal[&Outcome::from_bitstring("1").unwrap()] - p_one).abs() < 1e-12); + assert_eq!(sim.state.data(), original.as_slice()); + } + } + } +} + +#[test] +fn test_readout_combines_mapped_errors_without_reordering_full_probabilities() { + use crate::device::Outcome; + + let mut circuit = + Circuit::from_qubits(vec![Qubit::new(9), Qubit::new(2), Qubit::new(5)]).unwrap(); + circuit.x(Qubit::new(9)).unwrap(); + let receipt = circuit + .measure_bits([Qubit::new(2), Qubit::new(9)]) + .unwrap(); + let mut model = NoiseModel::new(); + for (id, p_0_given_1, p_1_given_0) in [(9, 0.1, 0.2), (2, 0.3, 0.4)] { + model + .add_readout_error( + Qubit::new(id), + ReadoutError { + p_0_given_1, + p_1_given_0, + }, + ) + .unwrap(); + } + let sim = DensityMatrixNoise::from_circuit(&circuit, Some(model)).unwrap(); + let expected = [0.06, 0.54, 0.04, 0.36, 0.0, 0.0, 0.0, 0.0]; + for positions in [&[0, 1][..], &[1, 0], &[2, 1, 0]] { + let actual = sim.probabilities_with_readout(positions).unwrap(); + assert_eq!(actual.len(), expected.len()); + for (actual, expected) in actual.iter().zip(expected) { + assert!((actual - expected).abs() < 1e-12); + } + } + let marginal = sim.probs_with_readout(&receipt).unwrap(); + for (bits, probability) in [("00", 0.06), ("10", 0.54), ("01", 0.04), ("11", 0.36)] { + assert!((marginal[&Outcome::from_bitstring(bits).unwrap()] - probability).abs() < 1e-12); + } + assert_eq!( + sim.probabilities(), + vec![0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + ); +} + +#[test] +fn test_readout_validates_positions_independently_of_noise_model() { + let mut model = NoiseModel::new(); + model + .add_readout_error( + Qubit::new(9), + ReadoutError { + p_0_given_1: 0.0, + p_1_given_0: 1.0, + }, + ) + .unwrap(); + for noise_model in [None, Some(NoiseModel::new()), Some(model)] { + for bound in [false, true] { + let mut sim = DensityMatrixNoise::new(2, noise_model.clone()); + if bound { + let circuit = Circuit::from_qubits(vec![Qubit::new(9), Qubit::new(2)]).unwrap(); + sim.apply_circuit(&circuit).unwrap(); + } + for index in [2, 9, usize::MAX] { + assert!(matches!( + sim.probabilities_with_readout(&[0, index]), + Err(QisError::IndexOutOfBounds { index: actual, max: 1 }) if actual == index + )); + } + assert_eq!( + sim.probabilities_with_readout(&[]).unwrap(), + sim.probabilities() + ); + } + let sim = DensityMatrixNoise::new(0, noise_model); + assert_eq!(sim.probabilities_with_readout(&[]).unwrap(), vec![1.0]); + assert!(matches!( + sim.probabilities_with_readout(&[0]), + Err(QisError::IndexOutOfBounds { index: 0, max: 0 }) + )); + } +} + #[test] fn test_apply_readout_noise_out_of_bounds() { let mut sim = DensityMatrixNoise::new(2, Some(NoiseModel::new())); @@ -1578,3 +1749,169 @@ fn identity_noise_respects_labels_and_repetition() { .is_err() ); } + +#[test] +fn controlled_gphase_circuit_evolution() { + for ids in [[0, 1, 2], [12, 5, 31]] { + let qubits = ids.map(Qubit::new); + for controls in [ + vec![qubits[2]], + vec![qubits[2], qubits[0]], + vec![qubits[0], qubits[2]], + ] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut preparation = Circuit::from_qubits(qubits.to_vec()).unwrap(); + for qubit in qubits { + preparation.h(qubit).unwrap(); + } + preparation.s(qubits[1]).unwrap(); + let mut phase = Circuit::from_qubits(qubits.to_vec()).unwrap(); + phase + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mut circuit = preparation.clone(); + circuit + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mask = controls.iter().fold(0usize, |mask, q| { + mask | (1 << qubits.iter().position(|candidate| candidate == q).unwrap()) + }); + let expected: Vec = (0..8) + .map(|index| { + let amplitude = if index & 2 == 0 { + Complex64::new(1.0, 0.0) + } else { + Complex64::new(0.0, 1.0) + } / 8.0_f64.sqrt(); + if index & mask == mask { + amplitude * Complex64::from_polar(1.0, theta) + } else { + amplitude + } + }) + .collect(); + let expected_data = expected + .iter() + .flat_map(|a| expected.iter().map(move |b| a * b.conj())) + .collect::>(); + let actual = DensityMatrixNoise::from_circuit(&circuit, None).unwrap(); + let mut incremental = DensityMatrixNoise::new(3, None); + incremental.apply_circuit(&preparation).unwrap(); + incremental.apply_circuit(&phase).unwrap(); + for values in [actual.state.data(), incremental.state.data()] { + for (actual, expected) in values.iter().zip(&expected_data) { + assert!( + (actual - expected).norm() < 1e-12, + "{actual:?} != {expected:?}" + ); + } + } + let mut inverse = Circuit::from_qubits(qubits.to_vec()).unwrap(); + inverse + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [(-theta).into()], + ) + .unwrap(); + incremental.apply_circuit(&inverse).unwrap(); + let baseline = DensityMatrixNoise::from_circuit(&preparation, None).unwrap(); + for (actual, expected) in incremental.state.data().iter().zip(baseline.state.data()) + { + assert!((actual - expected).norm() < 1e-12); + } + } + } + } +} + +#[test] +fn controlled_gphase_noise_applies_after_unitary_once() { + let qubits = [Qubit::new(12), Qubit::new(5)]; + let mut circuit = Circuit::from_qubits(qubits.to_vec()).unwrap(); + circuit.x(qubits[0]).unwrap(); + circuit.h(qubits[1]).unwrap(); + circuit + .multi_control( + StandardGate::GPhase, + [qubits[1]], + Vec::::new(), + [0.37.into()], + ) + .unwrap(); + let mut model = NoiseModel::new(); + model + .add_single_qubit_error( + StandardGate::GPhase, + qubits[1], + SingleQubitNoise::BitFlip(0.23), + ) + .unwrap(); + model + .add_single_qubit_error( + StandardGate::GPhase, + qubits[0], + SingleQubitNoise::BitFlip(1.0), + ) + .unwrap(); + let sim = DensityMatrixNoise::from_circuit(&circuit, Some(model)).unwrap(); + let amplitudes = [ + Complex64::new(0.0, 0.0), + Complex64::new(1.0, 0.0) / 2.0_f64.sqrt(), + Complex64::new(0.0, 0.0), + Complex64::from_polar(1.0 / 2.0_f64.sqrt(), 0.37), + ]; + for row in 0..4 { + for col in 0..4 { + let expected = 0.77 * amplitudes[row] * amplitudes[col].conj() + + 0.23 * amplitudes[row ^ 2] * amplitudes[col ^ 2].conj(); + assert!((sim.state.data()[row * 4 + col] - expected).norm() < 1e-12); + } + } +} + +#[test] +fn zero_control_gphase_circuit_evolution() { + use crate::circuit::gate::{Instruction, MCGate}; + for num_qubits in [0, 1, 3] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut circuit = Circuit::new(num_qubits); + if num_qubits > 0 { + circuit.h(Qubit::new(0)).unwrap(); + } + circuit + .append( + Instruction::McGate(Box::new(MCGate::new(0, StandardGate::GPhase))), + Vec::::new(), + [theta.into()], + None, + ) + .unwrap(); + let actual = DensityMatrixNoise::from_circuit(&circuit, None).unwrap(); + let mut amplitudes = vec![Complex64::new(0.0, 0.0); 1 << num_qubits]; + amplitudes[0] = Complex64::from_polar(1.0, theta); + if num_qubits > 0 { + amplitudes[0] /= 2.0_f64.sqrt(); + amplitudes[1] = amplitudes[0]; + } + let expected = amplitudes + .iter() + .flat_map(|a| amplitudes.iter().map(move |b| a * b.conj())) + .collect::>(); + for (actual, expected) in actual.state.data().iter().zip(&expected) { + assert!((actual - expected).norm() < 1e-12); + } + } + } +} diff --git a/crates/cqlib-core/src/qis/state/density_matrix_test.rs b/crates/cqlib-core/src/qis/state/density_matrix_test.rs index 91595c4b..95e89ce8 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_test.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_test.rs @@ -408,6 +408,98 @@ fn test_zeros_and_add_assign() { assert_relative_eq!(dm1.probabilities()[1], 0.5); } +#[test] +fn duplicate_two_qubit_targets_leave_density_matrix_unchanged() { + let mut initial = DensityMatrix::new(3); + initial.apply_h(0).unwrap(); + initial.apply_ry(1, 0.37).unwrap(); + initial.apply_s(0).unwrap(); + for gate in [ + StandardGate::CX, + StandardGate::CZ, + StandardGate::CY, + StandardGate::SWAP, + StandardGate::RXX, + StandardGate::RYY, + StandardGate::RZZ, + StandardGate::RZX, + StandardGate::CRX, + StandardGate::CRY, + StandardGate::CRZ, + StandardGate::FSIM, + ] { + for q in 0..3 { + let mut dm = initial.clone(); + let result = dm.apply_standard_gate(gate, &[q, q], &vec![0.37; gate.num_params()]); + assert!(matches!(result, Err(QisError::InvalidParameterValue(msg)) + if msg.contains(&format!("Duplicate qubit index {q} ")))); + assert_eq!(dm.data(), initial.data(), "{gate:?} on [{q}, {q}]"); + } + } + for targets in [[0, 0, 2], [0, 2, 0], [2, 0, 0]] { + let mut dm = initial.clone(); + assert!(matches!( + dm.apply_ccx(targets[0], targets[1], targets[2]), + Err(QisError::InvalidParameterValue(msg)) if msg.contains("Duplicate qubit index 0 ") + )); + assert_eq!(dm.data(), initial.data()); + } +} + +#[test] +fn duplicate_custom_gate_targets_leave_density_matrix_unchanged() { + let mut dm = DensityMatrix::new(3); + dm.apply_h(0).unwrap(); + dm.apply_s(0).unwrap(); + let before = dm.data().to_vec(); + let matrix = StandardGate::CX.matrix(&[]).unwrap(); + let matrix = std::array::from_fn(|r| std::array::from_fn(|c| matrix[[r, c]])); + for q in 0..3 { + assert!(matches!( + dm.apply_two_qubit_gate(q, q, matrix), + Err(QisError::InvalidParameterValue(msg)) + if msg.contains(&format!("Duplicate qubit index {q} ")) + )); + assert_eq!(dm.data(), before); + } +} + +#[test] +fn duplicate_kraus_targets_leave_density_matrix_unchanged() { + let amplitudes: Vec<_> = (1..=8) + .map(|i| Complex64::new(i as f64 / 204.0_f64.sqrt(), 0.0)) + .collect(); + let initial = DensityMatrix::from_state(3, amplitudes).unwrap(); + for (targets, duplicate) in [ + (&[0, 0][..], 0), + (&[1, 1][..], 1), + (&[2, 2][..], 2), + (&[0, 2, 0][..], 0), + (&[0, 0, 2][..], 0), + (&[2, 0, 0][..], 0), + ] { + let dim = 1 << targets.len(); + let norm: f64 = (1..=dim).map(|i| (i * i) as f64).sum(); + let operator = (0..dim) + .flat_map(|r| { + (0..dim).map(move |c| { + Complex64::new( + if r == c { 1.0 } else { 0.0 } - 2.0 * ((r + 1) * (c + 1)) as f64 / norm, + 0.0, + ) + }) + }) + .collect(); + let mut dm = initial.clone(); + assert!(matches!( + dm.apply_kraus(&[operator], targets), + Err(QisError::InvalidParameterValue(msg)) + if msg.contains(&format!("Duplicate qubit index {duplicate} ")) + )); + assert_eq!(dm.data(), initial.data()); + } +} + #[test] fn test_single_qubit_gates() { let mut dm = DensityMatrix::new(1); @@ -1240,3 +1332,123 @@ fn hermiticity_rejects_nonfinite_entries() { } } } + +#[test] +fn controlled_gphase_circuit_evolution() { + for ids in [[0, 1, 2], [12, 5, 31]] { + let qubits = ids.map(Qubit::new); + for controls in [ + vec![qubits[2]], + vec![qubits[2], qubits[0]], + vec![qubits[0], qubits[2]], + ] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut preparation = Circuit::from_qubits(qubits.to_vec()).unwrap(); + for qubit in qubits { + preparation.h(qubit).unwrap(); + } + preparation.s(qubits[1]).unwrap(); + let mut phase = Circuit::from_qubits(qubits.to_vec()).unwrap(); + phase + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mut circuit = preparation.clone(); + circuit + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mask = controls.iter().fold(0usize, |mask, q| { + mask | (1 << qubits.iter().position(|candidate| candidate == q).unwrap()) + }); + let expected: Vec = (0..8) + .map(|index| { + let amplitude = if index & 2 == 0 { + Complex64::new(1.0, 0.0) + } else { + Complex64::new(0.0, 1.0) + } / 8.0_f64.sqrt(); + if index & mask == mask { + amplitude * Complex64::from_polar(1.0, theta) + } else { + amplitude + } + }) + .collect(); + let expected_data = expected + .iter() + .flat_map(|a| expected.iter().map(move |b| a * b.conj())) + .collect::>(); + let actual = DensityMatrix::from_circuit(&circuit).unwrap(); + let mut incremental = DensityMatrix::new(3); + incremental.apply_circuit(&preparation).unwrap(); + incremental.apply_circuit(&phase).unwrap(); + for values in [actual.data(), incremental.data()] { + for (actual, expected) in values.iter().zip(&expected_data) { + assert!( + (actual - expected).norm() < 1e-12, + "{actual:?} != {expected:?}" + ); + } + } + let mut inverse = Circuit::from_qubits(qubits.to_vec()).unwrap(); + inverse + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [(-theta).into()], + ) + .unwrap(); + incremental.apply_circuit(&inverse).unwrap(); + let baseline = DensityMatrix::from_circuit(&preparation).unwrap(); + for (actual, expected) in incremental.data().iter().zip(baseline.data()) { + assert!((actual - expected).norm() < 1e-12); + } + } + } + } +} + +#[test] +fn zero_control_gphase_circuit_evolution() { + use crate::circuit::gate::{Instruction, MCGate}; + for num_qubits in [0, 1, 3] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut circuit = Circuit::new(num_qubits); + if num_qubits > 0 { + circuit.h(Qubit::new(0)).unwrap(); + } + circuit + .append( + Instruction::McGate(Box::new(MCGate::new(0, StandardGate::GPhase))), + Vec::::new(), + [theta.into()], + None, + ) + .unwrap(); + let actual = DensityMatrix::from_circuit(&circuit).unwrap(); + let mut amplitudes = vec![Complex64::new(0.0, 0.0); 1 << num_qubits]; + amplitudes[0] = Complex64::from_polar(1.0, theta); + if num_qubits > 0 { + amplitudes[0] /= 2.0_f64.sqrt(); + amplitudes[1] = amplitudes[0]; + } + let expected = amplitudes + .iter() + .flat_map(|a| amplitudes.iter().map(move |b| a * b.conj())) + .collect::>(); + for (actual, expected) in actual.data().iter().zip(&expected) { + assert!((actual - expected).norm() < 1e-12); + } + } + } +} diff --git a/crates/cqlib-core/src/qis/state/stabilizer.rs b/crates/cqlib-core/src/qis/state/stabilizer.rs index c62ec501..99d4640b 100644 --- a/crates/cqlib-core/src/qis/state/stabilizer.rs +++ b/crates/cqlib-core/src/qis/state/stabilizer.rs @@ -1662,20 +1662,16 @@ impl StabilizerState { Instruction::Standard(gate) => { state.apply_clifford_gate(*gate, &qubit_indices, ¶ms)?; } - Instruction::McGate(mc) if mc.num_ctrl_qubits() == 1 => { - let ctrl = qubit_indices[0]; - let tgt = qubit_indices[1]; - match mc.base_gate() { - StandardGate::X => state.apply_cx(ctrl, tgt)?, - StandardGate::Y => state.apply_cy(ctrl, tgt)?, - StandardGate::Z => state.apply_cz(ctrl, tgt)?, - g => { - return Err(QisError::UnsupportedOperation(format!( - "Controlled-{g} is not a Clifford gate" - ))); - } + Instruction::McGate(mc) if mc.num_ctrl_qubits() == 1 => match mc.base_gate() { + StandardGate::X => state.apply_cx(qubit_indices[0], qubit_indices[1])?, + StandardGate::Y => state.apply_cy(qubit_indices[0], qubit_indices[1])?, + StandardGate::Z => state.apply_cz(qubit_indices[0], qubit_indices[1])?, + g => { + return Err(QisError::UnsupportedOperation(format!( + "Controlled-{g} is not supported in stabilizer simulation" + ))); } - } + }, Instruction::McGate(mc) => { return Err(QisError::UnsupportedOperation(format!( "{}-controlled gates are not supported in stabilizer simulation", diff --git a/crates/cqlib-core/src/qis/state/stabilizer_test.rs b/crates/cqlib-core/src/qis/state/stabilizer_test.rs index 2633a1d2..8876a098 100644 --- a/crates/cqlib-core/src/qis/state/stabilizer_test.rs +++ b/crates/cqlib-core/src/qis/state/stabilizer_test.rs @@ -1434,3 +1434,58 @@ fn test_apply_circuit_rejects_classical_control_flow() { let error = StabilizerState::run_circuit(&c).unwrap_err(); assert!(matches!(error, QisError::UnsupportedOperation(_))); } + +#[test] +fn controlled_gphase_returns_unsupported_without_panicking() { + use crate::circuit::{Circuit, StandardGate}; + for theta in [0.0, 0.37, std::f64::consts::FRAC_PI_2] { + let mut circuit = Circuit::from_qubits(vec![Qubit::new(12)]).unwrap(); + circuit + .multi_control( + StandardGate::GPhase, + [Qubit::new(12)], + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mut state = StabilizerState::new(1); + for result in [ + StabilizerState::from_circuit(&circuit).map(|_| ()), + state.apply_circuit(&circuit), + StabilizerState::run_circuit(&circuit).map(|_| ()), + ] { + assert!( + matches!(result, Err(QisError::UnsupportedOperation(ref message)) + if message.contains("Controlled-GPhase") && message.contains("not supported")), + "unexpected result: {result:?}" + ); + } + } +} + +#[test] +fn controlled_paulis_match_standard_stabilizer_execution() { + use crate::circuit::{Circuit, StandardGate}; + for gate in [StandardGate::X, StandardGate::Y, StandardGate::Z] { + let mut circuit = Circuit::new(2); + circuit.h(Qubit::new(0)).unwrap(); + circuit.h(Qubit::new(1)).unwrap(); + let mut expected = StabilizerState::from_circuit(&circuit).unwrap(); + match gate { + StandardGate::X => expected.apply_cx(0, 1).unwrap(), + StandardGate::Y => expected.apply_cy(0, 1).unwrap(), + StandardGate::Z => expected.apply_cz(0, 1).unwrap(), + _ => unreachable!(), + } + circuit + .multi_control(gate, [Qubit::new(0)], [Qubit::new(1)], []) + .unwrap(); + for actual in [ + StabilizerState::from_circuit(&circuit).unwrap(), + StabilizerState::run_circuit(&circuit).unwrap().state, + ] { + assert_eq!(&*actual.tableau, &*expected.tableau); + assert_eq!(actual.phases, expected.phases); + } + } +} diff --git a/crates/cqlib-core/src/qis/state/statevector.rs b/crates/cqlib-core/src/qis/state/statevector.rs index 9b4f85da..21539d5e 100644 --- a/crates/cqlib-core/src/qis/state/statevector.rs +++ b/crates/cqlib-core/src/qis/state/statevector.rs @@ -576,44 +576,31 @@ impl Statevector { let num_controls = mc_gate.num_ctrl_qubits(); let base_gate = mc_gate.base_gate(); - if num_controls == 1 { - // Single-controlled gate - let control = qubit_indices[0]; - let target = qubit_indices[1]; - match base_gate { - StandardGate::X => sv.apply_cx(control, target)?, - StandardGate::Y => sv.apply_cy(control, target)?, - StandardGate::Z => sv.apply_cz(control, target)?, - StandardGate::RX => { - let theta = params[0]; - sv.apply_crx(control, target, theta)?; - } - StandardGate::RY => { - let theta = params[0]; - sv.apply_cry(control, target, theta)?; - } - StandardGate::RZ => { - let theta = params[0]; - sv.apply_crz(control, target, theta)?; - } - _ => { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError(CircuitError::NoMatrixRepresentation) - })?; - sv.apply_unitary_gate(&qubit_indices, &matrix)?; - } + match (num_controls, *base_gate) { + (0, StandardGate::GPhase) => { + sv.apply_standard_gate(*base_gate, &qubit_indices, ¶ms)? + } + (1, StandardGate::X) => sv.apply_cx(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::Y) => sv.apply_cy(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::Z) => sv.apply_cz(qubit_indices[0], qubit_indices[1])?, + (1, StandardGate::RX) => { + sv.apply_crx(qubit_indices[0], qubit_indices[1], params[0])? + } + (1, StandardGate::RY) => { + sv.apply_cry(qubit_indices[0], qubit_indices[1], params[0])? + } + (1, StandardGate::RZ) => { + sv.apply_crz(qubit_indices[0], qubit_indices[1], params[0])? + } + (2, StandardGate::X) => { + sv.apply_ccx(qubit_indices[0], qubit_indices[1], qubit_indices[2])? + } + _ => { + let matrix = mc_gate.matrix(¶ms).map_err(|_| { + QisError::CircuitError(CircuitError::NoMatrixRepresentation) + })?; + sv.apply_unitary_gate(&qubit_indices, &matrix)?; } - } else if num_controls == 2 && *base_gate == StandardGate::X { - // Toffoli gate - let c0 = qubit_indices[0]; - let c1 = qubit_indices[1]; - let target = qubit_indices[2]; - sv.apply_ccx(c0, c1, target)?; - } else { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError(CircuitError::NoMatrixRepresentation) - })?; - sv.apply_unitary_gate(&qubit_indices, &matrix)?; } } Instruction::UnitaryGate(u_gate) => { diff --git a/crates/cqlib-core/src/qis/state/statevector_test.rs b/crates/cqlib-core/src/qis/state/statevector_test.rs index 28f9d8a3..83fc75d1 100644 --- a/crates/cqlib-core/src/qis/state/statevector_test.rs +++ b/crates/cqlib-core/src/qis/state/statevector_test.rs @@ -2281,3 +2281,117 @@ fn test_apply_pauli_rotation_rejects_imaginary_global_phase() { pauli.phase = Phase::I; assert!(state.apply_pauli_rotation(&pauli, 0.2).is_err()); } + +#[test] +fn controlled_gphase_circuit_evolution() { + for ids in [[0, 1, 2], [12, 5, 31]] { + let qubits = ids.map(Qubit::new); + for controls in [ + vec![qubits[2]], + vec![qubits[2], qubits[0]], + vec![qubits[0], qubits[2]], + ] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut preparation = Circuit::from_qubits(qubits.to_vec()).unwrap(); + for qubit in qubits { + preparation.h(qubit).unwrap(); + } + preparation.s(qubits[1]).unwrap(); + let mut phase = Circuit::from_qubits(qubits.to_vec()).unwrap(); + phase + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mut circuit = preparation.clone(); + circuit + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [theta.into()], + ) + .unwrap(); + let mask = controls.iter().fold(0usize, |mask, q| { + mask | (1 << qubits.iter().position(|candidate| candidate == q).unwrap()) + }); + let expected: Vec = (0..8) + .map(|index| { + let amplitude = if index & 2 == 0 { + Complex64::new(1.0, 0.0) + } else { + Complex64::new(0.0, 1.0) + } / 8.0_f64.sqrt(); + if index & mask == mask { + amplitude * Complex64::from_polar(1.0, theta) + } else { + amplitude + } + }) + .collect(); + let expected_data = expected.clone(); + let actual = Statevector::from_circuit(&circuit).unwrap(); + let mut incremental = Statevector::new(3); + incremental.apply_circuit(&preparation).unwrap(); + incremental.apply_circuit(&phase).unwrap(); + for values in [actual.data(), incremental.data()] { + for (actual, expected) in values.iter().zip(&expected_data) { + assert!( + (actual - expected).norm() < 1e-12, + "{actual:?} != {expected:?}" + ); + } + } + let mut inverse = Circuit::from_qubits(qubits.to_vec()).unwrap(); + inverse + .multi_control( + StandardGate::GPhase, + controls.clone(), + Vec::::new(), + [(-theta).into()], + ) + .unwrap(); + incremental.apply_circuit(&inverse).unwrap(); + let baseline = Statevector::from_circuit(&preparation).unwrap(); + for (actual, expected) in incremental.data().iter().zip(baseline.data()) { + assert!((actual - expected).norm() < 1e-12); + } + } + } + } +} + +#[test] +fn zero_control_gphase_circuit_evolution() { + use crate::circuit::gate::{Instruction, MCGate}; + for num_qubits in [0, 1, 3] { + for theta in [0.0, 0.37, -0.51, PI] { + let mut circuit = Circuit::new(num_qubits); + if num_qubits > 0 { + circuit.h(Qubit::new(0)).unwrap(); + } + circuit + .append( + Instruction::McGate(Box::new(MCGate::new(0, StandardGate::GPhase))), + Vec::::new(), + [theta.into()], + None, + ) + .unwrap(); + let actual = Statevector::from_circuit(&circuit).unwrap(); + let mut amplitudes = vec![Complex64::new(0.0, 0.0); 1 << num_qubits]; + amplitudes[0] = Complex64::from_polar(1.0, theta); + if num_qubits > 0 { + amplitudes[0] /= 2.0_f64.sqrt(); + amplitudes[1] = amplitudes[0]; + } + let expected = amplitudes; + for (actual, expected) in actual.data().iter().zip(&expected) { + assert!((actual - expected).norm() < 1e-12); + } + } + } +} diff --git a/crates/cqlib-core/src/visualization/state/math.rs b/crates/cqlib-core/src/visualization/state/math.rs index 490fb7a8..c65890da 100644 --- a/crates/cqlib-core/src/visualization/state/math.rs +++ b/crates/cqlib-core/src/visualization/state/math.rs @@ -84,6 +84,8 @@ pub fn local_bloch_vectors( /// The returned labels enumerate tensor products over `I`, `X`, `Y`, and `Z`. /// Values are expectation coefficients `Tr(rho P)` in the same qubit indexing convention /// used by the statevector basis: qubit `q` corresponds to bit `1 << q`. +/// Labels follow [`crate::qis::pauli::PauliString`]: the highest-index qubit is on +/// the left by default. `reverse_bits` reverses the displayed characters only. pub(crate) fn pauli_coefficients( num_qubits: usize, rho: &[Complex64], @@ -97,6 +99,11 @@ pub(crate) fn pauli_coefficients( ))); } + let display_qubits: Vec = if reverse_bits { + (0..num_qubits).collect() + } else { + (0..num_qubits).rev().collect() + }; let mut out = Vec::with_capacity(total); for idx in 0..total { let mut code = idx; @@ -111,11 +118,6 @@ pub(crate) fn pauli_coefficients( ops.push((qubit, op)); code /= 4; } - let display_qubits: Vec = if reverse_bits { - (0..num_qubits).rev().collect() - } else { - (0..num_qubits).collect() - }; let label = display_qubits .iter() .map(|qubit| match ops[*qubit].1 { diff --git a/crates/cqlib-core/src/visualization/state/math_test.rs b/crates/cqlib-core/src/visualization/state/math_test.rs index 318a584b..1b14e9ad 100644 --- a/crates/cqlib-core/src/visualization/state/math_test.rs +++ b/crates/cqlib-core/src/visualization/state/math_test.rs @@ -14,7 +14,8 @@ use super::math::pauli_coefficients; use super::*; -use crate::qis::Statevector; +use crate::qis::pauli::PauliString; +use crate::qis::{DensityMatrix, Statevector}; use num_complex::Complex64; fn plus_state() -> Statevector { @@ -86,3 +87,81 @@ fn paulivec_contains_expected_plus_y() { .any(|(label, value)| label == "Y" && (*value - 1.0).abs() < 1e-10) ); } + +#[test] +fn paulivec_labels_identify_the_correct_qubit() { + // q0 = |+>, q1 = |0>: IX measures q0, XI measures q1. + let state = Statevector::from_state( + 2, + vec![ + Complex64::new(1.0 / 2.0_f64.sqrt(), 0.0), + Complex64::new(1.0 / 2.0_f64.sqrt(), 0.0), + Complex64::new(0.0, 0.0), + Complex64::new(0.0, 0.0), + ], + ) + .unwrap(); + let (_, rho) = state_to_density_matrix(&state).unwrap(); + for reverse_bits in [false, true] { + let coeffs = pauli_coefficients(2, &rho, reverse_bits).unwrap(); + for (label, expected) in [("IX", 1.0), ("XI", 0.0), ("IZ", 0.0), ("ZI", 1.0)] { + let displayed = if reverse_bits { + label.chars().rev().collect::() + } else { + label.to_string() + }; + let value = coeffs.iter().find(|(s, _)| *s == displayed).unwrap().1; + assert!((value - expected).abs() < 1e-10, "{displayed}: {value}"); + } + } +} + +fn assert_paulivec_matches_observables( + state: &impl StateVisualizationSource, + expectation: impl Fn(&PauliString) -> f64, +) { + let (num_qubits, rho) = state_to_density_matrix(state).unwrap(); + for reverse_bits in [false, true] { + let coeffs = pauli_coefficients(num_qubits, &rho, reverse_bits).unwrap(); + assert_eq!(coeffs.len(), 4usize.pow(num_qubits as u32)); + for (label, value) in coeffs { + let semantic_label = if reverse_bits { + label.chars().rev().collect::() + } else { + label.clone() + }; + let pauli: PauliString = semantic_label.parse().unwrap(); + let expected = expectation(&pauli); + assert!( + (value - expected).abs() < 1e-10, + "{num_qubits} qubits, reverse_bits={reverse_bits}, {label}: {value} != {expected}" + ); + } + } +} + +#[test] +fn paulivec_coefficients_match_observables_for_pure_and_mixed_states() { + for num_qubits in 1..=3 { + let dim = 1 << num_qubits; + // Asymmetric complex amplitudes exercise X, Y, Z and negative coefficients. + let mut amplitudes: Vec<_> = (0..dim) + .map(|i| Complex64::new(((i + 1) as f64 * 0.71).sin(), ((i + 2) as f64 * 0.43).cos())) + .collect(); + let norm = amplitudes.iter().map(|a| a.norm_sqr()).sum::().sqrt(); + for amplitude in &mut amplitudes { + *amplitude /= norm; + } + let pure = Statevector::from_state(num_qubits, amplitudes).unwrap(); + assert_paulivec_matches_observables(&pure, |p| pure.expectation(p).unwrap()); + + // A convex mixture with |11...1><11...1| is a genuine mixed state. + let (_, mut rho) = state_to_density_matrix(&pure).unwrap(); + for value in &mut rho { + *value *= 0.7; + } + rho[dim * dim - 1] += Complex64::new(0.3, 0.0); + let mixed = DensityMatrix::from_density_matrix_state(num_qubits, rho).unwrap(); + assert_paulivec_matches_observables(&mixed, |p| mixed.expectation(p).unwrap()); + } +} diff --git a/crates/cqlib-core/src/visualization/state/plot.rs b/crates/cqlib-core/src/visualization/state/plot.rs index 2243f3ed..cadd93d4 100644 --- a/crates/cqlib-core/src/visualization/state/plot.rs +++ b/crates/cqlib-core/src/visualization/state/plot.rs @@ -195,6 +195,11 @@ pub fn plot_state_city( /// Plot Pauli-basis expectation values as an SVG bar chart. /// +/// By default, labels follow [`crate::qis::pauli::PauliString`]: the rightmost +/// character acts on qubit 0. For example, `IX` measures X on qubit 0 of a +/// two-qubit state. Setting `options.reverse_bits` reverses the characters in +/// each label without changing the state or the corresponding expectation value. +/// /// # Arguments /// /// * `state` - Input core state object accepted by state visualization routines. diff --git a/crates/cqlib-core/src/visualization/state/plot_test.rs b/crates/cqlib-core/src/visualization/state/plot_test.rs index 256ebb04..25e286a6 100644 --- a/crates/cqlib-core/src/visualization/state/plot_test.rs +++ b/crates/cqlib-core/src/visualization/state/plot_test.rs @@ -172,6 +172,54 @@ fn paulivec_matches_reference_image() { assert_state_visual_match(&svg, "paulivec_plus.png"); } +#[test] +fn paulivec_svg_labels_match_asymmetric_state() { + let state = Statevector::from_state( + 2, + vec![ + Complex64::new(1.0 / 2.0_f64.sqrt(), 0.0), + Complex64::new(1.0 / 2.0_f64.sqrt(), 0.0), + Complex64::new(0.0, 0.0), + Complex64::new(0.0, 0.0), + ], + ) + .unwrap(); + let density = DensityMatrix::from_state(2, state.data().to_vec()).unwrap(); + for reverse_bits in [false, true] { + let options = StatePlotOptions { + reverse_bits, + ..StatePlotOptions::default() + }; + for svg in [ + plot_state_paulivec(&state, &options).unwrap(), + plot_state_paulivec(&density, &options).unwrap(), + ] { + let bar = |label: &str| { + let label_start = svg.find(&format!(">{label}")).unwrap(); + let bar_start = svg[..label_start].rfind(" 1.0); + for (label, expected) in [("IX", 1.0), ("XI", 0.0), ("IZ", 0.0), ("ZI", 1.0)] { + let displayed = if reverse_bits { + label.chars().rev().collect::() + } else { + label.to_string() + }; + let (y, height) = bar(&displayed); + if expected == 1.0 { + assert_eq!((y, height), (one_y, one_height), "{displayed}"); + } else { + assert_eq!(height, 1.0, "{displayed}"); + assert!((y - one_y - one_height).abs() < 1e-3, "{displayed}"); + } + } + } + } +} + #[test] fn density_matrix_uses_same_state_plot_api() { let statevector = plus_state(); diff --git a/crates/cqlib/Cargo.toml b/crates/cqlib/Cargo.toml index 650a1265..fbcc267d 100644 --- a/crates/cqlib/Cargo.toml +++ b/crates/cqlib/Cargo.toml @@ -26,4 +26,4 @@ categories = [ [dependencies] -cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.5" } +cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.6" } diff --git a/docs/api/cn/c/3_qis/3_density_matrix_noise.md b/docs/api/cn/c/3_qis/3_density_matrix_noise.md index 843d85fc..1f5ade40 100644 --- a/docs/api/cn/c/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/c/3_qis/3_density_matrix_noise.md @@ -250,7 +250,7 @@ int32_t density_matrix_noise_probabilities_with_readout(const struct CDensityMat - `ptr` (`const CDensityMatrixNoise*`)：噪声模拟器句柄。 - `buffer` (`double*`)：输出缓冲。 - `len` (`uintptr_t`)：缓冲长度，必须等于 `density_matrix_noise_probabilities_with_readout_len`。 -- `qubits` (`const uint32_t*`)：施加读出误差的比特数组；`qubit_len` 为 0 时可为 NULL，此时结果与 `density_matrix_noise_probabilities` 一致。 +- `qubits` (`const uint32_t*`)：施加读出误差的存储位置数组；绑定线路后按位置对应的线路比特 ID 查询噪声模型，未绑定时 ID 与位置相同。输出始终为按存储位序排列的完整概率数组。`qubit_len` 为 0 时可为 NULL，此时结果与 `density_matrix_noise_probabilities` 一致。 - `qubit_len` (`uintptr_t`)：比特数量。 返回：`0` 成功；`-1` NULL 指针；`-2` 比特越界；`-8` `len` 与实际数量不符。 diff --git a/docs/api/cn/c/5_visualization/3_state_plots.md b/docs/api/cn/c/5_visualization/3_state_plots.md index 2cef5286..61851474 100644 --- a/docs/api/cn/c/5_visualization/3_state_plots.md +++ b/docs/api/cn/c/5_visualization/3_state_plots.md @@ -76,6 +76,8 @@ char *plot_state_city_density_matrix(const struct CDensityMatrix *state, 把泡利基期望值绘制为柱状图。标签为 `I`、`X`、`Y`、`Z` 的张量积，共 `4^n` 条，其中 `n` 为量子比特数；非负柱与负值柱使用不同颜色。柱数不超过 `64` 时才绘制旋转的横轴标签。 +默认标签遵循 `PauliString`：最高位在左，qubit 0 在最右侧。例如，两比特态的 `IX` 表示 qubit 0 上的 X。启用 `reverse_bits` 仅反转每个标签中的字符顺序，对应的期望值不变。 + ```c char *plot_state_paulivec(const struct CStatevector *state, const struct CStatePlotOptions *options); diff --git a/docs/api/cn/python/3_qis/3_density_matrix_noise.md b/docs/api/cn/python/3_qis/3_density_matrix_noise.md index 07cd7924..331ac71c 100644 --- a/docs/api/cn/python/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/python/3_qis/3_density_matrix_noise.md @@ -171,7 +171,7 @@ assert sim.state[0, 0] != 999.0 **概率** - `probabilities() -> list[float]`：返回不含读出噪声的各基态概率，长度为 2^num_qubits。不修改状态。 -- `probabilities_with_readout(qubits) -> list[float]`：返回在指定比特上施加读出噪声后的概率分布，长度为 2^n，n 为 `qubits` 长度。 +- `probabilities_with_readout(qubits) -> list[float]`：`qubits` 为施加读出误差的存储位置；绑定线路后，按这些位置对应的线路比特 ID 查询噪声模型，未绑定时 ID 与位置相同。返回完整的 `2^num_qubits` 项概率数组，保持存储位序，不随选择的顺序或数量改变；空选择不施加读出误差。位置越界时抛出 `IndexError`，即使未配置噪声模型。 - `qubits` (`list[int]`)：参与读出噪声建模的比特下标列表。下标越界时抛出 `IndexError`。 **测量与采样** diff --git a/docs/api/cn/python/5_visualization/3_state_plots.md b/docs/api/cn/python/5_visualization/3_state_plots.md index 2bff0705..d563bd04 100644 --- a/docs/api/cn/python/5_visualization/3_state_plots.md +++ b/docs/api/cn/python/5_visualization/3_state_plots.md @@ -120,6 +120,8 @@ assert plot_state_paulivec(state) == plot_state_paulivec(density) 计算量子态在 Pauli 基下的期望值，并渲染成柱状图。 +默认标签遵循 `PauliString`：最高位在左，qubit 0 在最右侧。例如，两比特态的 `IX` 表示 qubit 0 上的 X。启用 `reverse_bits` 仅反转每个标签中的字符顺序，对应的期望值不变。 + 参数： - `state` (`QuantumState`)：待渲染的量子态。 diff --git a/docs/api/cn/rust/3_qis/3_density_matrix_noise.md b/docs/api/cn/rust/3_qis/3_density_matrix_noise.md index 47cc8b4b..bfab1bf8 100644 --- a/docs/api/cn/rust/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/rust/3_qis/3_density_matrix_noise.md @@ -87,7 +87,7 @@ pub struct DensityMatrixNoise { 模拟器把「量子态的噪声」与「读出的噪声」分开表达：前者作用于密度矩阵本身，后者只改变报告出来的概率分布。 - `fn probabilities(&self) -> Vec`：返回密度矩阵对角元，即各计算基态的测量概率。包含已经施加的门噪声，不含读出噪声。 -- `fn probabilities_with_readout(&self, qubits: &[usize]) -> Result, QisError>`：在 `probabilities` 的基础上，对 `qubits` 中列出的比特施加读出噪声。 +- `fn probabilities_with_readout(&self, qubits: &[usize]) -> Result, QisError>`：`qubits` 为施加读出误差的存储位置；绑定线路后，按这些位置对应的线路比特 ID 查询噪声模型，未绑定时 ID 与位置相同。返回完整的 `2^num_qubits` 项概率数组，保持存储位序；空选择不施加读出误差。位置越界时返回 `QisError::IndexOutOfBounds`，即使未配置噪声模型。 - `fn probs(&self, measurement: &Measurement) -> Result, QisError>`：返回 `Measurement` 选定的边缘概率分布，不含读出噪声。 - `fn probs_with_readout(&self, measurement: &Measurement) -> Result, QisError>`：在 `probs` 的基础上，对 `measurement.qubits()` 施加读出噪声。 diff --git a/docs/api/cn/rust/5_visualization/3_state_plots.md b/docs/api/cn/rust/5_visualization/3_state_plots.md index 1f35cca1..533209bf 100644 --- a/docs/api/cn/rust/5_visualization/3_state_plots.md +++ b/docs/api/cn/rust/5_visualization/3_state_plots.md @@ -118,6 +118,8 @@ assert!(svg.contains("Re[rho]")); 把泡利基期望值绘制为柱状图。标签为 `I`、`X`、`Y`、`Z` 的张量积，共 `4^n` 条，其中 `n` 为量子比特数。非负柱与负值柱使用不同颜色。 +默认标签遵循 `PauliString`：最高位在左，qubit 0 在最右侧。例如，两比特态的 `IX` 表示 qubit 0 上的 X。启用 `reverse_bits` 仅反转每个标签中的字符顺序，对应的期望值不变。 + 参数： - `state` (`&S`，`S: StateVisualizationSource + ?Sized`)：输入态。 diff --git a/docs/api/en/c/3_qis/3_density_matrix_noise.md b/docs/api/en/c/3_qis/3_density_matrix_noise.md index b542c22e..f6b24cc2 100644 --- a/docs/api/en/c/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/c/3_qis/3_density_matrix_noise.md @@ -250,7 +250,7 @@ Takes the ideal probability distribution and applies the readout errors register - `ptr` (`const CDensityMatrixNoise*`): noisy-simulator handle. - `buffer` (`double*`): output buffer. - `len` (`uintptr_t`): buffer length; must equal `density_matrix_noise_probabilities_with_readout_len`. -- `qubits` (`const uint32_t*`): array of qubits to apply readout errors to; may be NULL when `qubit_len` is 0, in which case the result matches `density_matrix_noise_probabilities`. +- `qubits` (`const uint32_t*`): storage positions on which to apply readout errors. Circuit-bound states look up noise using the circuit ID at each position; unbound states use identity IDs. The full output distribution retains storage order. May be NULL when `qubit_len` is 0, in which case the result matches `density_matrix_noise_probabilities`. - `qubit_len` (`uintptr_t`): number of qubits. Returns: `0` on success; `-1` NULL pointers; `-2` out-of-bounds qubit; `-8` when `len` does not match the actual count. diff --git a/docs/api/en/c/5_visualization/3_state_plots.md b/docs/api/en/c/5_visualization/3_state_plots.md index 8757ef43..fe296975 100644 --- a/docs/api/en/c/5_visualization/3_state_plots.md +++ b/docs/api/en/c/5_visualization/3_state_plots.md @@ -76,6 +76,8 @@ Failure scenarios: same as [plot_bloch_multivector](#plot_bloch_multivectorstate Plots Pauli-basis expectation values as a bar chart. Labels are tensor products of `I`, `X`, `Y`, `Z`, for `4^n` bars in total where `n` is the qubit count; non-negative and negative bars use different colors. Rotated x-axis labels are drawn only when the bar count does not exceed `64`. +By default, labels follow `PauliString`: the highest-index qubit is on the left and qubit 0 is on the right. For example, `IX` measures X on qubit 0 of a two-qubit state. Enabling `reverse_bits` reverses the characters in each label without changing the corresponding expectation value. + ```c char *plot_state_paulivec(const struct CStatevector *state, const struct CStatePlotOptions *options); diff --git a/docs/api/en/python/3_qis/3_density_matrix_noise.md b/docs/api/en/python/3_qis/3_density_matrix_noise.md index 1a318e58..0de4d8f6 100644 --- a/docs/api/en/python/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/python/3_qis/3_density_matrix_noise.md @@ -171,7 +171,7 @@ All gate application and state rewriting methods modify the current state in pla **Probabilities** - `probabilities() -> list[float]`: return the probability of each basis state without readout noise, with length 2^num_qubits. Does not modify the state. -- `probabilities_with_readout(qubits) -> list[float]`: return the probability distribution after readout noise is applied on the given qubits, with length 2^n, where n is the length of `qubits`. +- `probabilities_with_readout(qubits) -> list[float]`: `qubits` contains storage positions on which to apply readout errors. Circuit-bound states look up noise using the circuit ID at each position; unbound states use identity IDs. Returns all `2^num_qubits` probabilities in storage order, independent of the selection order or size; an empty selection applies no readout errors. Out-of-range positions raise `IndexError` even without a noise model. - `qubits` (`list[int]`): the list of qubit indices taking part in readout noise modeling. An `IndexError` is raised when an index is out of range. **Measurement and sampling** diff --git a/docs/api/en/python/5_visualization/3_state_plots.md b/docs/api/en/python/5_visualization/3_state_plots.md index 806f9aed..ccfcbf9c 100644 --- a/docs/api/en/python/5_visualization/3_state_plots.md +++ b/docs/api/en/python/5_visualization/3_state_plots.md @@ -120,6 +120,8 @@ Raises: Compute the expectation values of the quantum state in the Pauli basis and render them as a bar chart. +By default, labels follow `PauliString`: the highest-index qubit is on the left and qubit 0 is on the right. For example, `IX` measures X on qubit 0 of a two-qubit state. Enabling `reverse_bits` reverses the characters in each label without changing the corresponding expectation value. + Parameters: - `state` (`QuantumState`): the quantum state to render. diff --git a/docs/api/en/rust/3_qis/3_density_matrix_noise.md b/docs/api/en/rust/3_qis/3_density_matrix_noise.md index fbf788ba..a00c7174 100644 --- a/docs/api/en/rust/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/rust/3_qis/3_density_matrix_noise.md @@ -87,7 +87,7 @@ Methods: The simulator expresses "noise on the quantum state" and "noise on readout" separately: the former acts on the density matrix itself, while the latter only changes the reported probability distribution. - `fn probabilities(&self) -> Vec`: return the diagonal entries of the density matrix, that is the measurement probability of each computational basis state. Includes gate noise already applied and excludes readout noise. -- `fn probabilities_with_readout(&self, qubits: &[usize]) -> Result, QisError>`: based on `probabilities`, apply readout noise to the qubits listed in `qubits`. +- `fn probabilities_with_readout(&self, qubits: &[usize]) -> Result, QisError>`: `qubits` contains storage positions on which to apply readout errors. Circuit-bound states look up noise using the circuit ID at each position; unbound states use identity IDs. Returns all `2^num_qubits` probabilities in storage order; an empty selection applies no readout errors. Out-of-range positions return `QisError::IndexOutOfBounds` even without a noise model. - `fn probs(&self, measurement: &Measurement) -> Result, QisError>`: return the marginal probability distribution selected by `Measurement`, excluding readout noise. - `fn probs_with_readout(&self, measurement: &Measurement) -> Result, QisError>`: based on `probs`, apply readout noise to `measurement.qubits()`. diff --git a/docs/api/en/rust/5_visualization/3_state_plots.md b/docs/api/en/rust/5_visualization/3_state_plots.md index 6181fb22..b199bc6d 100644 --- a/docs/api/en/rust/5_visualization/3_state_plots.md +++ b/docs/api/en/rust/5_visualization/3_state_plots.md @@ -118,6 +118,8 @@ assert!(svg.contains("Re[rho]")); Draw Pauli basis expectation values as a bar chart. The labels are tensor products of `I`, `X`, `Y` and `Z`, `4^n` of them in total, where `n` is the number of qubits. Non-negative bars and negative bars use different colors. +By default, labels follow `PauliString`: the highest-index qubit is on the left and qubit 0 is on the right. For example, `IX` measures X on qubit 0 of a two-qubit state. Enabling `reverse_bits` reverses the characters in each label without changing the corresponding expectation value. + Parameters: - `state` (`&S`, `S: StateVisualizationSource + ?Sized`): the input state. diff --git a/docs/documentation/cn/1_cqlib/5_visualization/7_state_visualization.md b/docs/documentation/cn/1_cqlib/5_visualization/7_state_visualization.md index 539b52ea..7d8c5452 100644 --- a/docs/documentation/cn/1_cqlib/5_visualization/7_state_visualization.md +++ b/docs/documentation/cn/1_cqlib/5_visualization/7_state_visualization.md @@ -242,7 +242,9 @@ plot_state_paulivec( ![Asymmetric Pauli vector with reversed bits](assets/asymmetric_paulivec_reverse_bits.png) -在多比特状态图中，先确认 basis label 和 Pauli string 的显示顺序，再解释相干项、对角权重或 Pauli 期望值。 +默认 Pauli vector 标签遵循 `PauliString`：最高位在左，qubit 0 在最右侧。因此，`IX` 表示 qubit 0 上的 X，`XI` 表示 qubit 1 上的 X。上面的对角态在默认图中给出 `IZ = 0.6`、`ZI = -0.6`。设置 `reverse_bits=True` 时，每个标签中的字符顺序反转，对应的系数不变。 + +早期版本的多比特 Pauli 标签方向与此相反，现已修正。此前用 `reverse_bits=True` 对齐公共 `PauliString` 约定的代码，现在应使用默认设置。单比特图不受影响。 --- diff --git a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec.png b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec.png index 59102fa9..8675c36b 100644 Binary files a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec.png and b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec.png differ diff --git a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png index 7168b937..b2e9be6f 100644 Binary files a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png and b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png differ diff --git a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png index 7c4af68a..5878ed4a 100644 Binary files a/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png and b/docs/documentation/cn/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png differ diff --git a/docs/documentation/en/1_cqlib/5_visualization/7_state_visualization.md b/docs/documentation/en/1_cqlib/5_visualization/7_state_visualization.md index 6122d0a6..5307b0ba 100644 --- a/docs/documentation/en/1_cqlib/5_visualization/7_state_visualization.md +++ b/docs/documentation/en/1_cqlib/5_visualization/7_state_visualization.md @@ -242,7 +242,9 @@ Pauli vector after reversing the display order: ![Asymmetric Pauli vector with reversed bits](assets/asymmetric_paulivec_reverse_bits.png) -In multi-qubit state figures, confirm the display order of the basis labels and Pauli strings before interpreting coherence terms, diagonal weights or Pauli expectation values. +By default, Pauli-vector labels follow `PauliString`: the highest-index qubit is on the left and qubit 0 is on the right. Thus `IX` measures X on qubit 0 and `XI` measures X on qubit 1. For the diagonal state above, the default plot gives `IZ = 0.6` and `ZI = -0.6`. With `reverse_bits=True`, each label's characters are reversed while its coefficient stays the same. + +Earlier versions displayed multi-qubit Pauli labels in the opposite direction. This has been corrected; code that used `reverse_bits=True` to match the public `PauliString` convention should now use the default. Single-qubit plots are unaffected. --- diff --git a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec.png b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec.png index 59102fa9..8675c36b 100644 Binary files a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec.png and b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec.png differ diff --git a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png index 7168b937..b2e9be6f 100644 Binary files a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png and b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_report.png differ diff --git a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png index 7c4af68a..5878ed4a 100644 Binary files a/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png and b/docs/documentation/en/1_cqlib/5_visualization/assets/asymmetric_paulivec_reverse_bits.png differ diff --git a/docs/releasing.md b/docs/releasing.md index 2d2fa6a7..d828ba1a 100644 --- a/docs/releasing.md +++ b/docs/releasing.md @@ -39,9 +39,9 @@ Tags use the Python version to identify a release batch: - Stable releases: `vX.Y.Z`. - Prereleases: `vX.Y.Z-alpha.N`, `vX.Y.Z-beta.N`, or `vX.Y.Z-rc.N`, with N starting at 1. -- Current tag: `v2.0.0-beta.5`. -- Current Python source version: `2.0.0-beta5`, normalized to `2.0.0b5` for PEP 440 and wheel filenames. -- Current Rust/C version: `0.1.0-beta.5`; both inherit the workspace version. +- Planned release tag: `v2.0.0-beta.6` (create manually after validation). +- Current Python source version: `2.0.0-beta6`, normalized to `2.0.0b6` for PEP 440 and wheel filenames. +- Current Rust/C version: `0.1.0-beta.6`; both inherit the workspace version. Python and Rust/C versions evolve independently and do not need matching numbers. The Cargo.toml files at a given tag record the version mapping for that @@ -50,6 +50,25 @@ CI checks that the tag is equivalent to the Python version, core dependency versions are consistent across bindings, the minimum Python version is 3.11, and `abi3-py311` is enabled. Maintainers must update versions before tagging. +## Beta 6 Compatibility Notes + +Beta 6 fixes composite-gate parameter binding and qubit ordering, nested global +phase preservation, symbolic parameter substitution, density-matrix target +validation, readout-noise qubit lookup, NumPy complex Hamiltonian coefficients, +and Pauli-vector labels. Callers upgrading from beta 5 should account for: + +- `GPhase` matrices now have shape 1×1, consistent with a zero-qubit gate. +- Gate-local instruction and operation matrices put the first qubit argument + at the most-significant bit, including circuit-backed gates. Full circuit + matrices retain their least-significant-bit-first basis. +- Parameter substitution is simultaneous: replacement expressions are not + substituted again by other entries in the same binding map. +- Pauli-vector labels put qubit 0 at the right by default, matching `PauliString`. +- Readout-noise selections use storage positions; circuit-bound states resolve + the configured noise using the circuit qubit ID at each selected position. +- Density-matrix gate operations reject duplicate targets before changing the + state. + ## Platforms and Artifacts | Platform | Runner / Build Environment | diff --git a/tools/c-sdk/CMakeLists.txt b/tools/c-sdk/CMakeLists.txt index 4fc74dd4..3461149d 100644 --- a/tools/c-sdk/CMakeLists.txt +++ b/tools/c-sdk/CMakeLists.txt @@ -28,7 +28,7 @@ else() endif() foreach(kind static shared) - foreach(suite circuit qis_algebra) + foreach(suite circuit qis_algebra density_matrix_targets) set(test_target test_${suite}_${kind}) add_executable(${test_target} test_${suite}.c) set_target_properties(${test_target} PROPERTIES C_STANDARD 11 C_STANDARD_REQUIRED YES) diff --git a/tools/c-sdk/README.md b/tools/c-sdk/README.md index ae465515..64fbe60a 100644 --- a/tools/c-sdk/README.md +++ b/tools/c-sdk/README.md @@ -56,8 +56,10 @@ directory is needed by SDK consumers. Circuit construction, parameters, fermionic algebra, Jordan-Wigner mapping and checked sparse Pauli algebra are documented in the generated header and -illustrated in `tests/test_circuit.c` and `tests/test_qis_algebra.c`. Both suites -run against static and shared libraries. New checked algebra APIs return error +illustrated in `tests/test_circuit.c` and `tests/test_qis_algebra.c`. +`tests/test_density_matrix_targets.c` checks rejection of duplicate targets and +preservation of the density matrix. All three suites run against static and +shared libraries. New checked algebra APIs return error codes and deliver owned results through output pointers; failed calls leave outputs unchanged. Initialize object output slots to NULL and free results with their corresponding functions. diff --git a/tools/ci.py b/tools/ci.py index ad2257ce..fc5cafea 100644 --- a/tools/ci.py +++ b/tools/ci.py @@ -105,7 +105,11 @@ def build_and_test(platform=None): shutil.copy2(ROOT / "LICENSE.txt", sdk) shutil.copy2(ROOT / "tools/c-sdk/README.md", sdk / "README.md") (sdk / "tests").mkdir() - for test in ("test_circuit.c", "test_qis_algebra.c"): + for test in ( + "test_circuit.c", + "test_qis_algebra.c", + "test_density_matrix_targets.c", + ): shutil.copy2(ROOT / "crates/binding-c/tests" / test, sdk / "tests") shutil.copy2(ROOT / "tools/c-sdk/CMakeLists.txt", sdk / "tests") if sys.platform == "darwin":