【代码贡献】优化 QSEncode Walsh 变换与 Top-K 筛选性能

pyqpanda-algorithm/pyqpanda_alg/QSEncode/QSEncode.py

@@ -14,10 +14,55 @@
 import numpy as np
 from pyqpanda3.core import CPUQVM, QCircuit, QProg, Encode, H
 from scipy.fft import fft
-from sympy import fwht
 
 from .. plugin import *
 
+
+def _fast_walsh_hadamard_transform(values):
+    """Return the unnormalized Walsh-Hadamard transform of a 1-D array.
+
+    The input is zero-padded to the next power-of-two length, matching the
+    behavior of :func:`sympy.fwht`.  The butterfly is evaluated with NumPy
+    ufuncs, so numeric values stay in a native numeric dtype instead of being
+    converted to symbolic Python objects.
+
+    The input is never modified.  Runtime is ``O(N log N)`` and auxiliary
+    memory is ``O(N)``.
+    """
+    array = np.asarray(values)
+    if array.ndim != 1:
+        raise ValueError('values must be a 1D array')
+    if not np.issubdtype(array.dtype, np.number):
+        raise TypeError('values must contain numeric data')
+
+    size = array.size
+    dtype = np.result_type(array.dtype, np.float64)
+    if size == 0:
+        return np.empty(0, dtype=dtype)
+
+    padded_size = 1 << (size - 1).bit_length()
+    output = np.zeros(padded_size, dtype=dtype)
+    output[:size] = array
+    if padded_size == 1:
+        return output
+
+    # One reusable half-size buffer avoids allocating a Python/SymPy object
+    # for every coefficient and avoids a fresh temporary at every stage.
+    scratch = np.empty(padded_size // 2, dtype=dtype)
+    step = 1
+    while step < padded_size:
+        blocks = output.reshape(-1, 2 * step)
+        left = blocks[:, :step]
+        right = blocks[:, step:]
+        saved_left = scratch.reshape(-1, step)
+        np.copyto(saved_left, left)
+        np.add(saved_left, right, out=left)
+        np.subtract(saved_left, right, out=right)
+        step *= 2
+
+    return output
+
+
 class QSpare_Code:
 
     """
@@ -146,15 +191,50 @@ def select_top_n_complex_numbers(self, arr, n):
         -------
         np.ndarray
             A new array with only top-n magnitudes retained, rest set to zero.
+
+        Notes
+        -----
+        Selection uses ``numpy.argpartition`` and therefore takes average
+        ``O(N)`` time.  If equal magnitudes cross the selection boundary, the
+        larger original indices are retained to make the result deterministic.
         """
 
         if type(n) != int or n <= 0:
             raise ValueError('n must > 0 and with class int')
+        arr = np.asarray(arr)
+        if arr.ndim != 1:
+            raise ValueError('arr must be a 1D array')
+        if arr.size == 0:
+            return arr.copy()
+
         magnitudes = np.abs(arr)
-        top_n_indices = np.argsort(magnitudes)[-n:]
+        if not np.all(np.isfinite(magnitudes)):
+            raise ValueError('arr must contain only finite values')
+        if n >= arr.size:
+            return arr.copy()
+
+        # argpartition finds the top-n boundary in average O(N), whereas a
+        # complete argsort performs O(N log N) work.  Values tied at the
+        # boundary are resolved by index so results are deterministic.
+        partition_at = arr.size - n
+        partitioned = np.argpartition(magnitudes, partition_at)
+        threshold = magnitudes[partitioned[partition_at]]
+        top_n_indices = partitioned[partition_at:]
+
+        # The partition result is already complete for the overwhelmingly
+        # common unique-boundary case.  Only scan and rebuild the selection
+        # when equal magnitudes straddle the top-n boundary.
+        selected_ties = np.count_nonzero(magnitudes[top_n_indices] == threshold)
+        all_ties = np.count_nonzero(magnitudes == threshold)
+        if selected_ties != all_ties:
+            top_n_indices = np.flatnonzero(magnitudes > threshold)
+            slots_left = n - top_n_indices.size
+            tied_indices = np.flatnonzero(magnitudes == threshold)
+            top_n_indices = np.concatenate((top_n_indices, tied_indices[-slots_left:]))
+
         result = np.zeros_like(arr, dtype=arr.dtype)
         result[top_n_indices] = arr[top_n_indices]
-        return np.array(result)
+        return result
 
     def Transform(self, amp):
         """
@@ -169,9 +249,15 @@ def Transform(self, amp):
         -------
         np.ndarray
             Transformed amplitude vector in the selected basis.
+
+        Notes
+        -----
+        Walsh mode uses a native NumPy fast Walsh-Hadamard butterfly with
+        ``O(N log N)`` time and ``O(N)`` memory.
         """
         if self.mode == 'walsh':
-            transform = np.array(fwht(amp)) / np.sqrt(2 ** self.qubits_num)
+            transform = _fast_walsh_hadamard_transform(amp)
+            transform /= np.sqrt(2 ** self.qubits_num)
         elif self.mode == 'fourier':
             transform = fft(amp)
         else:

test/QAlgBase/Test_QSEncode_numeric_pipeline.py

@@ -0,0 +1,165 @@
+import numpy as np
+import pytest
+
+from pyqpanda_alg.QSEncode.QSEncode import (
+    QSpare_Code,
+    _fast_walsh_hadamard_transform,
+)
+
+
+def _hadamard_matrix(size):
+    """Build a small independent Sylvester Hadamard matrix."""
+    matrix = np.array([[1.0]])
+    while matrix.shape[0] < size:
+        matrix = np.block([[matrix, matrix], [matrix, -matrix]])
+    return matrix
+
+
+def _dense_fwht(values):
+    """Reference transform using an explicit matrix, not a butterfly."""
+    values = np.asarray(values)
+    if values.size == 0:
+        return np.empty(0, dtype=np.result_type(values.dtype, np.float64))
+    padded_size = 1 << (values.size - 1).bit_length()
+    padded = np.zeros(padded_size, dtype=np.result_type(values.dtype, np.float64))
+    padded[:values.size] = values
+    return _hadamard_matrix(padded_size) @ padded
+
+
+def _dense_sparse_walsh_probabilities(probabilities, cut):
+    """Independent end-to-end oracle for the Walsh sparse encoder."""
+    probabilities = np.asarray(probabilities, dtype=float)
+    amplitudes = np.sqrt(probabilities / probabilities.sum())
+    size = amplitudes.size
+    hadamard = _hadamard_matrix(size) / np.sqrt(size)
+    coefficients = hadamard @ amplitudes
+
+    order = np.argsort(np.abs(coefficients), kind='stable')
+    sparse = np.zeros_like(coefficients)
+    sparse[order[-cut:]] = coefficients[order[-cut:]]
+    sparse /= np.linalg.norm(sparse)
+    return np.abs(hadamard @ sparse) ** 2
+
+
+@pytest.mark.parametrize(
+    'values',
+    [
+        np.array([], dtype=float),
+        np.array([3.25]),
+        np.array([1.0, -2.0, 4.0]),
+        np.array([1.0, 2.0, 3.0, 4.0]),
+        np.array([1.0 + 2.0j, -3.0j, 0.5 - 0.25j]),
+    ],
+)
+def test_fast_walsh_transform_matches_dense_oracle(values):
+    actual = _fast_walsh_hadamard_transform(values)
+    expected = _dense_fwht(values)
+
+    assert actual.dtype != object
+    np.testing.assert_allclose(actual, expected, rtol=0.0, atol=1e-12)
+
+
+def test_fast_walsh_transform_matches_dense_oracle_for_lengths_1_to_64():
+    rng = np.random.default_rng(20260901)
+
+    for size in range(1, 65):
+        values = rng.normal(size=size)
+        if size % 2 == 0:
+            values = values + 1j * rng.normal(size=size)
+        actual = _fast_walsh_hadamard_transform(values)
+        expected = _dense_fwht(values)
+        np.testing.assert_allclose(actual, expected, rtol=1e-13, atol=1e-12)
+
+
+@pytest.mark.parametrize('qubits', range(13))
+def test_normalized_walsh_transform_preserves_norm_and_is_self_inverse(qubits):
+    size = 1 << qubits
+    values = np.random.default_rng(qubits).normal(size=size)
+    transformed = _fast_walsh_hadamard_transform(values) / np.sqrt(size)
+    restored = _fast_walsh_hadamard_transform(transformed) / np.sqrt(size)
+
+    np.testing.assert_allclose(np.linalg.norm(transformed), np.linalg.norm(values), rtol=1e-13)
+    np.testing.assert_allclose(restored, values, rtol=1e-13, atol=1e-13)
+
+
+def test_fast_walsh_transform_does_not_modify_input_and_promotes_numeric_dtype():
+    values = np.array([1, 2, 3], dtype=np.int16)
+    original = values.copy()
+
+    result = _fast_walsh_hadamard_transform(values)
+
+    np.testing.assert_array_equal(values, original)
+    assert result.dtype == np.float64
+
+
+def test_fast_walsh_transform_rejects_non_numeric_or_non_vector_input():
+    with pytest.raises(ValueError, match='1D'):
+        _fast_walsh_hadamard_transform(np.ones((2, 2)))
+    with pytest.raises(TypeError, match='numeric'):
+        _fast_walsh_hadamard_transform(['a', 'b'])
+
+
+@pytest.mark.parametrize(
+    ('size', 'cut'),
+    [(8, 1), (8, 7), (257, 17), (4096, 63)],
+)
+def test_linear_top_k_matches_full_stable_sort_for_unique_magnitudes(size, cut):
+    rng = np.random.default_rng(size + cut)
+    phases = np.exp(2j * np.pi * rng.random(size))
+    values = np.arange(1, size + 1, dtype=float) * phases
+    encoder = QSpare_Code([0.5, 0.5])
+
+    actual = encoder.select_top_n_complex_numbers(values, cut)
+    expected = np.zeros_like(values)
+    expected_indices = np.argsort(np.abs(values), kind='stable')[-cut:]
+    expected[expected_indices] = values[expected_indices]
+
+    np.testing.assert_array_equal(actual, expected)
+
+
+def test_linear_top_k_has_deterministic_boundary_ties():
+    values = np.array([-3.0, 3.0, 2.0, -2.0, 1.0])
+    encoder = QSpare_Code([0.5, 0.5])
+
+    actual = encoder.select_top_n_complex_numbers(values, 3)
+
+    np.testing.assert_array_equal(actual, [-3.0, 3.0, 0.0, -2.0, 0.0])
+
+
+def test_linear_top_k_keeps_input_unchanged_and_fast_paths_full_selection():
+    values = np.array([1.0 + 2.0j, -4.0j, 3.0])
+    original = values.copy()
+    encoder = QSpare_Code([0.5, 0.5])
+
+    result = encoder.select_top_n_complex_numbers(values, values.size)
+
+    np.testing.assert_array_equal(values, original)
+    np.testing.assert_array_equal(result, original)
+    assert result is not values
+
+
+def test_linear_top_k_rejects_invalid_shape_and_non_finite_values():
+    encoder = QSpare_Code([0.5, 0.5])
+    with pytest.raises(ValueError, match='1D'):
+        encoder.select_top_n_complex_numbers(np.ones((2, 2)), 1)
+    with pytest.raises(ValueError, match='finite'):
+        encoder.select_top_n_complex_numbers(np.array([1.0, np.nan]), 1)
+    with pytest.raises(ValueError, match='finite'):
+        encoder.select_top_n_complex_numbers(np.array([1.0, np.inf]), 2)
+
+
+@pytest.mark.parametrize(
+    ('probabilities', 'cut'),
+    [
+        ([0.50, 0.30, 0.15, 0.05], 2),
+        ([0.31, 0.19, 0.16, 0.12, 0.09, 0.06, 0.04, 0.03], 3),
+    ],
+)
+def test_walsh_quantum_result_matches_independent_dense_pipeline(probabilities, cut):
+    encoder = QSpare_Code(list(map(float, probabilities)), mode='walsh', cut_length=cut)
+
+    actual = encoder.Quantum_Res()
+    expected = _dense_sparse_walsh_probabilities(probabilities, cut)
+
+    np.testing.assert_allclose(actual, expected, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.sum(actual), 1.0, rtol=0.0, atol=1e-12)