Cache block fix

src/misc/warnings.hpp

@@ -30,9 +30,10 @@
   DISABLE_WARNING(-Wall)                                                       \
   DISABLE_WARNING(-Wextra)                                                     \
   DISABLE_WARNING(-Wshadow)                                                    \
-  DISABLE_WARNING(-Wfloat-equal)                                             \
+  DISABLE_WARNING(-Wfloat-equal)                                               \
   DISABLE_WARNING(-Wundef)                                                     \
   DISABLE_WARNING(-Wpedantic)                                                  \
+  DISABLE_WARNING(-Wredundant-decls)                                           \
   DISABLE_WARNING(-Wsign-compare)
 
 #define DISABLE_WARNING_POP DO_PRAGMA(GCC diagnostic pop)
@@ -43,4 +44,4 @@
 
 #endif
 
-#endif /* Guards */
+#endif /* Guards */

src/misc/wrap_thrust.hpp

@@ -0,0 +1,49 @@
+/**
+ * This code is part of Qiskit.
+ *
+ * (C) Copyright IBM 2020.
+ *
+ * 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.
+ */
+
+#ifndef _aer_misc_wrap_thrust_hpp_
+#define _aer_misc_wrap_thrust_hpp_
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC system_header
+#endif
+
+#include "misc/warnings.hpp"
+DISABLE_WARNING_PUSH
+#include <thrust/for_each.h>
+#include <thrust/complex.h>
+#include <thrust/inner_product.h>
+#include <thrust/transform.h>
+#include <thrust/transform_scan.h>
+#include <thrust/binary_search.h>
+#include <thrust/execution_policy.h>
+#include <thrust/functional.h>
+#include <thrust/tuple.h>
+#include <thrust/iterator/constant_iterator.h>
+#include <thrust/detail/vector_base.h>
+
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/fill.h>
+
+#ifdef AER_THRUST_CUDA
+#include <thrust/device_vector.h>
+#endif
+#include <thrust/host_vector.h>
+
+#include <thrust/system/omp/execution_policy.h>
+DISABLE_WARNING_POP
+
+#endif // inclusion guard

src/simulators/density_matrix/densitymatrix.hpp

