Pauli expectation values method (unitaryfund#1006)

* ExpectationPauliAll()

* PauliExpectation() wrapping

* Debug

* Self-review and bump version

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required (VERSION 3.9)
-project (Qrack VERSION 9.6.7 DESCRIPTION "High Performance Quantum Bit Simulation" LANGUAGES CXX)
+project (Qrack VERSION 9.7.0 DESCRIPTION "High Performance Quantum Bit Simulation" LANGUAGES CXX)
 
 # Installation commands
 include (GNUInstallDirs)

doxygen.config

@@ -38,7 +38,7 @@ PROJECT_NAME           = "Qrack"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 9.6
+PROJECT_NUMBER         = 9.7
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

include/pinvoke_api.hpp

@@ -72,6 +72,8 @@ MICROSOFT_QUANTUM_DECL double FactorizedExpectation(
     _In_ uintq sid, _In_ uintq n, _In_reads_(n) uintq* q, _In_ uintq m, uintq* c);
 MICROSOFT_QUANTUM_DECL double FactorizedExpectationRdm(
     _In_ uintq sid, _In_ uintq n, _In_reads_(n) uintq* q, _In_ uintq m, uintq* c, _In_ bool r);
+MICROSOFT_QUANTUM_DECL double PauliExpectation(
+    _In_ uintq sid, _In_ uintq n, _In_reads_(n) uintq* q, _In_reads_(n) uintq* b);
 #if FPPOW < 6
 MICROSOFT_QUANTUM_DECL double FactorizedExpectationFp(_In_ uintq sid, _In_ uintq n, _In_reads_(n) uintq* q, float* c);
 MICROSOFT_QUANTUM_DECL double FactorizedExpectationFpRdm(

include/qinterface.hpp

@@ -2470,6 +2470,16 @@ class QInterface : public ParallelFor {
         return ExpectationBitsFactorized(bits, perms, offset);
     }
 
+    /**
+     * Get permutation expectation value of bits
+     *
+     * The permutation expectation value of all included bits is returned, with bits valued from low to high as the
+     * order of the "bits" array parameter argument.
+     *
+     * \warning PSEUDO-QUANTUM
+     */
+    virtual real1_f ExpectationPauliAll(std::vector<bitLenInt> bits, std::vector<Pauli> paulis);
+
     /**
      * Get expectation value of bits, given an array of qubit weights
      *

include/wasm_api.hpp

@@ -190,6 +190,10 @@ real1_f FactorizedExpectation(quid sid, std::vector<QubitIntegerExpectation> q);
  * expectation value
  */
 real1_f FactorizedExpectationRdm(quid sid, std::vector<QubitIntegerExpectation> q, bool r);
+/**
+ * Pauli operator expectation value for the array of qubits and bases.
+ */
+real1_f PauliExpectation(quid sid, std::vector<bitLenInt> q, std::vector<Pauli> b);
 /**
  * Expectation value for bit-string integer from group of qubits with per-qubit real1 expectation value
  */

src/pinvoke_api.cpp

@@ -2592,6 +2592,26 @@ MICROSOFT_QUANTUM_DECL double FactorizedExpectationRdm(
     return _FactorizedExpectation(sid, n, q, m, c, NULL, r, true);
 }
 
+/**
+ * (External API) Get the Pauli operator expectation value for the array of qubits and bases.
+ */
+MICROSOFT_QUANTUM_DECL double PauliExpectation(
+    _In_ uintq sid, _In_ uintq n, _In_reads_(n) uintq* q, _In_reads_(n) uintq* b)
+{
+    SIMULATOR_LOCK_GUARD_DOUBLE(sid)
+
+    std::vector<bitLenInt> _q;
+    std::vector<Pauli> _b;
+    _q.reserve(n);
+    _b.reserve(n);
+    for (size_t i = 0U; i < n; ++i) {
+        _q.emplace_back((bitLenInt)q[i]);
+        _b.emplace_back((Pauli)b[i]);
+    }
+
+    return simulator->ExpectationPauliAll(_q, _b);
+}
+
 /**
  * (External API) Get the permutation expectation value, based upon the order of input qubits.
  */

src/qinterface/qinterface.cpp

@@ -487,6 +487,62 @@ real1_f QInterface::VarianceBitsAll(const std::vector<bitLenInt>& bits)
     return tot;
 }
 
+real1_f QInterface::ExpectationPauliAll(std::vector<bitLenInt> bits, std::vector<Pauli> paulis)
+{
+    for (size_t i = 0U; i < bits.size(); ++i) {
+        const size_t j = bits.size() - (i + 1U);
+        if (paulis[j] == PauliI) {
+            bits.erase(bits.begin() + j);
+            paulis.erase(paulis.begin() + j);
+        }
+    }
+
+    if (bits.empty()) {
+        return ONE_R1;
+    }
+
+    std::vector<real1> eigenVals;
+    eigenVals.reserve(bits.size() << 1U);
+    for (size_t i = 0U; i < bits.size(); ++i) {
+        eigenVals.push_back(ONE_R1);
+        eigenVals.push_back(-ONE_R1);
+
+        switch (paulis[i]) {
+        case PauliX:
+            H(bits[i]);
+            break;
+        case PauliY:
+            IS(bits[i]);
+            H(bits[i]);
+            break;
+        case PauliZ:
+        case PauliI:
+        default:
+            break;
+        }
+    }
+
+    const real1_f toRet = ExpectationFloatsFactorized(bits, eigenVals);
+
+    for (size_t i = 0U; i < bits.size(); ++i) {
+        switch (paulis[i]) {
+        case PauliX:
+            H(bits[i]);
+            break;
+        case PauliY:
+            H(bits[i]);
+            S(bits[i]);
+            break;
+        case PauliZ:
+        case PauliI:
+        default:
+            break;
+        }
+    }
+
+    return toRet;
+}
+
 real1_f QInterface::ExpectationBitsFactorized(
     const std::vector<bitLenInt>& bits, const std::vector<bitCapInt>& perms, bitCapInt offset)
 {
@@ -499,6 +555,10 @@ real1_f QInterface::ExpectationBitsFactorized(
         "QInterface::ExpectationBitsFactorized() parameter qubits vector values must be within allocated qubit "
         "bounds!");
 
+    if (bits.empty()) {
+        return ONE_R1;
+    }
+
     if (bits.size() == 1U) {
         const real1_f prob = Prob(bits[0]);
         return (
@@ -531,6 +591,10 @@ real1_f QInterface::ExpectationFloatsFactorized(const std::vector<bitLenInt>& bi
         "QInterface::ExpectationFloatsFactorized() parameter qubits vector values must be within allocated qubit "
         "bounds!");
 
+    if (bits.empty()) {
+        return ONE_R1;
+    }
+
     if (bits.size() == 1U) {
         const real1_f prob = Prob(bits[0]);
         return weights[0U] * (ONE_R1_F - prob) + weights[1U] * prob;

src/wasm_api.cpp

@@ -1863,6 +1863,22 @@ real1_f FactorizedExpectationRdm(quid sid, std::vector<QubitIntegerExpectation>
     return simulator->ExpectationBitsFactorizedRdm(r, _q, _c);
 }
 
+/**
+ * (External API) Get the Pauli operator expectation value for the array of qubits and bases.
+ */
+real1_f PauliExpectation(quid sid, std::vector<bitLenInt> q, std::vector<Pauli> b)
+{
+    SIMULATOR_LOCK_GUARD_REAL1_F(sid)
+
+    std::vector<bitLenInt> _q;
+    _q.reserve(q.size());
+    for (size_t i = 0U; i < q.size(); ++i) {
+        _q.emplace_back(shards[simulators[sid].get()][q[i]]);
+    }
+
+    return simulator->ExpectationPauliAll(_q, b);
+}
+
 /**
  * (External API) Get the permutation expectation value, based upon the order of input qubits.
  */