added ordering field to PauliStrSum

quest/include/paulis.h

@@ -66,11 +66,16 @@ typedef struct {
 
     qindex numTerms;
 
-    // arbitrarily-sized collection of Pauli strings and their
-    // coefficients are stored in heap memory.
+    // numTerms-sized collection of Pauli strings and their
+    // corresponding coefficients, stored in heap memory.
     PauliStr* strings;
     qcomp* coeffs;
 
+    // numTerms-sized list of term indices, consulted by some
+    // routines (such as Trotter functions) to effectively
+    // reorder the terms (strings and coeffs)
+    qindex* ordering;
+
     // whether the sum constitutes a Hermitian operator (0, 1, or -1 to indicate unknown),
     // which is lazily evaluated when a function validates Hermiticity them. The flag is 
     // stored in heap so even copies of structs are mutable, but the pointer is immutable;
@@ -101,7 +106,7 @@ typedef struct {
  * @brief Functions for printing Pauli data structures.
  *
  * @defgroup paulis_setters Setters
- * @brief Functions for overwriting the elements of Pauli data structures.
+ * @brief Functions for modifying existing Pauli data structures.
  */
 
 
@@ -435,27 +440,26 @@ extern "C" {
 
     /** @ingroup paulis_setters
      *
-     * Reorders the terms within a @p sum of weighted Pauli strings to sort Pauli
-     * strings into lexicographic (dictionary) ordering.
+     * Reorders the terms within a @p sum of weighted Pauli strings so that Pauli strings 
+     * are ordered lexicographically. 
+     * 
+     * This affects @p sum.ordering, and will change the behaviour of functions like
+     * applyTrotterizedPauliStrSumGadget() when randomised ordering is disabled.
      *
      * @formulae
-     * Let @f$ H = @f$ @p sum, which can be represented as
+     * 
+     * Let @f$ H = @f$ @p sum, satisfying
      * @f[
           H = \sum\limits_j c_j \, \hat{\sigma}_j
      * @f]
      * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$.
      *
-     * This function constructs and applies the permutation @f$ \pi @f$ to @f$ H @f$
-     * @f[
-          H = \sum\limits_j c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
-     * @f]
-     * such that
+     * This function modifies @p sum.ordering to the permutation @f$ \pi @f$ such that
      * @f[
      *   \hat{\sigma}_{\pi(i)} <_{lex} \hat{\sigma}_{\pi(j)} \  \forall \  \pi(i) < \pi(j).
      * @f]
      *
-     *
-     * @param[in,out] sum  a weighted sum of Pauli strings to reorder.
+     * @param[in,out] sum  a weighted sum of Pauli strings to reorder (via @p sum.ordering)
      *
      * @throws @validationerror
      * - if @p sum is not initialised.
@@ -469,21 +473,21 @@ extern "C" {
 
     /** @ingroup paulis_setters
      *
-     * Reorders the terms within a @p sum of weighted Pauli strings to sort Pauli
-     * strings into decreasing magnitude weights.
+     * Reorders the terms within a @p sum of weighted Pauli strings such that
+     * coefficients are ordered with decreasing magnitude.
+     * 
+     * This affects only @p sum.ordering, and will change the behaviour of functions like
+     * applyTrotterizedPauliStrSumGadget() when randomised ordering is disabled.
      *
      * @formulae
-     * Let @f$ H = @f$ @p sum, represented as the weighted sum
+     *
+     * Let @f$ H = @f$ @p sum, satisfying
      * @f[
           H = \sum\limits_j c_j \, \hat{\sigma}_j
      * @f]
      * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$.
      *
-     * This function constructs and applies the permutation @f$ \pi @f$ to @f$ H @f$
-     * @f[
-          H = \sum\limits_j c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
-     * @f]
-     * such that
+     * This function modifies @p sum.ordering to the permutation @f$ \pi @f$ such that
      * @f[
      *    |c_{\pi(i)}| > |c_{\pi(j)}| \, \forall \,  \pi(i) < \pi(j).
      * @f]

quest/include/trotterisation.h

@@ -109,14 +109,16 @@ extern "C" {
  * > These formulations are taken from 'Finding Exponential Product Formulas
  * > of Higher Orders', Naomichi Hatano and Masuo Suzuki (2005) (<a href="https://arxiv.org/abs/math-ph/0506007">arXiv</a>).
  * 
- * When @p permutePaulis=true the terms of @p sum are effected in a random order at each repetition. That is, each repetition of the Trotter-Suzuki decomposition is evaluated with the sum
+ * When @p permutePaulis=true, the terms of @p sum are effected in a random order at each repetition. 
+ * That is, each repetition of the Trotter-Suzuki decomposition is evaluated with the sum
  * @f[
       \hat{H} = \sum\limits_j^T c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
  * @f]
- * where @f$ \pi @f$  is a randomly selected permutation.
+ * where @f$ \pi @f$  is a randomly selected permutation. When @p permutePaulis=false, the fixed
+ * ordering @f$ \pi = @f$ @p sum.ordering is used.
  *
  * @important
- *   Using @p permutePaulis=true will cause @p sum to be mutated by the Trotterisation.
+ *   Using @p permutePaulis=true will mutate the ordering of @p sum (specifically @p sum.ordering).
  *
  * @equivalences
  * 
@@ -153,7 +155,8 @@ extern "C" {
  *   when all PauliStr in @p sum = @f$ \hat{H} @f$ commute, or @p reps @f$ \rightarrow \infty @f$.
  * 
  * @param[in,out] qureg                the state to modify.
- * @param[in]     sum                  a weighted sum of Pauli strings to approximately exponentiate.
+ * @param[in,out] sum                  a weighted sum of Pauli strings to approximately exponentiate,
+ *                                     with ordering mutated when @p permutePaulis=true.
  * @param[in]     angle                the prefactor of @p sum times @f$ i @f$ in the exponent.
  * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
  * @param[in]     reps                 the number of Trotter repetitions.
@@ -248,7 +251,8 @@ void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* con
  *   when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. 
  * 
  * @param[in,out] qureg                the state to modify.
- * @param[in]     sum                  a weighted sum of Pauli strings to approximately exponentiate.
+ * @param[in,out] sum                  a weighted sum of Pauli strings to approximately exponentiate,
+ *                                     with ordering mutated when @p permutePaulis=true.
  * @param[in]     angle                an effective prefactor of @p sum in the exponent.
  * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
  * @param[in]     reps                 the number of Trotter repetitions.
@@ -391,7 +395,8 @@ extern "C" {
  *  - applyTrotterizedNonUnitaryPauliStrSumGadget()
  * 
  * @param[in,out] qureg                the state to modify.
- * @param[in]     hamil                the Hamiltonian as a a weighted sum of Pauli strings.
+ * @param[in,out] hamil                the Hamiltonian as a a weighted sum of Pauli strings,
+ *                                     with ordering mutated when @p permutePaulis=true.
  * @param[in]     time                 the duration over which to simulate evolution.
  * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
  * @param[in]     reps                 the number of Trotter repetitions.
@@ -522,7 +527,8 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal
  *  - applyTrotterizedNonUnitaryPauliStrSumGadget()
  * 
  * @param[in,out] qureg                the state to modify.
- * @param[in]     hamil                the Hamiltonian as a a weighted sum of Pauli strings.
+ * @param[in,out] hamil                the Hamiltonian as a a weighted sum of Pauli strings,
+ *                                     with ordering mutated when @p permutePaulis=true.
  * @param[in]     tau                  the duration over which to simulate imaginary-time evolution.
  * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
  * @param[in]     reps                 the number of Trotter repetitions.
@@ -548,6 +554,11 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea
  * evolution approximated by symmetrized Trotterisation of the specified @p order and number of cycles
  * @p reps.
  * 
+ * Note the ordering of all passed PauliStrSum (through functions like sortPauliStrSumMagnitude()) will
+ * affect that of the internally created super-propagator and ergo the Trotter accuracy. This is overridden
+ * by passing @p permutePaulis=true, whereby the super-propagator order is randomised every Trotter repetition.
+ * This never mutates the ordering of all passed PauliStrSum.
+ * 
  * @formulae 
  * 
  * Let @f$ \rho = @f$ @p qureg, @f$ \hat{H} = @f$ @p hamil, @f$ t = @f$ @p time, and denote the @f$ i @f$-th

quest/src/api/paulis.cpp

@@ -20,6 +20,8 @@
 
 #include <vector>
 #include <string>
+#include <numeric>
+#include <algorithm>
 
 using std::string;
 using std::vector;
@@ -34,8 +36,9 @@ using std::vector;
 bool didAnyAllocsFailOnAnyNode(PauliStrSum sum) {
 
     bool anyFail = (
-        ! mem_isAllocated(sum.strings) || 
-        ! mem_isAllocated(sum.coeffs)  || 
+        ! mem_isAllocated(sum.strings)  || 
+        ! mem_isAllocated(sum.coeffs)   || 
+        ! mem_isAllocated(sum.ordering) ||
         ! mem_isAllocated(sum.isApproxHermitian) );
 
     if (comm_isInit())
@@ -50,6 +53,7 @@ void freePauliStrSum(PauliStrSum sum) {
     // these do not need to be allocated (freeing nullptr is legal)
     cpu_deallocPauliStrings(sum.strings);
     cpu_deallocArray(sum.coeffs);
+    cpu_deallocIndices(sum.ordering);
     util_deallocEpsilonSensitiveHeapFlag(sum.isApproxHermitian);
 }
 
@@ -173,6 +177,7 @@ extern "C" PauliStrSum createPauliStrSum(PauliStr* strings, qcomp* coeffs, qinde
     out.numTerms = numTerms;
     out.strings = cpu_allocPauliStrings(numTerms);                // nullptr if failed
     out.coeffs  = cpu_allocArray(numTerms);                       // nullptr if failed
+    out.ordering = cpu_allocIndices(numTerms);                    // nullptr if failed
     out.isApproxHermitian = util_allocEpsilonSensitiveHeapFlag(); // nullptr if failed
 
     // if either alloc failed, clear both before validation to avoid leak
@@ -182,6 +187,7 @@ extern "C" PauliStrSum createPauliStrSum(PauliStr* strings, qcomp* coeffs, qinde
     // otherwise copy given data into new heap structure, and set initial flags
     cpu_copyPauliStrSum(out, strings, coeffs);
     util_setFlagToUnknown(out.isApproxHermitian);
+    std::iota(out.ordering, out.ordering + out.numTerms, 0);
 
     return out;
 }
@@ -308,7 +314,8 @@ extern "C" void sortPauliStrSumLexicographic(PauliStrSum sum) {
         return std::tie(strI.highPaulis, strI.lowPaulis) < std::tie(strJ.highPaulis, strJ.lowPaulis);
     };
 
-    paulis_sortTermsViaComparator(sum, lexSort);
+    std::iota(sum.ordering, sum.ordering + sum.numTerms, 0);
+    std::stable_sort(sum.ordering, sum.ordering + sum.numTerms, lexSort);
 }
 
 
@@ -319,5 +326,6 @@ extern "C" void sortPauliStrSumMagnitude(PauliStrSum sum) {
         return std::norm(sum.coeffs[i]) > std::norm(sum.coeffs[j]);
     };
 
-    paulis_sortTermsViaComparator(sum, magSort);
+    std::iota(sum.ordering, sum.ordering + sum.numTerms, 0);
+    std::stable_sort(sum.ordering, sum.ordering + sum.numTerms, magSort);
 }

quest/src/api/trotterisation.cpp

@@ -18,6 +18,7 @@
 #include "quest/src/core/randomiser.hpp"
 
 #include <vector>
+#include <numeric>
 
 using std::vector;
 
@@ -33,9 +34,12 @@ void internal_applyFirstOrderTrotterRepetition(
 ) {
     // apply each sum term as a gadget, in forward or reverse order
     for (qindex i=0; i<sum.numTerms; i++) {
+
+        // identify the PauliStr which will become the gadget's generator
         int j = reverse? sum.numTerms - i - 1 : i;
-        qcomp coeff = sum.coeffs[j];
-        PauliStr str = sum.strings[j];
+        int k = sum.ordering[j];
+        qcomp coeff = sum.coeffs[k];
+        PauliStr str = sum.strings[k];
 
         // effect |psi> -> exp(i angle * coeff * term)|psi>
         qcomp arg = angle * coeff;
@@ -100,16 +104,18 @@ void internal_applyAllTrotterRepetitions(
 
     // perform carefully-ordered sequence of gadgets
     for (int r=0; r<reps; r++){
+
+        // optionally shuffle term ordering, else use fixed user-set order
         if (permutePaulis)
             rand_permutePauliStrSum(sum);
+
         internal_applyHigherOrderTrotterRepetition(
             qureg, ketCtrlsVec, braCtrlsVec, statesVec, sum, arg, order, onlyLeftApply);
     }
 
     /// @todo
-    /// the accuracy of Trotterisation is greatly improved by randomisation
-    /// or (even sub-optimal) grouping into commuting terms. Should we 
-    /// implement these above or into another function?
+    /// the accuracy of Trotterisation is greatly improved by 
+    /// (even sub-optimal) grouping into commuting terms. 
 }
 
 qindex internal_getNumTotalSuperPropagatorTerms(PauliStrSum hamil, PauliStrSum* jumps, int numJumps) {
@@ -299,17 +305,25 @@ void applyTrotterizedNoisyTimeEvolution(
     // validate memory allocations for all super-propagator terms
     vector<PauliStr> superStrings;
     vector<qcomp> superCoeffs;
-    auto callbackString = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(PauliStr), __func__); };
-    auto callbackCoeff  = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(qcomp),    __func__); };
-    util_tryAllocVector(superStrings, numSuperTerms, callbackString);
-    util_tryAllocVector(superCoeffs,  numSuperTerms, callbackCoeff);
-
+    vector<qindex> superOrdering;
+    auto callbackString   = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(PauliStr), __func__); };
+    auto callbackCoeff    = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(qcomp),    __func__); };
+    auto callbackOrdering = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(qindex),   __func__); };
+    util_tryAllocVector(superStrings,  numSuperTerms, callbackString);
+    util_tryAllocVector(superCoeffs,   numSuperTerms, callbackCoeff);
+    util_tryAllocVector(superOrdering, numSuperTerms, callbackOrdering);
+
+    // construct super-propagator term by term, in an order affected by
+    // both hamil.ordering and jumps[].ordering (inside paulis_setters())
     qindex superTermInd = 0;
 
     // collect -i[H,rho] terms
     for (qindex n=0; n<hamil.numTerms; n++) {
-        PauliStr oldStr = hamil.strings[n];
-        qcomp oldCoeff = hamil.coeffs[n];
+
+        // consult hamil ordering during construction
+        qindex m = hamil.ordering[n];
+        PauliStr oldStr = hamil.strings[m];
+        qcomp oldCoeff = hamil.coeffs[m];
 
         // term of -i Id (x) H
         superStrings[superTermInd] = oldStr;
@@ -324,8 +338,8 @@ void applyTrotterizedNoisyTimeEvolution(
 
     // below we bind superStrings/Coeffs to a spoofed PauliStrSum to pass to paulis functions
     PauliStrSum temp;
-    int flagForDebugSafety = -1;
-    temp.isApproxHermitian = &flagForDebugSafety;
+    temp.ordering          = nullptr; // not consulted
+    temp.isApproxHermitian = nullptr; // not consulted
 
     // collect jump terms
     for (int n=0; n<numJumps; n++) {
@@ -358,11 +372,15 @@ void applyTrotterizedNoisyTimeEvolution(
     if (superTermInd != numSuperTerms)
         error_unexpectedNumLindbladSuperpropTerms();
 
+    // initialise superOrdering = {0, 1, 2, ...}
+    std::iota(superOrdering.begin(), superOrdering.end(), 0);
+
     // pass superpropagator terms as temporary PauliStrSum
     PauliStrSum superSum; 
     superSum.numTerms = numSuperTerms;
     superSum.strings = superStrings.data();
     superSum.coeffs = superCoeffs.data();
+    superSum.ordering = superOrdering.data();
     superSum.isApproxHermitian = nullptr; // will not be queried
 
     // effect exp(t S) = exp(x i S) | x=-i*time, left-multiplying only

quest/src/core/memory.cpp

@@ -437,6 +437,7 @@ bool mem_isOuterAllocated(qcomp*** ptr) { return isNonNull(ptr); }
 bool mem_isAllocated(int* heapflag)   { return isNonNull(heapflag); }
 bool mem_isAllocated(qcomp* array )   { return isNonNull(array); }
 bool mem_isAllocated(PauliStr* array) { return isNonNull(array); }
+bool mem_isAllocated(qindex* array)   { return isNonNull(array); }
 
 bool mem_isAllocated(qcomp** matrix, qindex numRows) {
 

quest/src/core/memory.hpp

@@ -112,6 +112,7 @@ bool mem_canPauliStrSumFitInMemory(qindex numTerms, qindex numBytesPerNode);
 
 bool mem_isAllocated(int* heapflag);
 bool mem_isAllocated(PauliStr* array);
+bool mem_isAllocated(qindex* array);
 bool mem_isAllocated(qcomp* array);
 bool mem_isAllocated(qcomp** matrix, qindex numRows);
 bool mem_isAllocated(qcomp*** matrixList, qindex numRows, int numMatrices);