@@ -125,6 +125,13 @@ class DensityMatrix : public UnitaryMatrix<data_t> {
   // outcome in [0, 2^num_qubits - 1]
   virtual double probability(const uint_t outcome) const override;
 
+  //-----------------------------------------------------------------------
+  // Expectation Values
+  //-----------------------------------------------------------------------
+
+  // Return Pauli expectation value
+  double expval_pauli(const reg_t &qubits, const std::string &pauli,const complex_t initial_phase=1.0) const;
+
 protected:
 
   // Construct a vectorized superoperator from a vectorized matrix
@@ -344,6 +351,55 @@ void DensityMatrix<data_t>::apply_toffoli(const uint_t qctrl0,
   BaseVector::apply_permutation_matrix(qubits, pairs);
 }
 
+template <typename data_t>
+double DensityMatrix<data_t>::expval_pauli(const reg_t &qubits,
+                                           const std::string &pauli,const complex_t initial_phase) const {
+
+  uint_t x_mask, z_mask, num_y, x_max;
+  std::tie(x_mask, z_mask, num_y, x_max) = QV::pauli_masks_and_phase(qubits, pauli);
+
+  // Special case for only I Paulis
+  if (x_mask + z_mask == 0) {
+    return std::real(BaseMatrix::trace());
+  }
+
+  // Size of density matrix 
+  const size_t nrows = BaseMatrix::rows_;
+  const size_t diag_stride = 1 + nrows;
+
+  // specialize x_max == 0
+  if (!x_mask) {
+    auto lambda = [&](const int_t i, double &val_re, double &val_im)->void {
+      (void)val_im; // unused
+      auto val = std::real(BaseVector::data_[i * diag_stride]);
+      if (z_mask && (AER::Utils::popcount(i & z_mask) & 1)) {
+        val = -val;
+      }
+      val_re += val;
+    };
+    return std::real(BaseVector::apply_reduction_lambda(std::move(lambda), size_t(0), nrows));
+  }
+
+  auto phase = std::complex<data_t>(initial_phase);
+  QV::add_y_phase(num_y, phase);
+
+  const uint_t mask_u = ~MASKS[x_max + 1];
+  const uint_t mask_l = MASKS[x_max];
+  auto lambda = [&](const int_t i, double &val_re, double &val_im)->void {
+    (void)val_im; // unused
+    auto idx_vec = ((i << 1) & mask_u) | (i & mask_l);
+    auto idx_mat = idx_vec ^ x_mask + nrows * idx_vec;
+    // Since rho is hermitian rho[i, j] + rho[j, i] = 2 real(rho[i, j])
+    auto val = 2 * std::real(phase * BaseVector::data_[idx_mat]);
+    if (z_mask && (AER::Utils::popcount(idx_vec & z_mask) & 1)) {
+      val = - val;
+    }
+    val_re += val;
+  };
+  return std::real(BaseVector::apply_reduction_lambda(
+    std::move(lambda), size_t(0), nrows >> 1));
+}
+
 //-----------------------------------------------------------------------
 // Z-measurement outcome probabilities
 //-----------------------------------------------------------------------

src/simulators/density_matrix/densitymatrix_thrust.hpp

@@ -119,7 +119,6 @@ class DensityMatrixThrust : public UnitaryMatrixThrust<data_t> {
   // Apply a 3-qubit toffoli gate
   void apply_toffoli(const uint_t qctrl0, const uint_t qctrl1, const uint_t qtrgt);
 
-
   //-----------------------------------------------------------------------
   // Z-measurement outcome probabilities
   //-----------------------------------------------------------------------

src/simulators/statevector/chunk/chunk_container.hpp

@@ -16,34 +16,15 @@
 #ifndef _qv_chunk_container_hpp_
 #define _qv_chunk_container_hpp_
 
+#include "misc/warnings.hpp"
+DISABLE_WARNING_PUSH
 #ifdef AER_THRUST_CUDA
 #include <cuda.h>
 #include <cuda_runtime.h>
 #endif
+DISABLE_WARNING_POP
 
-#include <thrust/for_each.h>
-#include <thrust/complex.h>
-#include <thrust/inner_product.h>
-#include <thrust/transform.h>
-#include <thrust/transform_scan.h>
-#include <thrust/binary_search.h>
-#include <thrust/execution_policy.h>
-#include <thrust/functional.h>
-#include <thrust/tuple.h>
-#include <thrust/iterator/constant_iterator.h>
-#include <thrust/detail/vector_base.h>
-
-#include <thrust/iterator/counting_iterator.h>
-#include <thrust/iterator/transform_iterator.h>
-#include <thrust/iterator/permutation_iterator.h>
-#include <thrust/fill.h>
-
-#ifdef AER_THRUST_CUDA
-#include <thrust/device_vector.h>
-#endif
-#include <thrust/host_vector.h>
-
-#include <thrust/system/omp/execution_policy.h>
+#include "misc/wrap_thrust.hpp"
 
 #include <algorithm>
 #include <array>

src/simulators/statevector/qubitvector.hpp

@@ -308,14 +308,12 @@ class QubitVector {
 
   // Return the expectation value of an N-qubit Pauli matrix.
   // The Pauli is input as a length N string of I,X,Y,Z characters.
-  double expval_pauli(const reg_t &qubits, const std::string &pauli,
-                      const complex_t &coeff = 1) const;
+  double expval_pauli(const reg_t &qubits, const std::string &pauli,const complex_t initial_phase=1.0) const;
   //for multi-chunk inter chunk expectation
   double expval_pauli(const reg_t &qubits, const std::string &pauli,
                       const QubitVector<data_t>& pair_chunk, 
                       const uint_t z_count,
-                      const uint_t z_count_pair,
-                      const complex_t &coeff = 1) const;
+                      const uint_t z_count_pair,const complex_t initial_phase=1.0) const;
 
   //-----------------------------------------------------------------------
   // JSON configuration settings
@@ -1966,8 +1964,7 @@ void add_y_phase(uint_t num_y, std::complex<data_t>& coeff){
 
 template <typename data_t>
 double QubitVector<data_t>::expval_pauli(const reg_t &qubits,
-                                         const std::string &pauli,
-                                         const complex_t &coeff) const {
+                                         const std::string &pauli,const complex_t initial_phase) const {
 
   uint_t x_mask, z_mask, num_y, x_max;
   std::tie(x_mask, z_mask, num_y, x_max) = pauli_masks_and_phase(qubits, pauli);
@@ -1976,7 +1973,7 @@ double QubitVector<data_t>::expval_pauli(const reg_t &qubits,
   if (x_mask + z_mask == 0) {
     return norm();
   }
-  auto phase = std::complex<data_t>(coeff);
+  auto phase = std::complex<data_t>(initial_phase);
   add_y_phase(num_y, phase);
 
   // specialize x_max == 0
@@ -2018,14 +2015,13 @@ template <typename data_t>
 double QubitVector<data_t>::expval_pauli(const reg_t &qubits,
                                          const std::string &pauli,
                                          const QubitVector<data_t>& pair_chunk,
-                                         const uint_t z_count,const uint_t z_count_pair,
-                                         const complex_t &coeff) const 
+                                         const uint_t z_count,const uint_t z_count_pair,const complex_t initial_phase) const 
 {
 
   uint_t x_mask, z_mask, num_y, x_max;
   std::tie(x_mask, z_mask, num_y, x_max) = pauli_masks_and_phase(qubits, pauli);
 
-  auto phase = std::complex<data_t>(coeff);
+  auto phase = std::complex<data_t>(initial_phase);
   add_y_phase(num_y, phase);
 
   std::complex<data_t>* pair_ptr = pair_chunk.data();

src/simulators/statevector/qubitvector_thrust.hpp

@@ -313,13 +313,11 @@ class QubitVectorThrust {
 
   // Return the expectation value of an N-qubit Pauli matrix.
   // The Pauli is input as a length N string of I,X,Y,Z characters.
-  double expval_pauli(const reg_t &qubits, const std::string &pauli,
-                      const complex_t &coeff = 1) const;
+  double expval_pauli(const reg_t &qubits, const std::string &pauli,const complex_t initial_phase=1.0) const;
   //for multi-chunk inter chunk expectation
   double expval_pauli(const reg_t &qubits, const std::string &pauli,
                       const QubitVectorThrust<data_t>& pair_chunk,
-                      const uint_t z_count,const uint_t z_count_pair,
-                      const complex_t &coeff = 1) const;
+                      const uint_t z_count,const uint_t z_count_pair,const complex_t initial_phase=1.0) const;
 
 
   //-----------------------------------------------------------------------
@@ -3455,12 +3453,11 @@ class expval_pauli_Z_func : public GateFuncBase<data_t>
 {
 protected:
   uint_t z_mask_;
-  thrust::complex<data_t> phase_;
+
 public:
-  expval_pauli_Z_func(uint_t z,std::complex<data_t> p)
+  expval_pauli_Z_func(uint_t z)
   {
     z_mask_ = z;
-    phase_ = p;
   }
 
   bool is_diagonal(void)
@@ -3477,7 +3474,6 @@ class expval_pauli_Z_func : public GateFuncBase<data_t>
     vec = this->data_;
 
     q0 = vec[i];
-    q0 = phase_ * q0;
     ret = q0.real()*q0.real() + q0.imag()*q0.imag();
 
     if(z_mask_ != 0){
@@ -3559,8 +3555,7 @@ class expval_pauli_XYZ_func : public GateFuncBase<data_t>
 
 template <typename data_t>
 double QubitVectorThrust<data_t>::expval_pauli(const reg_t &qubits,
-                                               const std::string &pauli,
-                                               const complex_t &coeff) const 
+                                               const std::string &pauli,const complex_t initial_phase) const 
 {
   uint_t x_mask, z_mask, num_y, x_max;
   std::tie(x_mask, z_mask, num_y, x_max) = pauli_masks_and_phase(qubits, pauli);
@@ -3569,17 +3564,16 @@ double QubitVectorThrust<data_t>::expval_pauli(const reg_t &qubits,
   if (x_mask + z_mask == 0) {
     return norm();
   }
-
-  // Compute the overall phase of the operator.
-  // This is (-1j) ** number of Y terms modulo 4
-  auto phase = std::complex<data_t>(coeff);
-  add_y_phase(num_y, phase);
-
+
   // specialize x_max == 0
   if(x_mask == 0) {
-    return apply_function_sum( expval_pauli_Z_func<data_t>(z_mask, phase) );
+    return apply_function_sum( expval_pauli_Z_func<data_t>(z_mask) );
   }
 
+  // Compute the overall phase of the operator.
+  // This is (-1j) ** number of Y terms modulo 4
+  auto phase = std::complex<data_t>(initial_phase);
+  add_y_phase(num_y, phase);
   return apply_function_sum( expval_pauli_XYZ_func<data_t>(x_mask, z_mask, x_max, phase) );
 }
 
@@ -3651,8 +3645,7 @@ template <typename data_t>
 double QubitVectorThrust<data_t>::expval_pauli(const reg_t &qubits,
                                                const std::string &pauli,
                                                const QubitVectorThrust<data_t>& pair_chunk,
-                                               const uint_t z_count,const uint_t z_count_pair,
-                                               const complex_t &coeff) const 
+                                               const uint_t z_count,const uint_t z_count_pair,const complex_t initial_phase) const 
 {
   uint_t x_mask, z_mask, num_y, x_max;
   std::tie(x_mask, z_mask, num_y, x_max) = pauli_masks_and_phase(qubits, pauli);
@@ -3705,7 +3698,7 @@ double QubitVectorThrust<data_t>::expval_pauli(const reg_t &qubits,
 
   // Compute the overall phase of the operator.
   // This is (-1j) ** number of Y terms modulo 4
-  auto phase = std::complex<data_t>(coeff);
+  auto phase = std::complex<data_t>(initial_phase);
   add_y_phase(num_y, phase);
 
   ret = apply_function_sum( expval_pauli_inter_chunk_func<data_t>(x_mask, z_mask, phase, pair_ptr,z_count,z_count_pair) );

src/simulators/statevector/statevector_state_chunk.hpp

@@ -660,15 +660,15 @@ double State<statevec_t>::expval_pauli(const reg_t &qubits,
   }
 
   if(qubits_out_chunk.size() > 0){  //there are bits out of chunk
-    std::complex<double> coeff = 1.0;
+    std::complex<double> phase = 1.0;
 
     std::reverse(pauli_out_chunk.begin(),pauli_out_chunk.end());
     std::reverse(pauli_in_chunk.begin(),pauli_in_chunk.end());
 
     uint_t x_mask, z_mask, num_y, x_max;
     std::tie(x_mask, z_mask, num_y, x_max) = AER::QV::pauli_masks_and_phase(qubits_out_chunk, pauli_out_chunk);
 
-    AER::QV::add_y_phase(num_y,coeff);
+    AER::QV::add_y_phase(num_y,phase);
 
     if(x_mask != 0){    //pairing state is out of chunk
       bool on_same_process = true;
@@ -707,12 +707,12 @@ double State<statevec_t>::expval_pauli(const reg_t &qubits,
           z_count_pair = AER::Utils::popcount(pair_chunk & z_mask);
 
           if(iProc == BaseState::distributed_rank_){  //pair is on the same process
-            expval += BaseState::qregs_[iChunk-BaseState::global_chunk_index_].expval_pauli(qubits_in_chunk, pauli_in_chunk,BaseState::qregs_[pair_chunk - BaseState::global_chunk_index_],z_count,z_count_pair,coeff);
+            expval += BaseState::qregs_[iChunk-BaseState::global_chunk_index_].expval_pauli(qubits_in_chunk, pauli_in_chunk,BaseState::qregs_[pair_chunk - BaseState::global_chunk_index_],z_count,z_count_pair);
           }
           else{
             BaseState::recv_chunk(iChunk-BaseState::global_chunk_index_,pair_chunk);
             //refer receive buffer to calculate expectation value
-            expval += BaseState::qregs_[iChunk-BaseState::global_chunk_index_].expval_pauli(qubits_in_chunk, pauli_in_chunk,BaseState::qregs_[iChunk-BaseState::global_chunk_index_],z_count,z_count_pair,coeff);
+            expval += BaseState::qregs_[iChunk-BaseState::global_chunk_index_].expval_pauli(qubits_in_chunk, pauli_in_chunk,BaseState::qregs_[iChunk-BaseState::global_chunk_index_],z_count,z_count_pair);
           }
         }
         else if(iProc == BaseState::distributed_rank_){  //pair is on this process
@@ -727,7 +727,7 @@ double State<statevec_t>::expval_pauli(const reg_t &qubits,
         double sign = 1.0;
         if (z_mask && (AER::Utils::popcount((i + BaseState::global_chunk_index_) & z_mask) & 1))
           sign = -1.0;
-        expval += sign * BaseState::qregs_[i].expval_pauli(qubits_in_chunk, pauli_in_chunk,coeff);
+        expval += sign * BaseState::qregs_[i].expval_pauli(qubits_in_chunk, pauli_in_chunk);
       }
     }
   }

test/terra/backends/qasm_simulator/qasm_save_expval.py

@@ -16,11 +16,15 @@
 from ddt import ddt, data
 from numpy import allclose
 import qiskit.quantum_info as qi
+from qiskit import QuantumCircuit
 from qiskit.circuit.library import QuantumVolume
 from qiskit.compiler import transpile, assemble
 
 from qiskit.providers.aer import QasmSimulator
 
+PAULI2 = ['II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
+          'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ']
+
 
 @ddt
 class QasmSaveExpectationValueTests:
@@ -29,8 +33,7 @@ class QasmSaveExpectationValueTests:
     SIMULATOR = QasmSimulator()
     BACKEND_OPTS = {}
 
-    @data('II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
-          'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ')
+    @data(*PAULI2)
     def test_save_expval_stabilizer_pauli(self, pauli):
         """Test Pauli expval for stabilizer circuit"""
 
@@ -62,8 +65,7 @@ def test_save_expval_stabilizer_pauli(self, pauli):
             value = result.data(0)['expval']
             self.assertAlmostEqual(value, target)
 
-    @data('II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
-          'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ')
+    @data(*PAULI2)
     def test_save_expval_var_stabilizer_pauli(self, pauli):
         """Test Pauli expval_var for stabilizer circuit"""
 
@@ -163,8 +165,7 @@ def test_save_expval_var_stabilizer_hermitian(self, qubits):
             value = result.data(0)['expval']
             self.assertTrue(allclose(value, target))
 
-    @data('II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
-          'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ')
+    @data(*PAULI2)
     def test_save_expval_nonstabilizer_pauli(self, pauli):
         """Test Pauli expval for non-stabilizer circuit"""
 
@@ -195,8 +196,7 @@ def test_save_expval_nonstabilizer_pauli(self, pauli):
             value = result.data(0)['expval']
             self.assertAlmostEqual(value, target)
 
-    @data('II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
-          'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ')
+    @data(*PAULI2)
     def test_save_expval_var_nonstabilizer_pauli(self, pauli):
         """Test Pauli expval_var for non-stabilizer circuit"""
 
@@ -292,3 +292,69 @@ def test_save_expval_var_nonstabilizer_hermitian(self, qubits):
             self.assertTrue(result.success)
             value = result.data(0)['expval']
             self.assertTrue(allclose(value, target))
+
+    @data(*PAULI2)
+    def test_save_expval_cptp_pauli(self, pauli):
+        """Test Pauli expval for stabilizer circuit"""
+
+        SUPPORTED_METHODS = [
+            'density_matrix', 'density_matrix_gpu', 'density_matrix_thrust'
+        ]
+        SEED = 5832
+
+        opts = self.BACKEND_OPTS.copy()
+        if opts.get('method') in SUPPORTED_METHODS:
+
+            oper = qi.Pauli(pauli)
+
+            # CPTP channel test circuit
+            channel = qi.random_quantum_channel(4, seed=SEED)
+            state_circ = QuantumCircuit(2)
+            state_circ.append(channel, range(2))
+
+            state = qi.DensityMatrix(state_circ)
+            target = state.expectation_value(oper).real.round(10)
+
+            # Snapshot circuit
+            circ = transpile(state_circ, self.SIMULATOR)
+            circ.save_expectation_value(oper, [0, 1], label='expval')
+            qobj = assemble(circ)
+            result = self.SIMULATOR.run(qobj, **opts).result()
+
+            self.assertTrue(result.success)
+            value = result.data(0)['expval']
+            self.assertAlmostEqual(value, target)
+
+    @data(*PAULI2)
+    def test_save_expval_var_cptp_pauli(self, pauli):
+        """Test Pauli expval_var for stabilizer circuit"""
+
+        SUPPORTED_METHODS = [
+            'density_matrix', 'density_matrix_gpu', 'density_matrix_thrust'
+        ]
+        SEED = 5832
+
+        opts = self.BACKEND_OPTS.copy()
+        if opts.get('method') in SUPPORTED_METHODS:
+
+            oper = qi.Pauli(pauli)
+
+            # CPTP channel test circuit
+            channel = qi.random_quantum_channel(4, seed=SEED)
+            state_circ = QuantumCircuit(2)
+            state_circ.append(channel, range(2))
+
+            state = qi.DensityMatrix(state_circ)
+            expval = state.expectation_value(oper).real
+            variance = state.expectation_value(oper ** 2).real - expval ** 2
+            target = [expval, variance]
+
+            # Snapshot circuit
+            circ = transpile(state_circ, self.SIMULATOR)
+            circ.save_expectation_value_variance(oper, [0, 1], label='expval')
+            qobj = assemble(circ)
+            result = self.SIMULATOR.run(qobj, **opts).result()
+
+            self.assertTrue(result.success)
+            value = result.data(0)['expval']
+            self.assertTrue(allclose(value, target))