PR #2683 Follow up

.github/workflows/publishing.yml

@@ -839,11 +839,10 @@ jobs:
       - name: Basic validation (GPU backends)
         shell: bash
         run: |
-          # We don't test the 'dynamics' target in these tests
           backends_to_test=`\
               for file in $(ls $CUDA_QUANTUM_PATH/targets/*.yml); \
               do
-                  if [[ "$file" != *"dynamics"* ]] && [ -n "$(cat $file | grep "gpu-requirements")" ]; then \
+                  if [ -n "$(cat $file | grep "gpu-requirements")" ]; then \
                     basename $file | cut -d "." -f 1; \
                   elif [ -n "$(basename $file | grep mqpu)" ]; then \
                     echo remote-mqpu; \

docs/sphinx/examples/cpp/dynamics/cavity_qed.cpp

@@ -108,7 +108,7 @@ int main() {
   auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
     std::vector<double> expectations;
 
-    auto all_exps = result.get_expectation_values().value();
+    auto all_exps = result.expectation_values.value();
     for (auto exp_vals : all_exps) {
       expectations.push_back((double)exp_vals[idx]);
     }

docs/sphinx/examples/cpp/dynamics/cross_resonance.cpp

@@ -119,7 +119,7 @@ int main() {
   auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
     std::vector<double> expectations;
 
-    auto all_exps = result.get_expectation_values().value();
+    auto all_exps = result.expectation_values.value();
     for (auto exp_vals : all_exps) {
       expectations.push_back((double)exp_vals[idx]);
     }

docs/sphinx/examples/cpp/dynamics/heisenberg_model.cpp

@@ -121,7 +121,7 @@ int main() {
     auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
       std::vector<double> expectations;
 
-      auto all_exps = result.get_expectation_values().value();
+      auto all_exps = result.expectation_values.value();
       for (auto exp_vals : all_exps) {
         expectations.push_back((double)exp_vals[idx]);
       }

docs/sphinx/examples/cpp/dynamics/qubit_control.cpp

@@ -46,9 +46,8 @@ int main() {
   // 2. A time-dependent driving term: omega_x * cos(omega_drive * t) * `Sx`_0,
   // which induces rotations about the X-axis. The scalar_operator(mod_`func`)
   // allows the drive term to vary in time according to mod_`func`.
-  auto hamiltonian =
-      0.5 * omega_z * cudaq::spin_operator::z(0) +
-      omega_x * cudaq::scalar_operator(mod_func) * cudaq::spin_operator::x(0);
+  auto hamiltonian = 0.5 * omega_z * cudaq::spin_operator::z(0) +
+                     mod_func * cudaq::spin_operator::x(0) * omega_x;
 
   // A single qubit with dimension 2.
   cudaq::dimension_map dimensions = {{0, 2}};
@@ -100,7 +99,7 @@ int main() {
   auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
     std::vector<double> expectations;
 
-    auto all_exps = result.get_expectation_values().value();
+    auto all_exps = result.expectation_values.value();
     for (auto exp_vals : all_exps) {
       expectations.push_back(exp_vals[idx].expectation());
     }

docs/sphinx/examples/cpp/dynamics/qubit_dynamics.cpp

@@ -54,7 +54,7 @@ int main() {
   auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
     std::vector<double> expectations;
 
-    auto all_exps = result.get_expectation_values().value();
+    auto all_exps = result.expectation_values.value();
     for (auto exp_vals : all_exps) {
       expectations.push_back((double)exp_vals[idx]);
     }

python/runtime/common/py_EvolveResult.cpp

@@ -42,17 +42,15 @@ void bindEvolveResult(py::module &mod) {
           "detail.\n")
       .def(
           "intermediate_states",
-          [](evolve_result &self) { return self.get_intermediate_states(); },
+          [](evolve_result &self) { return self.intermediate_states; },
           "Stores all intermediate states, meaning the state after each step "
           "in a defined schedule, produced by a call to :func:`evolve`, "
           "including the final state. This property is only populated if "
           "saving intermediate results was requested in the call to "
           ":func:`evolve`.\n")
       .def(
           "final_expectation_values",
-          [](evolve_result &self) {
-            return self.get_final_expectation_values();
-          },
+          [](evolve_result &self) { return self.final_expectation_values; },
           "Stores the final expectation values, that is the results produced "
           "by "
           "calls to :func:`observe`, triggered by a call to :func:`evolve`. "
@@ -63,7 +61,7 @@ void bindEvolveResult(py::module &mod) {
           "call.\n")
       .def(
           "expectation_values",
-          [](evolve_result &self) { return self.get_expectation_values(); },
+          [](evolve_result &self) { return self.expectation_values; },
           "Stores the expectation values, that is the results from the calls "
           "to "
           ":func:`observe`, at each step in the schedule produced by a call to "

runtime/common/EvolveResult.h

@@ -24,20 +24,19 @@ class evolve_result {
   // The final state after the time evolution.
   state final_state;
 
+public:
   // The state after each step in the evolution, if computed.
-  std::optional<std::vector<state>> intermediate_states = {};
+  std::optional<std::vector<state>> intermediate_states;
 
   // The computed expectation values at the end of the evolution.
-  std::optional<std::vector<observe_result>> final_expectation_values = {};
+  std::optional<std::vector<observe_result>> final_expectation_values;
 
   // The computed expectation values for each step of the evolution.
-  std::optional<std::vector<std::vector<observe_result>>> expectation_values =
-      {};
+  std::optional<std::vector<std::vector<observe_result>>> expectation_values;
 
   // The result of sampling of an analog Hamiltonian simulation on a QPU
-  std::optional<sample_result> sampling_result = {};
+  std::optional<sample_result> sampling_result;
 
-public:
   evolve_result(state state) : final_state(state) {}
 
   evolve_result(state state, const std::vector<observe_result> &expectations)
@@ -100,24 +99,6 @@ class evolve_result {
     return final_state;
   }
 
-  const std::optional<std::vector<state>> &get_intermediate_states() {
-    return intermediate_states;
-  }
-
-  const std::optional<std::vector<observe_result>> &
-  get_final_expectation_values() {
-    return final_expectation_values;
-  }
-
-  const std::optional<std::vector<std::vector<observe_result>>> &
-  get_expectation_values() {
-    return expectation_values;
-  }
-
-  const std::optional<sample_result> &get_sampling_result() {
-    return sampling_result;
-  }
-
 private:
   state getLastStateIfValid(const std::vector<state> &states) {
     if (states.empty())

runtime/cudaq/dynamics/schedule.cpp

@@ -71,9 +71,9 @@ void schedule::reset() {
 
 // Get the current step
 std::optional<std::complex<double>> schedule::current_step() const {
-  if (current_idx >= 0 && current_idx < static_cast<int>(steps.size())) {
+  if (current_idx >= 0 && current_idx < static_cast<int>(steps.size()))
     return steps[current_idx];
-  }
+
   return std::nullopt;
 }
 

runtime/nvqir/cudensitymat/CuDensityMatEvolution.cpp

@@ -59,9 +59,8 @@ evolve_result evolveSingle(
   cudmState->initialize_cudm(handle, dims);
 
   state initial_State = [&]() {
-    if (!collapseOperators.empty() && !cudmState->is_density_matrix()) {
+    if (!collapseOperators.empty() && !cudmState->is_density_matrix())
       return state(new CuDensityMatState(cudmState->to_density_matrix()));
-    }
     return initialState;
   }();
 

runtime/nvqir/cudensitymat/CuDensityMatOpConverter.cpp

@@ -118,9 +118,8 @@ cudaq::product_operator<cudaq::matrix_operator> computeDagger(
   }
   std::reverse(daggerOps.begin(), daggerOps.end());
 
-  if (daggerOps.empty()) {
+  if (daggerOps.empty())
     throw std::runtime_error("Empty product operator");
-  }
   cudaq::product_operator<cudaq::matrix_operator> daggerProduct = daggerOps[0];
   for (std::size_t i = 1; i < daggerOps.size(); ++i) {
     daggerProduct *= daggerOps[i];
@@ -346,9 +345,9 @@ cudaq::dynamics::CuDensityMatOpConverter::createProductOperatorTerm(
       int degree = sub_degrees[j];
       int modality = modalities[j];
 
-      if (sub_degrees[i] < 0) {
+      if (sub_degrees[i] < 0)
         throw std::out_of_range("Degree cannot be negative!");
-      }
+
       allDegrees.emplace_back(degree);
       allModeActionDuality.emplace_back(modality);
     }
@@ -367,9 +366,8 @@ cudaq::dynamics::CuDensityMatOpConverter::convertToCudensitymatOperator(
     const std::unordered_map<std::string, std::complex<double>> &parameters,
     const operator_sum<cudaq::matrix_operator> &op,
     const std::vector<int64_t> &modeExtents) {
-  if (op.get_terms().empty()) {
+  if (op.get_terms().empty())
     throw std::invalid_argument("Operator sum cannot be empty.");
-  }
 
   cudensitymatOperator_t cudmOperator;
   HANDLE_CUDM_ERROR(cudensitymatCreateOperator(

runtime/nvqir/cudensitymat/CuDensityMatState.cpp

@@ -20,10 +20,9 @@ CuDensityMatState::overlap(const cudaq::SimulationState &other) {
     throw std::runtime_error("[CuDensityMatState] overlap error - other state "
                              "dimension not equal to this state dimension.");
 
-  if (other.getPrecision() != getPrecision()) {
+  if (other.getPrecision() != getPrecision())
     throw std::runtime_error(
         "[CuDensityMatState] overlap error - precision mismatch.");
-  }
 
   if (!isDensityMatrix) {
     Eigen::VectorXcd state(dimension);
@@ -112,10 +111,10 @@ CuDensityMatState::createFromSizeAndPtr(std::size_t size, void *dataPtr,
 // for an invalid tensor index.
 cudaq::SimulationState::Tensor
 CuDensityMatState::getTensor(std::size_t tensorIdx) const {
-  if (tensorIdx != 0) {
+  if (tensorIdx != 0)
     throw std::runtime_error(
         "CuDensityMatState state only supports a single tensor");
-  }
+
   const std::vector<std::size_t> extents =
       isDensityMatrix ? std::vector<std::size_t>{dimension, dimension}
                       : std::vector<std::size_t>{dimension};
@@ -156,10 +155,10 @@ CuDensityMatState::operator()(std::size_t tensorIdx,
 // Copy the state device data to the user-provided host data pointer.
 void CuDensityMatState::toHost(std::complex<double> *userData,
                                std::size_t numElements) const {
-  if (numElements != dimension * (isDensityMatrix ? dimension : 1)) {
+  if (numElements != dimension * (isDensityMatrix ? dimension : 1))
     throw std::runtime_error("Number of elements in user data does not match "
                              "the size of the state");
-  }
+
   HANDLE_CUDA_ERROR(cudaMemcpy(userData, devicePtr,
                                numElements * sizeof(std::complex<double>),
                                cudaMemcpyDeviceToHost));
@@ -205,9 +204,9 @@ CuDensityMatState::CuDensityMatState(
     const std::vector<int64_t> &dims)
     : cudmHandle(handle), dimension(rawData.size()), cudmState(nullptr),
       hilbertSpaceDims(dims) {
-  if (rawData.empty()) {
+  if (rawData.empty())
     throw std::invalid_argument("Raw data cannot be empty.");
-  }
+
   // Allocate device memory
   size_t dataSize = rawData.size() * sizeof(std::complex<double>);
   HANDLE_CUDA_ERROR(
@@ -255,10 +254,9 @@ CuDensityMatState::CuDensityMatState(
   // Validate device memory
   size_t totalSize = std::accumulate(componentBufferSizes.begin(),
                                      componentBufferSizes.end(), 0);
-  if (totalSize > rawData.size() * sizeof(std::complex<double>)) {
+  if (totalSize > rawData.size() * sizeof(std::complex<double>))
     throw std::invalid_argument(
         "Device memory size is insufficient to cover all components.");
-  }
 
   // Attach storage for using device memory (devicePtr)
   std::vector<void *> componentBuffers(numStateComponents);
@@ -426,12 +424,11 @@ CuDensityMatState &
 CuDensityMatState::operator=(CuDensityMatState &&other) noexcept {
   if (this != &other) {
     // Free existing resources
-    if (cudmState) {
+    if (cudmState)
       cudensitymatDestroyState(cudmState);
-    }
-    if (devicePtr) {
+
+    if (devicePtr)
       cudaFree(devicePtr);
-    }
 
     // Move data from other
     isDensityMatrix = other.isDensityMatrix;
@@ -456,21 +453,18 @@ CuDensityMatState::~CuDensityMatState() { destroyState(); }
 bool CuDensityMatState::is_initialized() const { return cudmState != nullptr; }
 
 bool cudaq::CuDensityMatState::is_density_matrix() const {
-  if (!is_initialized()) {
+  if (!is_initialized())
     return false;
-  }
 
   return dimension == calculate_density_matrix_size(hilbertSpaceDims);
 }
 
 CuDensityMatState cudaq::CuDensityMatState::to_density_matrix() const {
-  if (!is_initialized()) {
+  if (!is_initialized())
     throw std::runtime_error("State is not initialized.");
-  }
 
-  if (is_density_matrix()) {
+  if (is_density_matrix())
     throw std::runtime_error("State is already a density matrix.");
-  }
 
   size_t vectorSize = calculate_state_vector_size(hilbertSpaceDims);
   std::vector<std::complex<double>> stateVecData(vectorSize);
@@ -543,9 +537,8 @@ void CuDensityMatState::initialize_cudm(cudensitymatHandle_t handleToSet,
 
 CuDensityMatState
 cudaq::CuDensityMatState::operator+(const CuDensityMatState &other) const {
-  if (dimension != other.dimension) {
+  if (dimension != other.dimension)
     throw std::invalid_argument("State size mismatch for addition.");
-  }
 
   CuDensityMatState result = CuDensityMatState::clone(*this);
 
@@ -556,11 +549,10 @@ cudaq::CuDensityMatState::operator+(const CuDensityMatState &other) const {
 
 CuDensityMatState &
 cudaq::CuDensityMatState::operator+=(const CuDensityMatState &other) {
-  if (dimension != other.dimension) {
+  if (dimension != other.dimension)
     throw std::invalid_argument(
         fmt::format("State size mismatch for addition ({} vs {}).", dimension,
                     other.dimension));
-  }
 
   cuDoubleComplex scalar{1.0, 0.0};
   HANDLE_CUBLAS_ERROR(cublasZaxpy(

runtime/nvqir/cudensitymat/CuDensityMatTimeStepper.cpp

@@ -17,9 +17,8 @@ CuDensityMatTimeStepper::CuDensityMatTimeStepper(
 state CuDensityMatTimeStepper::compute(
     const state &inputState, double t, double step_size,
     const std::unordered_map<std::string, std::complex<double>> &parameters) {
-  if (step_size == 0.0) {
+  if (step_size == 0.0)
     throw std::runtime_error("Step size cannot be zero.");
-  }
 
   auto *simState =
       cudaq::state_helper::getSimulationState(const_cast<state *>(&inputState));
@@ -34,15 +33,13 @@ state CuDensityMatTimeStepper::compute(
   // Create a new state for the next step
   auto next_state = CuDensityMatState::zero_like(state);
 
-  if (!next_state.is_initialized()) {
+  if (!next_state.is_initialized())
     throw std::runtime_error("Next state failed to initialize.");
-  }
 
-  if (state.get_hilbert_space_dims() != next_state.get_hilbert_space_dims()) {
+  if (state.get_hilbert_space_dims() != next_state.get_hilbert_space_dims())
     throw std::runtime_error("As the dimensions of both the old and the new "
                              "state do no match, the "
                              "operator cannot act on the states.");
-  }
 
   // Prepare the operator for action
   HANDLE_CUDM_ERROR(cudensitymatOperatorPrepareAction(

runtime/nvqir/cudensitymat/CuDensityMatUtils.h

@@ -13,6 +13,7 @@
 #include <vector>
 
 namespace cudaq::dynamics {
+// Adapted from cuquantum team
 // GPU memory management
 template <std::floating_point T>
 void *createArrayGpu(const std::vector<std::complex<T>> &cpuArray) {
@@ -26,6 +27,8 @@ void *createArrayGpu(const std::vector<std::complex<T>> &cpuArray) {
   }
   return gpuArray;
 }
+
+// Adapted from cuquantum team
 inline void destroyArrayGpu(void *gpuArray) {
   if (gpuArray)
     HANDLE_CUDA_ERROR(cudaFree(gpuArray));