diff --git a/qiskit_dynamics/perturbation/dyson_magnus.py b/qiskit_dynamics/perturbation/dyson_magnus.py index b6de8338a..7642bcd25 100644 --- a/qiskit_dynamics/perturbation/dyson_magnus.py +++ b/qiskit_dynamics/perturbation/dyson_magnus.py @@ -804,13 +804,13 @@ def _get_q_term_list(complete_index_multisets: List[Multiset]) -> List: def _get_dyson_lmult_rule( complete_index_multisets: List[Multiset], perturbation_labels: Optional[List[Multiset]] = None ) -> List: - """Given a complete list of index multisets, return - the lmult rule in the format required for ``CustomProduct``. - Note, the generator :math:`G(t)` is encoded as index ``-1``, as - it will be prepended to the list of A matrices. + """Given a complete list of index multisets, return the lmult rule in the format required for + ``CustomProduct``. Note, the generator :math:`G(t)` is encoded as index ``-1``, as it will be + prepended to the list of A matrices. Similarly the index of the solution to the LMDE with + generator :math:`G(t)` is encoded as ``-1``. - While not required within the logic of this function, the input - should be canonically ordered according to ``_get_all_submultisets``. + While not required within the logic of this function, the input should be canonically ordered + according to ``_get_all_submultisets``. Args: complete_index_multisets: List of complete multisets. @@ -834,8 +834,15 @@ def _get_dyson_lmult_rule( lmult_rule = [(np.array([1.0]), np.array([[-1, -1]]))] for term_idx, term in enumerate(complete_index_multisets): - if len(term) == 1: - lmult_rule.append((np.array([1.0, 1.0]), np.array([[-1, term_idx], [term_idx, -1]]))) + if len(term) == 1 and term in perturbation_labels: + # if term not in perturbation_labels for len(term) == 1 the corresponding Dyson integral + # will always be 0 + lmult_rule.append( + ( + np.array([1.0, 1.0]), + np.array([[-1, term_idx], [perturbation_labels.index(term), -1]]), + ) + ) else: # self multiplied by base generator lmult_indices = [[-1, term_idx]] diff --git a/qiskit_dynamics/perturbation/multiset_utils.py b/qiskit_dynamics/perturbation/multiset_utils.py index 79f83db18..04fa1b235 100644 --- a/qiskit_dynamics/perturbation/multiset_utils.py +++ b/qiskit_dynamics/perturbation/multiset_utils.py @@ -15,7 +15,7 @@ Utility functions for working with multisets. """ -from typing import List, Optional, Tuple +from typing import List, Optional, Tuple, Iterable import itertools from multiset import Multiset @@ -46,15 +46,51 @@ def _multiset_to_sorted_list(multiset: Multiset) -> List: return sorted_list -def _sorted_multisets(multisets: List[Multiset]) -> List[Multiset]: - """Sort in non-decreasing order according to: - - ms1 <= ms2 if len(ms1) < len(ms2), or if - len(ms1) == len(ms2) and if - str(_multiset_to_sorted_list(ms1)) <= str(_multiset_to_sorted_list(ms2)). +class _MultisetSortKey: + """Dummy class for usage as a key when sorting Multiset instances. This assumes the elements + of the multisets can themselves be sorted. """ - return sorted(multisets, key=lambda x: str(len(x)) + ", " + str(_multiset_to_sorted_list(x))) + __slots__ = ("multiset",) + + def __init__(self, multiset: Multiset): + self.multiset = multiset + + def __lt__(self, other: Multiset) -> bool: + """Implements an ordering on multisets. + + This orders first according to length (the number of elements in each multiset). If ``self`` + and ``other`` are the same length, ``self < other`` if, when written as fully expanded and + sorted lists, ``self < other`` in lexicographic ordering. E.g. it holds that ``Multiset({0: + 2, 1: 1}) < Multiset({0: 1, 1: 2})``, as the list versions are ``x = [0, 0, 1]``, and ``y = + [0, 1, 1]``. Here ``x[0] == y[0]``, but ``x[1] < y[1]``, and hence ``x < y`` in this + ordering. + """ + if len(self.multiset) < len(other.multiset): + return True + + if len(other.multiset) < len(self.multiset): + return False + + unique_entries = set(self.multiset.distinct_elements()) + unique_entries.update(other.multiset.distinct_elements()) + unique_entries = sorted(unique_entries) + + for element in unique_entries: + self_count = self.multiset[element] + other_count = other.multiset[element] + + if self_count != other_count: + return self_count > other_count + + return False + + +def _sorted_multisets(multisets: Iterable[Multiset]) -> List[Multiset]: + """Sort in non-decreasing order according to the ordering described in the dummy class + _MultisetSort. + """ + return sorted(multisets, key=_MultisetSortKey) def _clean_multisets(multisets: List[Multiset]) -> List[Multiset]: diff --git a/releasenotes/notes/multiset-order-bug-fix-1f1603ee1e230cba.yaml b/releasenotes/notes/multiset-order-bug-fix-1f1603ee1e230cba.yaml new file mode 100644 index 000000000..eb1110ce1 --- /dev/null +++ b/releasenotes/notes/multiset-order-bug-fix-1f1603ee1e230cba.yaml @@ -0,0 +1,6 @@ +--- +fixes: + - | + Fixes a bug in the perturbation module with internal sorting of ``Multiset`` instances, which + caused incorrect computation of perturbation theory terms when ``>10`` perturbations are + present. diff --git a/test/dynamics/perturbation/test_dyson_magnus.py b/test/dynamics/perturbation/test_dyson_magnus.py index dd0943b26..3d812cba4 100644 --- a/test/dynamics/perturbation/test_dyson_magnus.py +++ b/test/dynamics/perturbation/test_dyson_magnus.py @@ -555,6 +555,31 @@ def test__get_dyson_lmult_rule_power_series_case1(self): perturbation_labels=perturbation_labels, ) + def test__get_dyson_lmult_rule_power_series_case1_missing(self): + """Test _get_dyson_lmult_rule for higher order terms in the generator + decomposition case 1 where there is no corresponding perturbation_label term for + a desired expansion term. + """ + + expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]] + expansion_labels = [Multiset(label) for label in expansion_labels] + perturbation_labels = [[0], [0, 1], [1, 1]] + perturbation_labels = [Multiset(label) for label in perturbation_labels] + expected_lmult_rule = [ + (np.ones(1, dtype=float), np.array([[-1, -1]])), + (np.ones(2, dtype=float), np.array([[-1, 0], [0, -1]])), + (np.ones(1, dtype=float), np.array([[-1, 1]])), + (np.ones(3, dtype=float), np.array([[-1, 2], [0, 1], [1, -1]])), + (np.ones(2, dtype=float), np.array([[-1, 3], [2, -1]])), + (np.ones(4, dtype=float), np.array([[-1, 4], [0, 3], [1, 1], [2, 0]])), + ] + + self._test__get_dyson_lmult_rule( + expansion_labels, + expected_lmult_rule, + perturbation_labels=perturbation_labels, + ) + def test__get_dyson_lmult_rule_power_series_case2(self): """Test _get_dyson_lmult_rule for higher order terms in the generator decomposition case 2. @@ -579,6 +604,36 @@ def test__get_dyson_lmult_rule_power_series_case2(self): perturbation_labels=perturbation_labels, ) + def test__get_dyson_lmult_rule_power_series_case2_unordered(self): + """Test _get_dyson_lmult_rule for higher order terms in the generator + decomposition case 2, with the perturbation_labels being non-canonically ordered. + + Note that the conversion of case2 to case2_unordered requires relabelling 1 <-> 3 + in expected_lmult_rule, but also changing the ordering of the pairs in the matrix-product + description of expected_lmult_rule, as the rule is constructed by iterating through the + entries of perturbation_labels. The matrix-product description must be ordered by the first + index. + """ + + expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]] + expansion_labels = [Multiset(label) for label in expansion_labels] + perturbation_labels = [[0], [0, 1], [2], [1]] + perturbation_labels = [Multiset(label) for label in perturbation_labels] + expected_lmult_rule = [ + (np.ones(1, dtype=float), np.array([[-1, -1]])), + (np.ones(2, dtype=float), np.array([[-1, 0], [0, -1]])), + (np.ones(2, dtype=float), np.array([[-1, 1], [3, -1]])), + (np.ones(4, dtype=float), np.array([[-1, 2], [0, 1], [1, -1], [3, 0]])), + (np.ones(2, dtype=float), np.array([[-1, 3], [3, 1]])), + (np.ones(4, dtype=float), np.array([[-1, 4], [0, 3], [1, 1], [3, 2]])), + ] + + self._test__get_dyson_lmult_rule( + expansion_labels, + expected_lmult_rule, + perturbation_labels=perturbation_labels, + ) + def test__get_dyson_lmult_rule_case1(self): """Test __get_dyson_lmult_rule case 1.""" expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]] @@ -652,11 +707,10 @@ def test__get_dyson_lmult_rule_case2(self): self._test__get_dyson_lmult_rule(expansion_labels, expected_lmult_rule) def _test__get_dyson_lmult_rule( - self, complete_symmetric_dyson_term_list, expected, perturbation_labels=None + self, complete_index_multisets, expected, perturbation_labels=None ): - """Run a test case for _get_symmetric_dyson_mult_rules.""" - lmult_rule = _get_dyson_lmult_rule(complete_symmetric_dyson_term_list, perturbation_labels) - + """Run a test case for _get_dyson_lmult_rule.""" + lmult_rule = _get_dyson_lmult_rule(complete_index_multisets, perturbation_labels) self.assertMultRulesEqual(lmult_rule, expected) def assertMultRulesEqual(self, rule1, rule2): diff --git a/test/dynamics/perturbation/test_multiset_utils.py b/test/dynamics/perturbation/test_multiset_utils.py index 8be47b188..38526a6e6 100644 --- a/test/dynamics/perturbation/test_multiset_utils.py +++ b/test/dynamics/perturbation/test_multiset_utils.py @@ -59,6 +59,15 @@ def test_case1(self): expected = [Multiset([1]), Multiset([0, 2]), Multiset([0, 0, 1]), Multiset([0, 1, 1])] self.assertTrue(output == expected) + def test_case2(self): + """Test case which would have caught a bug introduced by a previous implementation of the + ordering. + """ + multisets = [Multiset([10]), Multiset([1])] + output = _sorted_multisets(multisets) + expected = [Multiset([1]), Multiset([10])] + self.assertTrue(output == expected) + class TestToSortedList(QiskitDynamicsTestCase): """Test conversion to sorted list.""" diff --git a/test/dynamics/perturbation/test_solve_lmde_perturbation.py b/test/dynamics/perturbation/test_solve_lmde_perturbation.py index 724061619..1e70a556d 100644 --- a/test/dynamics/perturbation/test_solve_lmde_perturbation.py +++ b/test/dynamics/perturbation/test_solve_lmde_perturbation.py @@ -501,6 +501,81 @@ def A1(t): self.assertAllClose(expected_D0001, results.perturbation_data.get_item([0, 0, 0, 1])[-1]) self.assertAllClose(expected_D0011, results.perturbation_data.get_item([0, 0, 1, 1])[-1]) + def test_dyson_analytic_case1_1d_relabeled(self): + """This is the same numerical test as test_dyson_analytic_case1_1d, however it relabels + the perturbation indices as 0 -> 1, and 1 -> 10. This is an integration test that would + have caught multiset sorting bug that broke solve_lmde_perturbation results computation. + """ + + def generator(t): + return Array([[1, 0], [0, 1]], dtype=complex).data + + def A0(t): + return Array([[0, t], [t**2, 0]], dtype=complex).data + + def A1(t): + return Array([[t, 0], [0, t**2]], dtype=complex).data + + T = np.pi * 1.2341 + + results = solve_lmde_perturbation( + perturbations=[A0, A1], + t_span=[0, T], + generator=generator, + y0=np.array([0.0, 1.0], dtype=complex), + expansion_method="dyson", + expansion_order=2, + expansion_labels=[[1, 1, 10], [1, 1, 1, 10], [1, 1, 10, 10]], + perturbation_labels=[[1], [10]], + dyson_in_frame=False, + integration_method=self.integration_method, + atol=1e-13, + rtol=1e-13, + ) + + T2 = T**2 + T3 = T * T2 + T4 = T * T3 + T5 = T * T4 + T6 = T * T5 + T7 = T * T6 + T8 = T * T7 + T9 = T * T8 + T10 = T * T9 + T11 = T * T10 + + U = expm(np.array(generator(0)) * T) + + expected_D0 = U @ np.array([[T2 / 2], [0]], dtype=complex) + expected_D1 = U @ np.array([[0], [T3 / 3]], dtype=complex) + expected_D00 = U @ np.array([[0], [T5 / 10]], dtype=complex) + expected_D01 = U @ ( + np.array([[T5 / 15], [0]], dtype=complex) + np.array([[T4 / 8], [0]], dtype=complex) + ) + expected_D11 = U @ np.array([[0], [T6 / 18]], dtype=complex) + expected_D001 = U @ ( + np.array([[0], [T8 / 120]], dtype=complex) + + np.array([[0], [T7 / 56]], dtype=complex) + + np.array([[0], [T8 / 80]], dtype=complex) + ) + expected_D0001 = U @ np.array([[(T10 / 480) + (T9 / 280)], [0]], dtype=complex) + expected_D0011 = U @ np.array( + [ + [0], + [(T11 / 396) + 23 * (T10 / 8400) + (T9 / 432)], + ], + dtype=complex, + ) + + self.assertAllClose(expected_D0, results.perturbation_data.get_item([1])[-1]) + self.assertAllClose(expected_D1, results.perturbation_data.get_item([10])[-1]) + self.assertAllClose(expected_D00, results.perturbation_data.get_item([1, 1])[-1]) + self.assertAllClose(expected_D01, results.perturbation_data.get_item([1, 10])[-1]) + self.assertAllClose(expected_D11, results.perturbation_data.get_item([10, 10])[-1]) + self.assertAllClose(expected_D001, results.perturbation_data.get_item([1, 1, 10])[-1]) + self.assertAllClose(expected_D0001, results.perturbation_data.get_item([1, 1, 1, 10])[-1]) + self.assertAllClose(expected_D0011, results.perturbation_data.get_item([1, 1, 10, 10])[-1]) + def test_dyson_analytic_case1(self): """Analytic test of computing dyson terms.