diff --git a/docs/sphinx/targets/cpp/infleqtion.cpp b/docs/sphinx/targets/cpp/infleqtion.cpp index 4a61b71571..418252cffc 100644 --- a/docs/sphinx/targets/cpp/infleqtion.cpp +++ b/docs/sphinx/targets/cpp/infleqtion.cpp @@ -1,30 +1,30 @@ // Compile and run with: // ``` -// nvq++ --target infleqtion infleqtion.cpp -o out.x && ./out.x +// nvq++ --target `infleqtion` `infleqtion`.cpp -o out.x && ./out.x // ``` -// This will submit the job to the Infleqtion's ideal simulator, -// cq_sqale_simulator (default). Alternatively, we can enable hardware noise -// model simulation by specifying `noise-sim` to the flag `--infleqtion-method`, +// This will submit the job to the `Infleqtion's` ideal simulator, +// `cq`_`sqale`_simulator (default). Alternatively, we can enable hardware noise +// model simulation by specifying `noise-sim` to the flag --`infleqtion`-method, // e.g., // ``` -// nvq++ --target infleqtion --infleqtion-machine cq_sqale_qpu -// --infleqtion-method noise-sim infleqtion.cpp -o out.x && ./out.x +// nvq++ --target `infleqtion` --`infleqtion`-machine `cq`_`sqale`_qpu +// --`infleqtion`-method noise-sim `infleqtion`.cpp -o out.x && ./out.x // ``` -// where "noise-sim" instructs Superstaq to perform a noisy emulation of the +// where "noise-sim" instructs `Superstaq` to perform a noisy emulation of the // QPU. An ideal dry-run execution on the QPU may be performed by passing -// `dry-run` to the `--infleqtion-method` flag, e.g., +// `dry-run` to the --`infleqtion`-method flag, e.g., // ``` -// nvq++ --target infleqtion --infleqtion-machine cq_sqale_qpu -// --infleqtion-method dry-run infleqtion.cpp -o out.x && ./out.x +// nvq++ --target `infleqtion` --`infleqtion`-machine `cq`_`sqale`_qpu +// --`infleqtion`-method dry-run `infleqtion`.cpp -o out.x && ./out.x // ``` -// Note: If targeting ideal cloud simulation, `--infleqtion-machine -// cq_sqale_simulator` is optional since it is the default configuration if not -// provided. +// Note: If targeting ideal cloud simulation, --`infleqtion`-machine +// `cq`_`sqale`_simulator is optional since it is the default configuration if +// not provided. #include #include -// Define a simple quantum kernel to execute on Infleqtion backends. +// Define a simple quantum kernel to execute on `infleqtion` backends. struct ghz { // Maximally entangled state between 5 qubits. auto operator()() __qpu__ { @@ -38,7 +38,7 @@ struct ghz { }; int main() { - // Submit to infleqtion asynchronously (e.g., continue executing + // Submit to `infleqtion` asynchronously (e.g., continue executing // code in the file until the job has been returned). auto future = cudaq::sample_async(ghz{}); // ... classical code to execute in the meantime ... @@ -58,7 +58,7 @@ int main() { auto async_counts = readIn.get(); async_counts.dump(); - // OR: Submit to infleqtion synchronously (e.g., wait for the job + // OR: Submit to `infleqtion` synchronously (e.g., wait for the job // result to be returned before proceeding). auto counts = cudaq::sample(ghz{}); counts.dump(); diff --git a/runtime/cudaq/algorithms/observe.h b/runtime/cudaq/algorithms/observe.h index 407a528013..dfaa7c9341 100644 --- a/runtime/cudaq/algorithms/observe.h +++ b/runtime/cudaq/algorithms/observe.h @@ -59,9 +59,9 @@ concept ObserveCallValid = /// @param shots number of shots to run for the given kernel, or -1 if not /// applicable. /// @param noise noise model to use for the sample operation -/// @param num_trajectories is the optional number of trajectories to be used when -/// computing the expectation values in the presence of noise. This parameter is -/// only applied to simulation backends that support noisy +/// @param num_trajectories is the optional number of trajectories to be used +/// when computing the expectation values in the presence of noise. This +/// parameter is only applied to simulation backends that support noisy /// simulation of trajectories. struct observe_options { int shots = -1; diff --git a/runtime/cudaq/dynamics.h b/runtime/cudaq/dynamics.h new file mode 100644 index 0000000000..2697d17f03 --- /dev/null +++ b/runtime/cudaq/dynamics.h @@ -0,0 +1,11 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cudaq/dynamics/operators.h" diff --git a/runtime/cudaq/dynamics/CMakeLists.txt b/runtime/cudaq/dynamics/CMakeLists.txt new file mode 100644 index 0000000000..c08bbf7d23 --- /dev/null +++ b/runtime/cudaq/dynamics/CMakeLists.txt @@ -0,0 +1,55 @@ +# ============================================================================ # +# Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. # +# All rights reserved. # +# # +# This source code and the accompanying materials are made available under # +# the terms of the Apache License 2.0 which accompanies this distribution. # +# ============================================================================ # + +if(NOT CUDENSITYMAT_ROOT) + SET(CUDENSITYMAT_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}") +endif() + +set(LIBRARY_NAME cudaq-operators) +set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-ctad-maybe-unsupported") +set(INTERFACE_POSITION_INDEPENDENT_CODE ON) + +# FIXME: add boost_integrator.cpp rungakutta_integrator.cpp back to fix compile errors +set(CUDAQ_OPS_SRC + scalar_operator.cpp elementary_operator.cpp product_operator.cpp operator_sum.cpp operator_helpers.cpp matrix_arithmetics.cpp +) + +add_library(${LIBRARY_NAME} SHARED ${CUDAQ_OPS_SRC}) +set_property(GLOBAL APPEND PROPERTY CUDAQ_RUNTIME_LIBS ${LIBRARY_NAME}) +target_include_directories(${LIBRARY_NAME} + PUBLIC + $ + $ + $ + ${CUDENSITYMAT_ROOT}/include + ${CUDAToolkit_INCLUDE_DIRS} + $ + PRIVATE .) + +set(OPERATOR_DEPENDENCIES "") +list(APPEND OPERATOR_DEPENDENCIES fmt::fmt-header-only) +add_openmp_configurations(${LIBRARY_NAME} OPERATOR_DEPENDENCIES) + +find_library(CUDENSITYMAT_LIB REQUIRED + NAMES cudensitymat + HINTS + ${CUDENSITYMAT_ROOT}/lib +) + +get_filename_component(CUDENSITYMAT_LIB_DIR ${CUDENSITYMAT_LIB} DIRECTORY) +SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_RPATH}:${CUDENSITYMAT_LIB_DIR}") + +target_link_libraries(${LIBRARY_NAME} PRIVATE ${OPERATOR_DEPENDENCIES}) +target_link_libraries(${LIBRARY_NAME} PRIVATE ${CUDENSITYMAT_LIB}) + +install(TARGETS ${LIBRARY_NAME} EXPORT cudaq-operator-targets DESTINATION lib) + +install(EXPORT cudaq-operator-targets + FILE CUDAQSpinTargets.cmake + NAMESPACE cudaq:: + DESTINATION lib/cmake/cudaq) diff --git a/runtime/cudaq/dynamics/elementary_operator.cpp b/runtime/cudaq/dynamics/elementary_operator.cpp new file mode 100644 index 0000000000..ca47a3c9c9 --- /dev/null +++ b/runtime/cudaq/dynamics/elementary_operator.cpp @@ -0,0 +1,561 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "common/EigenDense.h" +#include "common/Logger.h" +#include "operators.h" +#include +#include + +namespace cudaq { + +/// Elementary Operator constructor. +elementary_operator::elementary_operator(std::string operator_id, + std::vector degrees) + : id(operator_id), degrees(degrees) {} +elementary_operator::elementary_operator(const elementary_operator &other) + : m_ops(other.m_ops), expected_dimensions(other.expected_dimensions), + degrees(other.degrees), id(other.id) {} +elementary_operator::elementary_operator(elementary_operator &other) + : m_ops(other.m_ops), expected_dimensions(other.expected_dimensions), + degrees(other.degrees), id(other.id) {} + +elementary_operator elementary_operator::identity(int degree) { + std::string op_id = "identity"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + int degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + // Build up the identity matrix. + for (std::size_t i = 0; i < dimension; i++) { + mat(i, i) = 1.0 + 0.0 * 'j'; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "done\n\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::zero(int degree) { + std::string op_id = "zero"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + // Need to set the degree via the op itself because the + // argument to the outer function goes out of scope when + // the user invokes this later on via, e.g, `to_matrix()`. + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + mat.set_zero(); + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::annihilate(int degree) { + std::string op_id = "annihilate"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + for (std::size_t i = 0; i + 1 < dimension; i++) { + mat(i, i + 1) = std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::create(int degree) { + std::string op_id = "create"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + for (std::size_t i = 0; i + 1 < dimension; i++) { + mat(i + 1, i) = std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::position(int degree) { + std::string op_id = "position"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = cudaq::complex_matrix(dimension, dimension); + // position = 0.5 * (create + annihilate) + for (std::size_t i = 0; i + 1 < dimension; i++) { + mat(i + 1, i) = 0.5 * std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + mat(i, i + 1) = 0.5 * std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::momentum(int degree) { + std::string op_id = "momentum"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = cudaq::complex_matrix(dimension, dimension); + // momentum = 0.5j * (create - annihilate) + for (std::size_t i = 0; i + 1 < dimension; i++) { + mat(i + 1, i) = + (0.5j) * std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + mat(i, i + 1) = + -1. * (0.5j) * std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::number(int degree) { + std::string op_id = "number"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + for (std::size_t i = 0; i < dimension; i++) { + mat(i, i) = static_cast(i) + 0.0j; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "done\n\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator elementary_operator::parity(int degree) { + std::string op_id = "parity"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto mat = complex_matrix(dimension, dimension); + for (std::size_t i = 0; i < dimension; i++) { + mat(i, i) = std::pow(-1., static_cast(i)) + 0.0j; + } + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "done\n\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + return op; +} + +elementary_operator +elementary_operator::displace(int degree, std::complex amplitude) { + std::string op_id = "displace"; + std::vector degrees = {degree}; + auto op = elementary_operator(op_id, degrees); + // A dimension of -1 indicates this operator can act on any dimension. + op.expected_dimensions[degree] = -1; + if (op.m_ops.find(op_id) == op.m_ops.end()) { + auto func = [&](std::map dimensions, + std::map> _none) { + auto degree = op.degrees[0]; + int dimension = dimensions[degree]; + auto temp_mat = cudaq::complex_matrix(dimension, dimension); + // displace = exp[ (amplitude * create) - (conj(amplitude) * annihilate) ] + for (std::size_t i = 0; i + 1 < dimension; i++) { + temp_mat(i + 1, i) = + amplitude * std::sqrt(static_cast(i + 1)) + 0.0 * 'j'; + temp_mat(i, i + 1) = + -1. * std::conj(amplitude) * std::sqrt(static_cast(i + 1)) + + 0.0 * 'j'; + } + // Not ideal that our method of computing the matrix exponential + // requires copies here. Maybe we can just use eigen directly here + // to limit to one copy, but we can address that later. + auto mat = temp_mat.exp(); + std::cout << "dumping the complex mat: \n"; + mat.dump(); + std::cout << "\ndone\n"; + return mat; + }; + op.define(op_id, op.expected_dimensions, func); + } + throw std::runtime_error("currently have a bug in implementation."); + return op; +} + +elementary_operator +elementary_operator::squeeze(int degree, std::complex amplitude) { + throw std::runtime_error("Not yet implemented."); +} + +complex_matrix elementary_operator::to_matrix( + const std::map dimensions, + const std::map> parameters) const { + if (m_ops.find(id) == m_ops.end()) + throw std::runtime_error( + fmt::format("No operator found with this ID '{}'", id)); + + return m_ops.at(id).generator(dimensions, parameters); +} + +template +TEval elementary_operator::_evaluate( + operator_arithmetics &arithmetics) const { + std::cout << "In ElementaryOp _evaluate" << std::endl; + return arithmetics.evaluate(); +} +/// Elementary Operator Arithmetic. + +operator_sum elementary_operator::operator+(scalar_operator other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + return operator_sum({product_operator({*this}), product_operator({other})}); +} + +operator_sum elementary_operator::operator-(scalar_operator other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + return operator_sum( + {product_operator({*this}), product_operator({-1. * other})}); +} + +product_operator elementary_operator::operator*(scalar_operator other) { + return product_operator({*this, other}); +} + +product_operator elementary_operator::operator/(scalar_operator other) { + return product_operator({*this, (1. / other)}); +} + +operator_sum elementary_operator::operator+(std::complex other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + auto other_scalar = scalar_operator(other); + return operator_sum( + {product_operator({*this}), product_operator({other_scalar})}); +} + +operator_sum elementary_operator::operator-(std::complex other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + auto other_scalar = scalar_operator((-1. * other)); + return operator_sum( + {product_operator({*this}), product_operator({other_scalar})}); +} + +product_operator elementary_operator::operator*(std::complex other) { + auto other_scalar = scalar_operator(other); + return product_operator({*this, other_scalar}); +} + +product_operator elementary_operator::operator/(std::complex other) { + auto other_scalar = scalar_operator((1. / other)); + return product_operator({*this, other_scalar}); +} + +operator_sum elementary_operator::operator+(double other) { + std::complex value(other, 0.0); + return *this + value; +} + +operator_sum elementary_operator::operator-(double other) { + std::complex value(other, 0.0); + return *this - value; +} + +product_operator elementary_operator::operator*(double other) { + std::complex value(other, 0.0); + return *this * value; +} + +product_operator elementary_operator::operator/(double other) { + std::complex value(other, 0.0); + return *this / value; +} + +operator_sum operator+(std::complex other, elementary_operator self) { + auto other_scalar = scalar_operator(other); + return operator_sum( + {product_operator({other_scalar}), product_operator({self})}); +} + +operator_sum operator-(std::complex other, elementary_operator self) { + auto other_scalar = scalar_operator(other); + return operator_sum({product_operator({other_scalar}), (-1. * self)}); +} + +product_operator operator*(std::complex other, + elementary_operator self) { + auto other_scalar = scalar_operator(other); + return product_operator({other_scalar, self}); +} + +product_operator operator/(std::complex other, + elementary_operator self) { + auto other_scalar = scalar_operator((1. / other)); + return product_operator({other_scalar, self}); +} + +operator_sum operator+(double other, elementary_operator self) { + auto other_scalar = scalar_operator(other); + return operator_sum( + {product_operator({other_scalar}), product_operator({self})}); +} + +operator_sum operator-(double other, elementary_operator self) { + auto other_scalar = scalar_operator(other); + return operator_sum({product_operator({other_scalar}), (-1. * self)}); +} + +product_operator operator*(double other, elementary_operator self) { + auto other_scalar = scalar_operator(other); + return product_operator({other_scalar, self}); +} + +product_operator operator/(double other, elementary_operator self) { + auto other_scalar = scalar_operator((1. / other)); + return product_operator({other_scalar, self}); +} + +product_operator elementary_operator::operator*(elementary_operator other) { + return product_operator({*this, other}); +} + +operator_sum elementary_operator::operator+(elementary_operator other) { + return operator_sum({product_operator({*this}), product_operator({other})}); +} + +operator_sum elementary_operator::operator-(elementary_operator other) { + return operator_sum({product_operator({*this}), (-1. * other)}); +} + +/// FIXME: +// product_operator elementary_operator::operator/(elementary_operator other) { +// return product_operator({*this, (1./other)}); +// } + +// /// temporary helper +// std::vector reverseArange(int start, int stop, int step = 1) { +// std::vector values; +// if (start == stop) { +// values = {start}; +// } else { +// for (int value = stop; value > start; value -= step) +// values.push_back(value); +// } +// return values; +// } + +// /// Temporary helper function. +// /// FIXME: Needs a bunch of performance improvement but the initial goal +// /// is just to get something working. +// complex_matrix _cast_to_full_hilbert(complex_matrix &mat, int degree, +// std::vector allDegreesReverseOrder, +// std::map dimensions) { +// std::cout << "\nL343\n"; +// std::vector degreesSubset; +// for (auto degree_index : allDegreesReverseOrder) { +// // If we don't have a dimension defined for this degree of freedom, +// // skip it. +// if (dimensions.find(degree_index) == dimensions.end()) { +// continue; +// } else { +// degreesSubset.push_back(degree_index); +// } +// } +// std::cout << "\nL354\n"; + +// // Create an identity matrix appropriately sized for the initial degree of +// // freedom in our reveresed subset. +// complex_matrix returnMatrix = complex_matrix::identity( +// dimensions[degreesSubset[0]], dimensions[degreesSubset[0]]); +// // Now we can just loop through the degrees of freedom that we have levels +// // defined for. +// // If the initial degree of freedom is for the current `mat`, then we will +// // start with that instead of an identity matrix. +// if (degreesSubset[0] == degree) { +// returnMatrix = mat; +// } +// for (auto degree_index : degreesSubset) { +// if (degree != degree_index) { +// auto identity = complex_matrix::identity(dimensions[degree_index], +// dimensions[degree_index]); +// returnMatrix = returnMatrix.kronecker(identity); +// } +// // If this is the degree of our operator, no identity padding is needed. +// else { +// continue; +// } +// } +// std::cout << "\nL379\n"; +// return returnMatrix; +// } + +/// FIXME: Assuming each operator only acts on a single degree of freedom +/// for right now. +// product_operator elementary_operator::operator*(elementary_operator other) { +// /// FIXME: See above about assuming only 1 DOF per operator for right now. +// int thisDegree = this->degrees[0]; +// int otherDegree = other.degrees[0]; +// auto allDegrees = reverseArange(std::min(thisDegree, otherDegree), +// std::max(thisDegree, otherDegree)); + +// auto returnOperator = product_operator({*this, other}); + +// if (thisDegree == otherDegree) { +// auto func = [&](std::map dimensions, +// std::map> parameters) { +// /// FIXME: Making everything its own explicit variable for now while +// /// I work out testing bugs. +// auto thisMatrix = this->to_matrix(dimensions, parameters); +// auto otherMatrix = other.to_matrix(dimensions, parameters); +// auto returnMatrix = thisMatrix * otherMatrix; +// return returnMatrix; +// }; +// returnOperator.m_generator = CallbackFunction(func); +// } + +// else if (thisDegree != otherDegree) { +// std::cout << "defining function.\n"; +// auto func = [&](std::map dimensions, +// std::map> parameters) { +// std::cout << "inside function.\n"; +// auto thisMatrix = this->to_matrix(dimensions, parameters); +// std::cout << "calling _to_full_hilbert.\n"; +// // auto thisFullSpace = +// // _cast_to_full_hilbert(thisMatrix, thisDegree, allDegrees, +// dimensions); +// // std::cout << "back from _to_full_hilbert.\n"; +// // thisFullSpace.dump(); + +// // auto otherMatrix = other.to_matrix(dimensions, parameters); +// // auto otherFullSpace = +// // _cast_to_full_hilbert(otherMatrix, otherDegree, allDegrees, +// dimensions); +// // otherFullSpace.dump(); + +// // auto returnMatrix = thisFullSpace * otherFullSpace; +// // return returnMatrix; +// return thisFullSpace; +// }; +// std::cout << "function defined.\n"; +// returnOperator.m_generator = CallbackFunction(func); +// } + +// // else { +// // throw std::runtime_error("unexpected error encountered."); +// // } + +// return returnOperator; +// } + +// operator_sum +// elementary_operator::operator+(const elementary_operator &other) const { +// std::map merged_params = +// aggregate_parameters(this->m_ops, other.m_ops); + +// auto merged_func = +// [this, other](std::vector degrees, +// std::vector parameters) -> complex_matrix { +// complex_matrix result1 = this->func(degrees, parameters); +// complex_matrix result2 = other.func(degrees, parameters); + +// for (size_t i = 0; i < result1.rows(); i++) { +// for (size_t j = 0; j < result1.cols(); j++) { +// result1(i, j) += result2(i, j); +// } +// } + +// return result1; +// }; + +// return operator_sum(merged_func); +// } + +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/matrix_arithmetics.cpp b/runtime/cudaq/dynamics/matrix_arithmetics.cpp new file mode 100644 index 0000000000..c0c3f84ac3 --- /dev/null +++ b/runtime/cudaq/dynamics/matrix_arithmetics.cpp @@ -0,0 +1,165 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include +#include +#include +#include + +namespace cudaq { +Evaluated::Evaluated(const std::vector °rees, + const cudaq::complex_matrix &matrix) + : _degrees(degrees), _matrix(matrix) {} + +const std::vector &Evaluated::get_degrees() const { return _degrees; } + +const cudaq::complex_matrix &Evaluated::get_matrix() const { return _matrix; } + +matrix_arithmetics::matrix_arithmetics(const std::map &dimensions) + : dimensions(dimensions) {} + +std::vector +generate_all_states(const std::vector °rees, + const std::map &dimensions) { + // Determine total combinations based on degrees dimensions + int total_states = 1; + for (int deg : degrees) { + total_states *= dimensions.at(deg); + } + + std::vector states(total_states); + int num_states = total_states; + + for (size_t i = 0; i < degrees.size(); i++) { + int dim = dimensions.at(degrees[i]); + int repeats = num_states / dim; + int index = 0; + + for (int j = 0; j < total_states; j++) { + states[j] += std::to_string(index); + if ((j + 1) % repeats == 0) { + index = (index + 1) % dim; + } + } + } + + return states; +} + +std::vector +matrix_arithmetics::compute_permutation(const std::vector &op_degrees, + const std::vector &canon_degrees) { + // Generate all states for canonical degrees + std::vector states = + generate_all_states(canon_degrees, dimensions); + + // Determine the reordering indices + std::vector reordering; + for (int deg : op_degrees) { + auto it = std::find(canon_degrees.begin(), canon_degrees.end(), deg); + if (it != canon_degrees.end()) { + reordering.push_back(std::distance(canon_degrees.begin(), it)); + } + } + + // Generate the permutation vector + std::vector permutation; + for (const auto &state : states) { + std::string reordered_state; + for (int index : reordering) { + reordered_state += state[index]; + } + + auto perm_index = std::find(states.begin(), states.end(), reordered_state); + permutation.push_back(std::distance(states.begin(), perm_index)); + } + return permutation; +} + +std::pair> +matrix_arithmetics::canonicalize(const cudaq::complex_matrix &op_matrix, + const std::vector &op_degrees) { + std::vector canon_degrees = op_degrees; + std::sort(canon_degrees.begin(), canon_degrees.end()); + + if (op_degrees == canon_degrees) { + return std::make_pair(op_matrix, canon_degrees); + } + + // Compute permutation + std::vector permutation = compute_permutation(op_degrees, canon_degrees); + + cudaq::complex_matrix permuted_matrix = op_matrix; + for (int i = 0; i < permutation.size(); i++) { + auto row = op_matrix.get_row(permutation[i]); + permuted_matrix.set_row(i, row); + } + + return std::make_pair(permuted_matrix, canon_degrees); +} + +Evaluated matrix_arithmetics::evaluate(const Operator &op) const { + // cudaq::complex_matrix matrix = op.to_matrix(dimensions, op.parameters); + // // Replace with degrees from operator + // std::vector degrees = op.degrees; + // return Evaluated(degrees, matrix); + + return std::visit( + [this](auto &&arg) -> Evaluated { + using T = std::decay_t; + + if constexpr (std::is_same_v) { + cudaq::complex_matrix matrix = arg.to_matrix(dimensions, {}); + std::vector degrees = arg.degrees; + return Evaluated(degrees, matrix); + } else if constexpr (std::is_same_v) { + // Handle the scalar operator + if (arg.has_val()) { + // double scalar_value = arg.get_val(); + int dim = dimensions.at(arg.degrees[0]); + cudaq::complex_matrix matrix = + cudaq::complex_matrix::identity(dim, dim); + return Evaluated(arg.degrees, matrix); + } else { + throw std::runtime_error( + "Dynamic scalar operator evaluation not implemented."); + } + } else { + throw std::runtime_error( + "Unknown operator type in CudaqOperatorVariant."); + } + }, + op); +} + +// std::shared_ptr matrix_arithmetics::add(const +// std::shared_ptr &op1, const std::shared_ptr &op2) { +// assert(op1.degrees() == op2.degrees()); +// Matrix result = op1.matrix() + op2.matrix(); +// return std::make_shared(op1.degrees(), result); +// } + +// matrix_arithmetics::Evaluated matrix_arithmetics::mul(const Evaluated &op1, +// const Evaluated &op2) { +// assert(op1.degrees() == op2.degrees()); +// Matrix result = op1.matrix() * op2.matrix(); +// return Evaluated(op1.degrees(), result); +// } + +// matrix_arithmetics::Evaluated matrix_arithmetics::tensor(const Evaluated +// &op1, const Evaluated &op2) { +// std::vector new_degrees(op1.degrees()); +// new_degrees.insert(new_degrees.end(), op2.degrees().begin(), +// op2.degrees().end()); + +// Matrix result = Eigen::kroneckerProduct(op1.matrix(), +// op2.matrix()).eval(); auto [canonical_matrix, canonical_degrees] = +// canonicalize(result, new_degrees); return Evaluated(new_degrees, result); +// } + +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/matrix_arithmetics.h b/runtime/cudaq/dynamics/matrix_arithmetics.h new file mode 100644 index 0000000000..5e0260f3db --- /dev/null +++ b/runtime/cudaq/dynamics/matrix_arithmetics.h @@ -0,0 +1,50 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cudaq.h" +#include +#include +#include +#include + +namespace cudaq { +class Evaluated { +public: + Evaluated(const std::vector °rees, + const cudaq::complex_matrix &matrix); + const std::vector &get_degrees() const; + const cudaq::complex_matrix &get_matrix() const; + +private: + std::vector _degrees; + cudaq::complex_matrix _matrix; +}; + +class matrix_arithmetics : public operator_arithmetics { +public: + matrix_arithmetics(const std::map &dimensions); + + Evaluated evaluate(const Operator &op) const override; + // std::shared_ptr add(const std::shared_ptr &op1, const + // std::shared_ptr &op2) override; std::shared_ptr + // mul(const Evaluated &op1, const Evaluated &op2) override; + // std::shared_ptr tensor(const Evaluated &op1, const Evaluated + // &op2) override; + + std::vector compute_permutation(const std::vector &op_degrees, + const std::vector &canon_degrees); + std::pair> + canonicalize(const cudaq::complex_matrix &op_matrix, + const std::vector &op_degrees); + +private: + std::map dimensions; +}; +} // namespace cudaq diff --git a/runtime/cudaq/dynamics/operator_arithmetics.h b/runtime/cudaq/dynamics/operator_arithmetics.h new file mode 100644 index 0000000000..1470daecfe --- /dev/null +++ b/runtime/cudaq/dynamics/operator_arithmetics.h @@ -0,0 +1,34 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace cudaq { +class elementary_operator; +class scalar_operator; +class product_operator; +class operator_sum; + +template +class operator_arithmetics { +public: + virtual ~operator_arithmetics() = default; + + virtual TEval evaluate(const Operator &op) const = 0; + // virtual TEval add(const TEval &val1, const TEval &val2) = 0; + // virtual TEval mul(const TEval &val1, const TEval &val2) = 0; + // virtual TEval tensor(const TEval &val1, const TEval &val2) = 0; +}; +} // namespace cudaq diff --git a/runtime/cudaq/dynamics/operator_helpers.cpp b/runtime/cudaq/dynamics/operator_helpers.cpp new file mode 100644 index 0000000000..00dc28a305 --- /dev/null +++ b/runtime/cudaq/dynamics/operator_helpers.cpp @@ -0,0 +1,118 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "operator_helpers.h" +#include + +namespace cudaq { +std::map operator_helpers::aggregateParameters( + const std::vector> ¶meterMappings) { + std::map paramDescriptions; + for (const auto &descriptions : parameterMappings) { + for (const auto &[key, newDesc] : descriptions) { + std::string &existingDesc = paramDescriptions[key]; + if (!existingDesc.empty() && !newDesc.empty()) { + existingDesc += "\n---\n" + newDesc; + } else if (!newDesc.empty()) { + existingDesc = newDesc; + } + } + } + return paramDescriptions; +} + +std::string operator_helpers::parameterDocs(const std::string ¶mName, + const std::string &docs) { + if (paramName.empty() || docs.empty()) + return ""; + + try { + std::regex keywordPattern(R"(^\\s*(Arguments|Args):\\s*\\n)", + std::regex::multiline); + std::smatch match; + if (std::regex_search(docs, match, keywordPattern)) { + std::string paramDocs = match.suffix(); + std::regex paramPattern( + R"(^\s*" + paramName + R"(\s*(\(.*\))?:)\s*(.*)$)", + std::regex::multiline); + if (std::regex_search(paramDocs, match, paramPattern)) { + return match.str(2); + } + } + } catch (const std::exception &) { + return ""; + } + return ""; +} + +std::pair, std::map> +operator_helpers::argsFromKwargs(std::function fct, + std::map &kwargs) { + // FIXME implement later + return {}; +} + +std::vector +operator_helpers::generateAllStates(const std::vector °rees, + const std::map &dimensions) { + if (degrees.empty()) + return std::vector(); + + std::vector> states; + for (int state = 0; state < dimensions.at(degrees[0]); state++) { + states.push_back({std::to_string(state)}); + } + + for (size_t i = 1; i < degrees.size(); i++) { + std::vector> newStates; + int d = degrees[i]; + for (const auto ¤t : states) { + for (int state = 0; state < dimensions.at(d); state++) { + std::vector newState = current; + newState.push_back(std::to_string(state)); + newStates.push_back(newState); + } + } + states = newStates; + } + + std::vector result; + for (const auto &state : states) { + std::string combined; + for (const auto &s : state) { + combined += s; + } + result.push_back(combined); + } + return result; +} + +// void operator_helpers::permuteMatrix(Eigen::MatrixXcd &matrix, const +// std::vector &permutation) { +// Eigen::MatrixXcd permuted = matrix(permutation, +// Eigen::AllAtOnceTraversal); matrix = permuted; +// } + +Eigen::MatrixXcd +operator_helpers::cmatrixToNparray(const cudaq::complex_matrix &cmatrix) { + Eigen::MatrixXcd matrix(cmatrix.rows(), cmatrix.cols()); + for (int row = 0; row < cmatrix.rows(); row++) { + for (int col = 0; col < cmatrix.cols(); col++) { + matrix(row, col) = cmatrix(row, col); + } + } + return matrix; +} + +std::vector +operator_helpers::canonicalizeDegrees(const std::vector °rees) { + std::vector sortedDegrees = degrees; + std::sort(sortedDegrees.rbegin(), sortedDegrees.rend()); + return sortedDegrees; +} +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/operator_helpers.h b/runtime/cudaq/dynamics/operator_helpers.h new file mode 100644 index 0000000000..cdff4d4fed --- /dev/null +++ b/runtime/cudaq/dynamics/operator_helpers.h @@ -0,0 +1,54 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cudaq.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cudaq { +class operator_helpers { +public: + // Aggregate parameters from a list of mappings + static std::map aggregateParameters( + const std::vector> ¶meterMappings); + + // Extract Documentation for a specific parameter + static std::string parameterDocs(const std::string ¶mName, + const std::string &docs); + + // Extract arguments from keyword arguments + static std::pair, std::map> + argsFromKwargs(std::function fct, + std::map &kwargs); + + // Generate all possible states for the given degrees and dimensions + static std::vector + generateAllStates(const std::vector °rees, + const std::map &dimensions); + + // Permutes a matrix according to the given permutation + static void permuteMatrix(Eigen::MatrixXcd &matrix, + const std::vector &permutation); + + // Converts a cudaq::ComplexMatrix to Eigen::MatrixXcd + static Eigen::MatrixXcd cmatrixToNparray(const cudaq::complex_matrix &matrix); + + // Returns the degrees in canonical order (sorted in reverse) + static std::vector canonicalizeDegrees(const std::vector °ress); +}; +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/operator_sum.cpp b/runtime/cudaq/dynamics/operator_sum.cpp new file mode 100644 index 0000000000..a6991fda75 --- /dev/null +++ b/runtime/cudaq/dynamics/operator_sum.cpp @@ -0,0 +1,284 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "common/EigenDense.h" +#include "operator_arithmetics.h" +#include "operators.h" + +#include +#include + +namespace cudaq { + +/// Operator sum constructor given a vector of product operators. +operator_sum::operator_sum(const std::vector &terms) + : m_terms(terms) {} + +/// Canonicalize the terms of the operator sum. +std::vector> +operator_sum::_canonical_terms() const { + std::vector> + canonicalized_terms; + + for (const auto &term : m_terms) { + std::vector> term_ops = + term.get_terms(); + std::vector scalars; + std::vector non_scalars; + + for (const auto &op_variant : term_ops) { + std::visit( + [&scalars, &non_scalars](auto &op) { + if constexpr (std::is_same_v) { + scalars.push_back(op); + } else if constexpr (std::is_same_v) { + non_scalars.push_back(op); + } + }, + op_variant); + } + + // Sort scalars and non-scalars for canonical order + std::sort(non_scalars.begin(), non_scalars.end(), + [](const elementary_operator &a, const elementary_operator &b) { + return a.degrees < b.degrees; + }); + + if (scalars.empty()) { + scalars.push_back(scalar_operator(1.0)); + } + + for (const auto &scalar : scalars) { + for (const auto &non_scalar : non_scalars) { + canonicalized_terms.emplace_back(scalar, non_scalar); + } + } + } + + return canonicalized_terms; +} + +/// Canonicalize the operator sum into a new instance +operator_sum operator_sum::canonicalize() const { + std::vector canonical_terms; + for (const auto &[scalar, elementary] : _canonical_terms()) { + canonical_terms.emplace_back(product_operator({scalar, elementary})); + } + + return operator_sum(canonical_terms); +} + +/// Equality operator to compare canonicalized terms. +bool operator_sum::operator==(const operator_sum &other) const { + return _canonical_terms() == other._canonical_terms; +} + +/// Get unique degrees of freedom in canonical order. +std::vector operator_sum::degrees() const { + std::set unique_degrees; + for (const auto &term : m_terms) { + for (const auto &op_variant : term.get_terms()) { + std::visit( + [&unique_degrees](auto &op) { + unique_degrees.insert(op.degrees.begin(), op.degrees.end()); + }, + op_variant); + } + } + + return std::vector(unique_degrees.begin(), unique_degrees.end()); +} + +/// Aggregate parameters of all terms. +std::map> operator_sum::parameters() const { + std::map> param_map; + + for (const auto &term : m_terms) { + for (const auto &op_variant : term.get_terms()) { + std::visit( + [¶m_map](auto &op) { + const auto &op_params = op.parameters; + param_map.insert(op_params.begin(), op_params.end()); + }, + op_variant); + } + } + + return param_map; +} + +/// Check if the operator sum acts as a spin operator. +bool operator_sum::_is_spinop() const { + return std::all_of( + m_terms.begin(), m_terms.end(), [](const product_operator &term) { + return std::all_of(term.get_terms().begin(), term.get_terms().end(), + [](const auto &op_variant) { + return std::visit( + [](auto &op) { return op._is_spinop(); }, + op_variant); + }); + }); +} + +/// Evaluate the opetrator sum using provided arithmetic rules. +template +TEval operator_sum::_evaluate(operator_arithmetics &arithmetics) const { + // Collect all degrees of freedom from the terms + std::set degrees; + for (const auto &term : m_terms) { + for (auto &op_variant : term.get_terms()) { + std::visit( + [°rees](auto &op) { + degrees.insert(op.degrees.begin(), op.degrees.end()); + }, + op_variant); + } + } + + // Function to pad a product operator with identities for missing degrees + auto padded_term = [&](product_operator term) -> product_operator { + std::set term_degrees; + for (const auto &op_variant : term.get_terms()) { + std::visit( + [&term_degrees](auto &op) { + term.degrees.insert(op.degrees.begin(), op.degrees.end()); + }, + op_variant); + } + + // Add identity operators for missing degrees + for (int degree : degrees) { + if (term_degrees.find(degree) == term_degrees.end()) { + term *= elementary_operator::identity(degree); + } + } + return term; + }; + + // Evaluate the first term and initialize the sum + TEval sum = padded_term(m_terms[0])._evaluate(arithmetics); + + // Evaluate remaining terms and add them to the sum + for (size_t i = 1; i < m_terms.size(); i++) { + sum = arithmetics.add(sum, padded_term(m_terms[i])._evaluate(arithmetics)); + } + + return sum; +} + +// Arithmetic operators +operator_sum operator_sum::operator+(const operator_sum &other) const { + std::vector combined_terms = m_terms; + combined_terms.insert(combined_terms.end(), other.m_terms.begin(), + other.m_terms.end()); + return operator_sum(combined_terms).canonicalize(); +} + +/// FIXME: +// operator_sum operator_sum::operator-(const operator_sum &other) const { +// return *this + (-1 * other); +// } + +/// FIXME: +// operator_sum operator_sum::operator-=(const operator_sum &other) { +// *this = *this - other; +// return *this; +// } + +operator_sum operator_sum::operator+=(const operator_sum &other) { + *this = *this + other; + return *this; +} + +/// FIXME: +// operator_sum operator_sum::operator*(const operator_sum &other) const { +// std::vector product_terms; +// for (const auto &self_term : m_terms) { +// for (const auto &other_term : other.m_terms) { +// product_terms.push_back(self_term * other_term); +// } +// } +// return operator_sum(product_terms); +// } + +/// FIXME: +// operator_sum operator_sum::operator/(const operator_sum &other) const { +// std::vector divided_terms; +// for (const auto &term : m_terms) { +// divided_terms.push_back(term / other); +// } + +// return operator_sum(divided_terms); +// } + +operator_sum operator_sum::operator+(const scalar_operator &other) const { + std::vector combined_terms = m_terms; + combined_terms.push_back(product_operator({other})); + return operator_sum(combined_terms); +} + +operator_sum operator_sum::operator-(const scalar_operator &other) const { + return *this + (-1.0 * other); +} + +operator_sum operator_sum::operator+=(const scalar_operator &other) { + *this = *this + other; + return *this; +} + +operator_sum operator_sum::operator-=(const scalar_operator &other) { + *this = *this - other; + return *this; +} + +operator_sum operator_sum::operator+(const product_operator &other) const { + std::vector combined_terms = m_terms; + combined_terms.push_back(other); + return operator_sum(combined_terms); +} + +operator_sum operator_sum::operator+=(const product_operator &other) { + *this = *this + other; + return *this; +} + +/// FIXME: +// operator_sum operator_sum::operator-(const product_operator &other) const { +// return *this + (-1. * other); +// } + +/// FIXME: +// operator_sum operator_sum::operator-=(const product_operator &other) { +// *this = *this - other; +// return *this; +// } + +complex_matrix +operator_sum::to_matrix(const std::map &dimensions, + const std::map ¶ms) const { + complex_matrix result = m_terms[0].to_matrix(dimensions, params); + for (size_t i = 1; i < m_terms.size(); i++) { + result += m_terms[i].to_matrix(dimensions, params); + } + return result; +} + +/// Convert the operator sum to a string representation. +std::string operator_sum::to_string() const { + std::string result; + for (const auto &term : m_terms) { + result += term.to_string() + " + "; + } + // Remove last " + " + if (!result.empty()) + result.erase(result.size() - 3); + return result; +} + +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/operators.h b/runtime/cudaq/dynamics/operators.h new file mode 100644 index 0000000000..e1994bdbf4 --- /dev/null +++ b/runtime/cudaq/dynamics/operators.h @@ -0,0 +1,495 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cudaq/matrix.h" +#include "cudaq/qis/state.h" +#include "definition.h" + +#include +#include +#include +#include + +namespace cudaq { + +template +class operator_arithmetics; + +class operator_sum; +class product_operator; +class scalar_operator; +class elementary_operator; + +/// @brief Represents an operator expression consisting of a sum of terms, where +/// each term is a product of elementary and scalar operators. Operator +/// expressions cannot be used within quantum kernels, but they provide methods +/// to convert them to data types that can. +class operator_sum { +private: + std::vector m_terms; + + std::vector> + canonicalize_product(product_operator &prod) const; + + std::vector> + _canonical_terms() const; + +public: + /// @brief Empty constructor that a user can aggregate terms into. + operator_sum() = default; + + /// @brief Construct a `cudaq::operator_sum` given a sequence of + /// `cudaq::product_operator`'s. + /// This operator expression represents a sum of terms, where each term + /// is a product of elementary and scalar operators. + operator_sum(const std::vector &terms); + + operator_sum canonicalize() const; + + /// @brief The degrees of freedom that the operator acts on in canonical + /// order. + std::vector degrees() const; + + std::map> parameters() const; + + bool _is_spinop() const; + + /// TODO: implement + template + TEval _evaluate(operator_arithmetics &arithmetics) const; + + /// @brief Return the `operator_sum` as a matrix. + /// @arg `dimensions` : A mapping that specifies the number of levels, + /// that is, the dimension of each degree of freedom + /// that the operator acts on. Example for two, 2-level + /// degrees of freedom: `{0:2, 1:2}`. + /// @arg `parameters` : A map of the parameter names to their concrete, + /// complex values. + complex_matrix + to_matrix(const std::map &dimensions, + const std::map ¶ms = {}) const; + + // Arithmetic operators + operator_sum operator+(const operator_sum &other) const; + operator_sum operator-(const operator_sum &other) const; + operator_sum operator*(const operator_sum &other) const; + operator_sum operator/(const operator_sum &other) const; + operator_sum operator+=(const operator_sum &other); + operator_sum operator-=(const operator_sum &other); + + operator_sum operator+(const scalar_operator &other) const; + operator_sum operator-(const scalar_operator &other) const; + operator_sum operator+=(const scalar_operator &other); + operator_sum operator-=(const scalar_operator &other); + + operator_sum operator+(const product_operator &other) const; + operator_sum operator-(const product_operator &other) const; + operator_sum operator+=(const product_operator &other); + operator_sum operator-=(const product_operator &other); + + operator_sum operator+(const elementary_operator &other) const; + operator_sum operator-(const elementary_operator &other) const; + operator_sum operator+=(const elementary_operator &other); + operator_sum operator-=(const elementary_operator &other); + + /// @brief Return the operator_sum as a string. + std::string to_string() const; + + /// @brief True, if the other value is an operator_sum with equivalent terms, + /// and False otherwise. The equality takes into account that operator + /// addition is commutative, as is the product of two operators if they + /// act on different degrees of freedom. + /// The equality comparison does *not* take commutation relations into + /// account, and does not try to reorder terms block-wise; it may hence + /// evaluate to False, even if two operators in reality are the same. + /// If the equality evaluates to True, on the other hand, the operators + /// are guaranteed to represent the same transformation for all arguments. + bool operator==(const operator_sum &other) const; + + const std::vector &get_terms() const { return m_terms; } + + std::complex + evaluate(std::map> parameters); +}; + +/// @brief Represents an operator expression consisting of a product of +/// elementary and scalar operators. Operator expressions cannot be used within +/// quantum kernels, but they provide methods to convert them to data types +/// that can. +class product_operator : public operator_sum { +private: + /// FIXME: Not sure if this is quite correct but using it as a dummy + /// type for now. + std::vector> m_terms; + +public: + product_operator() = default; + ~product_operator() = default; + + /// @brief Constructor for an operator expression that represents a product + /// of scalar and elementary operators. + /// @arg atomic_operators : The operators of which to compute the product when + /// evaluating the operator expression. + product_operator( + std::vector> + atomic_operators); + + // Arithmetic overloads against all other operator types. + operator_sum operator+(std::complex other); + operator_sum operator-(std::complex other); + product_operator operator*(std::complex other); + product_operator operator*=(std::complex other); + product_operator operator/(const std::complex scalar) const; + operator_sum operator+(operator_sum other); + operator_sum operator-(operator_sum other); + product_operator operator*(operator_sum other); + product_operator operator*=(operator_sum other); + operator_sum operator+(const scalar_operator &other) const; + operator_sum operator-(const scalar_operator &other) const; + product_operator operator*(const scalar_operator &other) const; + product_operator &operator*=(const scalar_operator &other); + operator_sum operator+(const product_operator &other) const; + operator_sum operator-(const product_operator &other) const; + product_operator operator*(const product_operator &other) const; + product_operator operator*=(product_operator other); + operator_sum operator+(elementary_operator other); + operator_sum operator-(elementary_operator other); + product_operator operator*(elementary_operator other); + product_operator operator*=(elementary_operator other); + /// @brief True, if the other value is an operator_sum with equivalent terms, + /// and False otherwise. The equality takes into account that operator + /// addition is commutative, as is the product of two operators if they + /// act on different degrees of freedom. + /// The equality comparison does *not* take commutation relations into + /// account, and does not try to reorder terms block-wise; it may hence + /// evaluate to False, even if two operators in reality are the same. + /// If the equality evaluates to True, on the other hand, the operators + /// are guaranteed to represent the same transformation for all arguments. + bool operator==(product_operator other); + + /// @brief Converts a Pauli word to a product operator + product_operator _from_word(const std::string &word); + + /// @brief Return the `product_operator` as a string. + std::string to_string() const; + + /// @brief Return the `operator_sum` as a matrix. + /// @arg `dimensions` : A mapping that specifies the number of levels, + /// that is, the dimension of each degree of freedom + /// that the operator acts on. Example for two, 2-level + /// degrees of freedom: `{0:2, 1:2}`. + /// @arg `parameters` : A map of the parameter names to their concrete, + /// complex values. + complex_matrix + to_matrix(const std::map dimensions, + const std::map> parameters) const; + + /// @brief Creates a representation of the operator as a `cudaq::pauli_word` + /// that can be passed as an argument to quantum kernels. + // pauli_word to_pauli_word(); + + /// @brief The degrees of freedom that the operator acts on in canonical + /// order. + std::vector degrees() const; + + std::complex + evaluate(std::map> parameters); + + template + TEval _evaluate(operator_arithmetics &arithmetics) const; + + /// @brief A map of the parameter names to their concrete, complex values. + std::map> parameters; + + const std::vector> & + get_terms() const { + return m_terms; + }; +}; + +class elementary_operator : public product_operator { +private: + std::map m_ops; + +public: + // The constructor should never be called directly by the user: + // Keeping it internally documentd for now, however. + /// @brief Constructor. + /// @arg operator_id : The ID of the operator as specified when it was + /// defined. + /// @arg degrees : the degrees of freedom that the operator acts upon. + elementary_operator(std::string operator_id, std::vector degrees); + + // Copy constructor. + elementary_operator(const elementary_operator &other); + elementary_operator(elementary_operator &other); + + // Arithmetic overloads against all other operator types. + /// NOTE: All of below arithmetic implemented except for the division between + /// two elementary ops. Further testing at the matrix level is needed. + operator_sum operator+(std::complex other); + operator_sum operator-(std::complex other); + product_operator operator*(std::complex other); + product_operator operator/(std::complex other); + operator_sum operator+(double other); + operator_sum operator-(double other); + product_operator operator*(double other); + product_operator operator/(double other); + operator_sum operator+(scalar_operator other); + operator_sum operator-(scalar_operator other); + product_operator operator*(scalar_operator other); + product_operator operator/(scalar_operator other); + operator_sum operator+(elementary_operator other); + operator_sum operator-(elementary_operator other); + product_operator operator*(elementary_operator other); + product_operator operator/(elementary_operator other); + + /// TODO: implement and test the below + operator_sum operator+(operator_sum other); + operator_sum operator-(operator_sum other); + operator_sum operator+=(operator_sum other); + operator_sum operator-=(operator_sum other); + /// FIXME: Correct signature? + product_operator operator*(operator_sum other); + // + operator_sum operator+(product_operator other); + operator_sum operator-(product_operator other); + product_operator operator*(product_operator other); + product_operator operator/(product_operator other); + + /// @brief True, if the other value is an elementary operator with the same id + /// acting on the same degrees of freedom, and False otherwise. + bool operator==(elementary_operator other); + + /// @brief Return the `elementary_operator` as a string. + std::string to_string() const; + + /// @brief Return the `elementary_operator` as a matrix. + /// @arg `dimensions` : A map specifying the number of levels, + /// that is, the dimension of each degree of freedom + /// that the operator acts on. Example for two, 2-level + /// degrees of freedom: `{0 : 2, 1 : 2}`. + complex_matrix + to_matrix(const std::map dimensions, + const std::map> parameters) const; + + // Predefined operators. + static elementary_operator identity(int degree); + static elementary_operator zero(int degree); + static elementary_operator annihilate(int degree); + static elementary_operator create(int degree); + static elementary_operator momentum(int degree); + static elementary_operator number(int degree); + static elementary_operator parity(int degree); + static elementary_operator position(int degree); + /// FIXME: + static elementary_operator squeeze(int degree, + std::complex amplitude); + static elementary_operator displace(int degree, + std::complex amplitude); + + /// @brief Adds the definition of an elementary operator with the given id to + /// the class. After definition, an the defined elementary operator can be + /// instantiated by providing the operator id as well as the degree(s) of + /// freedom that it acts on. An elementary operator is a parameterized object + /// acting on certain degrees of freedom. To evaluate an operator, for example + /// to compute its matrix, the level, that is the dimension, for each degree + /// of freedom it acts on must be provided, as well as all additional + /// parameters. Additional parameters must be provided in the form of keyword + /// arguments. Note: The dimensions passed during operator evaluation are + /// automatically validated against the expected dimensions specified during + /// definition - the `create` function does not need to do this. + /// @arg operator_id : A string that uniquely identifies the defined operator. + /// @arg expected_dimensions : Defines the number of levels, that is the + /// dimension, + /// for each degree of freedom in canonical (that is sorted) order. A + /// negative or zero value for one (or more) of the expected dimensions + /// indicates that the operator is defined for any dimension of the + /// corresponding degree of freedom. + /// @arg create : Takes any number of complex-valued arguments and returns the + /// matrix representing the operator in canonical order. If the matrix + /// can be defined for any number of levels for one or more degree of + /// freedom, the `create` function must take an argument called + /// `dimension` (or `dim` for short), if the operator acts on a single + /// degree of freedom, and an argument called `dimensions` (or `dims` for + /// short), if the operator acts + /// on multiple degrees of freedom. + template + void define(std::string operator_id, std::map expected_dimensions, + Func create) { + if (m_ops.find(operator_id) != m_ops.end()) { + // todo: make a nice error message to say op already exists + throw; + } + auto defn = Definition(); + defn.create_definition(operator_id, expected_dimensions, create); + m_ops[operator_id] = defn; + } + + // Attributes. + + /// @brief The number of levels, that is the dimension, for each degree of + /// freedom in canonical order that the operator acts on. A value of zero or + /// less indicates that the operator is defined for any dimension of that + /// degree. + std::map expected_dimensions; + /// @brief The degrees of freedom that the operator acts on in canonical + /// order. + std::vector degrees; + /// @brief A map of the parameter names to their concrete, complex values. + /// This will be enabled once we can handle generalized callback function + /// arguments. + /// @FIXME: Not needed until generalizing the function arguments. + std::map> parameters; + std::string id; + + // /// @brief Creates a representation of the operator as `pauli_word` that + // can be passed as an argument to quantum kernels. + // pauli_word to_pauli_word ovveride(); + + std::complex + evaluate(std::map> parameters); + + template + TEval _evaluate(operator_arithmetics &arithmetics) const; +}; +// Reverse order arithmetic for elementary operators against pure scalars. +operator_sum operator+(std::complex other, elementary_operator self); +operator_sum operator-(std::complex other, elementary_operator self); +product_operator operator*(std::complex other, + elementary_operator self); +product_operator operator/(std::complex other, + elementary_operator self); +operator_sum operator+(double other, elementary_operator self); +operator_sum operator-(double other, elementary_operator self); +product_operator operator*(double other, elementary_operator self); +product_operator operator/(double other, elementary_operator self); + +class scalar_operator : public product_operator { +private: + // If someone gave us a constant value, we will just return that + // directly to them when they call `evaluate`. + std::complex m_constant_value; + +public: + /// @brief Constructor that just takes a callback function with no + /// arguments. + + scalar_operator(ScalarCallbackFunction &&create) { + generator = ScalarCallbackFunction(create); + } + + /// @brief Constructor that just takes and returns a complex double value. + /// @NOTE: This replicates the behavior of the python `scalar_operator::const` + /// without the need for an extra member function. + scalar_operator(std::complex value); + scalar_operator(double value); + + // Arithmetic overloads against other operator types. + scalar_operator operator+(scalar_operator other); + scalar_operator operator-(scalar_operator other); + scalar_operator operator*(scalar_operator other); + scalar_operator operator/(scalar_operator other); + /// TODO: + scalar_operator pow(scalar_operator other); + + /// TODO: All of the below need deeper testing but are implemented. + operator_sum operator+(elementary_operator other); + operator_sum operator-(elementary_operator other); + product_operator operator*(elementary_operator other); + product_operator operator/(elementary_operator other); + // TODO: + operator_sum operator+(operator_sum other); + operator_sum operator-(operator_sum other); + operator_sum operator*(operator_sum other); + operator_sum operator/(operator_sum other); + operator_sum operator+(product_operator other); + operator_sum operator-(product_operator other); + product_operator operator*(product_operator other); + product_operator operator/(product_operator other); + + /// @brief Return the scalar operator as a concrete complex value. + std::complex + evaluate(const std::map> parameters) const; + + // Return the scalar operator as a 1x1 matrix. This is needed for + // compatability with the other inherited classes. + complex_matrix + to_matrix(const std::map dimensions, + const std::map> parameters) const; + + // /// @brief Returns true if other is a scalar operator with the same + // /// generator. + // bool operator==(scalar_operator other); + + /// @brief The function that generates the value of the scalar operator. + /// The function can take a vector of complex-valued arguments + /// and returns a number. + ScalarCallbackFunction generator; + + // Only populated when we've performed arithmetic between various + // scalar operators. + std::vector _operators_to_compose; + + /// NOTE: We should revisit these constructors and remove any that have + /// become unnecessary as the implementation improves. + scalar_operator() = default; + // Copy constructor. + scalar_operator(const scalar_operator &other); + scalar_operator(scalar_operator &other); + + ~scalar_operator() = default; + + template + TEval _evaluate(operator_arithmetics &arithmetics) const; + + // Need this property for consistency with other inherited types. + // Particularly, to be used when the scalar operator is held within + // a variant type next to elementary operators. + std::vector degrees = {-1}; + + // REMOVEME: just using this as a temporary patch: + std::complex get_val() { return m_constant_value; }; + + bool has_val() { + return m_constant_value.real() != 0.0 && m_constant_value.imag() != 0.0; + }; +}; + +scalar_operator operator+(scalar_operator self, std::complex other); +scalar_operator operator-(scalar_operator self, std::complex other); +scalar_operator operator*(scalar_operator self, std::complex other); +scalar_operator operator/(scalar_operator self, std::complex other); +scalar_operator operator+(std::complex other, scalar_operator self); +scalar_operator operator-(std::complex other, scalar_operator self); +scalar_operator operator*(std::complex other, scalar_operator self); +scalar_operator operator/(std::complex other, scalar_operator self); +scalar_operator operator+(scalar_operator self, double other); +scalar_operator operator-(scalar_operator self, double other); +scalar_operator operator*(scalar_operator self, double other); +scalar_operator operator/(scalar_operator self, double other); +scalar_operator operator+(double other, scalar_operator self); +scalar_operator operator-(double other, scalar_operator self); +scalar_operator operator*(double other, scalar_operator self); +scalar_operator operator/(double other, scalar_operator self); +void operator+=(scalar_operator &self, std::complex other); +void operator-=(scalar_operator &self, std::complex other); +void operator*=(scalar_operator &self, std::complex other); +void operator/=(scalar_operator &self, std::complex other); +void operator+=(scalar_operator &self, scalar_operator other); +void operator-=(scalar_operator &self, scalar_operator other); +void operator*=(scalar_operator &self, scalar_operator other); +void operator/=(scalar_operator &self, scalar_operator other); + +using Operator = std::variant, + std::unique_ptr, + std::unique_ptr, + std::unique_ptr>; +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/product_operator.cpp b/runtime/cudaq/dynamics/product_operator.cpp new file mode 100644 index 0000000000..76891daea4 --- /dev/null +++ b/runtime/cudaq/dynamics/product_operator.cpp @@ -0,0 +1,202 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "common/EigenDense.h" +#include "operators.h" + +#include +#include +#include +#include + +namespace cudaq { + +/// Product Operator constructors. +product_operator::product_operator( + std::vector> + atomic_operators) + : m_terms(atomic_operators) {} + +complex_matrix kroneckerHelper(std::vector &matrices) { + // essentially we pass in the list of elementary operators to + // this function -- with lowest degree being leftmost -- then it computes the + // kronecker product of all of them. + auto kronecker = [](complex_matrix &self, complex_matrix &other) { + return self.kronecker(other); + }; + + return std::accumulate(begin(matrices), end(matrices), + complex_matrix::identity(1, 1), kronecker); +} + +/// Convert the product_operator to a matrix representation. +complex_matrix product_operator::to_matrix( + const std::map dimensions, + const std::map> parameters) const { + std::vector matricesFullVectorSpace; + + for (auto &term : m_terms) { + std::vector op_degrees; + complex_matrix operator_matrix; + + // Retrieve degrees and the matrix for the current operator + std::visit( + [&](auto &&op) { + op_degrees = op.degrees; + operator_matrix = op.to_matrix(dimensions, parameters); + }, + term); + + std::vector matrixWithIdentities; + + // Add identity matrices for missing degrees of freedom + for (const auto &[degree, level] : dimensions) { + if (std::find(op_degrees.begin(), op_degrees.end(), degree) == + op_degrees.end()) { + matrixWithIdentities.push_back(complex_matrix::identity(level, level)); + } else { + matrixWithIdentities.push_back(operator_matrix); + } + } + + // Compute the Kronecker product for this term + matricesFullVectorSpace.push_back(kroneckerHelper(matrixWithIdentities)); + } + + // Sum all Kronecker product to form the final matrix + if (matricesFullVectorSpace.empty()) { + throw std::runtime_error( + "No matrices to sum in product_operator::to_matrix."); + } + + complex_matrix result(dimensions.begin()->second, dimensions.end()->second); + result.set_zero(); + + for (const auto &matrix : matricesFullVectorSpace) { + result = result + matrix; + } + + return result; +} + +// Degrees property +std::vector product_operator::degrees() const { + std::set unique_degrees; + // The variant type makes it difficult + auto beginFunc = [](auto &&t) { return t.degrees.begin(); }; + auto endFunc = [](auto &&t) { return t.degrees.end(); }; + for (const auto &term : m_terms) { + unique_degrees.insert(std::visit(beginFunc, term), + std::visit(endFunc, term)); + } + // Erase any `-1` degree values that may have come from scalar operators. + auto it = unique_degrees.find(-1); + if (it != unique_degrees.end()) { + unique_degrees.erase(it); + } + return std::vector(unique_degrees.begin(), unique_degrees.end()); +} + +template +TEval elementary_operator::_evaluate( + operator_arithmetics &arithmetics) const { + std::cout << "In ProductOp _evaluate" << std::endl; + return arithmetics.evaluate(*this); +} + +operator_sum product_operator::operator+(const product_operator &other) const { + return operator_sum({*this, other}); +} + +operator_sum product_operator::operator-(const product_operator &other) const { + return *this + (-1.0 * other); +} + +product_operator +product_operator::operator*(const product_operator &other) const { + std::vector> + combined_terms = m_terms; + combined_terms.insert(combined_terms.end(), other.m_terms.begin(), + other.m_terms.end()); + return product_operator(combined_terms); +} + +product_operator +product_operator::operator/(const std::complex scalar) const { + return *this * (1.0 / scalar); +} + +product_operator product_operator::_from_word(const std::string &word) { + if (word.empty()) { + throw std::invalid_argument("Empty Pauli word!"); + } + + std::vector> ops; + + for (size_t i = 0; i < word.size(); i++) { + char c = std::tolower(word[i]); + switch (c) { + case 'x': + ops.push_back(elementary_operator("pauli_x", {static_cast(i)})); + break; + case 'y': + ops.push_back(elementary_operator("pauli_y", {static_cast(i)})); + break; + case 'z': + ops.push_back(elementary_operator("pauli_z", {static_cast(i)})); + break; + case 'i': + ops.push_back(elementary_operator::identity(static_cast(i))); + break; + default: + throw std::invalid_argument("Invalid character in Pauli word: " + c); + } + } + + return product_operator(ops); +} + +operator_sum product_operator::operator+(const scalar_operator &other) const { + return operator_sum({*this, product_operator({other})}); +} + +operator_sum product_operator::operator-(const scalar_operator &other) const { + return *this + (-1.0 * other); +} + +product_operator +product_operator::operator*(const scalar_operator &other) const { + std::vector> + combined_terms = m_terms; + combined_terms.insert(combined_terms.begin(), other); + return product_operator(combined_terms); +} + +product_operator &product_operator::operator*=(const scalar_operator &other) { + m_terms.insert(m_terms.begin(), other); + return *this; +} + +// product_operator product_operator::operator*=(std::complex other); +// operator_sum product_operator::operator+(operator_sum other); +// operator_sum product_operator::operator-(operator_sum other); +// product_operator product_operator::operator*(operator_sum other); +// product_operator product_operator::operator*=(operator_sum other); +// operator_sum product_operator::operator-(scalar_operator other); +// product_operator product_operator::operator*(scalar_operator other); +// product_operator product_operator::operator*=(scalar_operator other); +// operator_sum product_operator::operator+(product_operator other); +// operator_sum product_operator::operator-(product_operator other); +// product_operator product_operator::operator*(product_operator other); +// product_operator product_operator::operator*=(product_operator other); +// operator_sum product_operator::operator+(elementary_operator other); +// operator_sum product_operator::operator-(elementary_operator other); +// product_operator product_operator::operator*(elementary_operator other); +// product_operator product_operator::operator*=(elementary_operator other); + +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/dynamics/scalar_operator.cpp b/runtime/cudaq/dynamics/scalar_operator.cpp new file mode 100644 index 0000000000..db9866ea21 --- /dev/null +++ b/runtime/cudaq/dynamics/scalar_operator.cpp @@ -0,0 +1,253 @@ +/******************************************************************************* + * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "common/EigenDense.h" +#include "operators.h" + +#include +#include + +namespace cudaq { + +/// Constructors. +scalar_operator::scalar_operator(const scalar_operator &other) + : generator(other.generator), m_constant_value(other.m_constant_value) {} +scalar_operator::scalar_operator(scalar_operator &other) + : generator(other.generator), m_constant_value(other.m_constant_value) {} + +/// @brief Constructor that just takes and returns a complex double value. +scalar_operator::scalar_operator(std::complex value) { + m_constant_value = value; + auto func = [&](std::map> _none) { + return m_constant_value; + }; + generator = ScalarCallbackFunction(func); +} + +/// @brief Constructor that just takes a double and returns a complex double. +scalar_operator::scalar_operator(double value) { + std::complex castValue(value, 0.0); + m_constant_value = castValue; + auto func = [&](std::map> _none) { + return m_constant_value; + }; + generator = ScalarCallbackFunction(func); +} + +std::complex scalar_operator::evaluate( + const std::map> parameters) const { + return generator(parameters); +} + +complex_matrix scalar_operator::to_matrix( + const std::map dimensions, + const std::map> parameters) const { + complex_matrix returnOperator(1, 1); + returnOperator.set_zero(); + returnOperator(0, 0) = evaluate(parameters); + return returnOperator; +} + +template +TEval elementary_operator::_evaluate( + operator_arithmetics &arithmetics) const { + std::cout << "In ScalarOp _evaluate" << std::endl; + return arithmetics.evaluate(*this); +} + +#define ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES(op) \ + scalar_operator operator op(std::complex other, \ + scalar_operator self) { \ + /* Create an operator for the complex double value. */ \ + auto otherOperator = scalar_operator(other); \ + /* Create an operator that we will store the result in and return to the \ + * user. */ \ + scalar_operator returnOperator; \ + /* Store the previous generator functions in the new operator. This is \ + * needed as the old generator functions would effectively be lost once we \ + * leave this function scope. */ \ + returnOperator._operators_to_compose.push_back(self); \ + returnOperator._operators_to_compose.push_back(otherOperator); \ + auto newGenerator = \ + [&](std::map> parameters) { \ + /* FIXME: I have to use this hacky `.get_val()` on the newly created \ + * operator for the given complex double -- because calling the \ + * evaluate function returns 0.0 . I have no clue why??? */ \ + return returnOperator._operators_to_compose[0] \ + .evaluate(parameters) op returnOperator._operators_to_compose[1] \ + .get_val(); \ + }; \ + returnOperator.generator = ScalarCallbackFunction(newGenerator); \ + return returnOperator; \ + } + +#define ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_REVERSE(op) \ + scalar_operator operator op(scalar_operator self, \ + std::complex other) { \ + /* Create an operator for the complex double value. */ \ + auto otherOperator = scalar_operator(other); \ + /* Create an operator that we will store the result in and return to the \ + * user. */ \ + scalar_operator returnOperator; \ + /* Store the previous generator functions in the new operator. This is \ + * needed as the old generator functions would effectively be lost once we \ + * leave this function scope. */ \ + returnOperator._operators_to_compose.push_back(self); \ + returnOperator._operators_to_compose.push_back(otherOperator); \ + auto newGenerator = \ + [&](std::map> parameters) { \ + /* FIXME: I have to use this hacky `.get_val()` on the newly created \ + * operator for the given complex double -- because calling the \ + * evaluate function returns 0.0 . I have no clue why??? */ \ + return returnOperator._operators_to_compose[1] \ + .get_val() op returnOperator._operators_to_compose[0] \ + .evaluate(parameters); \ + }; \ + returnOperator.generator = ScalarCallbackFunction(newGenerator); \ + return returnOperator; \ + } + +#define ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_ASSIGNMENT(op) \ + void operator op(scalar_operator &self, std::complex other) { \ + /* Create an operator for the complex double value. */ \ + auto otherOperator = scalar_operator(other); \ + /* Need to move the existing generating function to a new operator so that \ + * we can modify the generator in `self` in-place. */ \ + scalar_operator copy(self); \ + /* Store the previous generator functions in the new operator. This is \ + * needed as the old generator functions would effectively be lost once we \ + * leave this function scope. */ \ + self._operators_to_compose.push_back(copy); \ + self._operators_to_compose.push_back(otherOperator); \ + auto newGenerator = \ + [&](std::map> parameters) { \ + /* FIXME: I have to use this hacky `.get_val()` on the newly created \ + * operator for the given complex double -- because calling the \ + * evaluate function returns 0.0 . I have no clue why??? */ \ + return self._operators_to_compose[0] \ + .evaluate(parameters) op self._operators_to_compose[1] \ + .get_val(); \ + }; \ + self.generator = ScalarCallbackFunction(newGenerator); \ + } + +#define ARITHMETIC_OPERATIONS_DOUBLES(op) \ + scalar_operator operator op(double other, scalar_operator self) { \ + std::complex value(other, 0.0); \ + return self op value; \ + } + +#define ARITHMETIC_OPERATIONS_DOUBLES_REVERSE(op) \ + scalar_operator operator op(scalar_operator self, double other) { \ + std::complex value(other, 0.0); \ + return value op self; \ + } + +#define ARITHMETIC_OPERATIONS_DOUBLES_ASSIGNMENT(op) \ + void operator op(scalar_operator &self, double other) { \ + std::complex value(other, 0.0); \ + self op value; \ + } + +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES(+); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES(-); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES(*); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES(/); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_REVERSE(+); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_REVERSE(-); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_REVERSE(*); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_REVERSE(/); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_ASSIGNMENT(+=); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_ASSIGNMENT(-=); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_ASSIGNMENT(*=); +ARITHMETIC_OPERATIONS_COMPLEX_DOUBLES_ASSIGNMENT(/=); +ARITHMETIC_OPERATIONS_DOUBLES(+); +ARITHMETIC_OPERATIONS_DOUBLES(-); +ARITHMETIC_OPERATIONS_DOUBLES(*); +ARITHMETIC_OPERATIONS_DOUBLES(/); +ARITHMETIC_OPERATIONS_DOUBLES_REVERSE(+); +ARITHMETIC_OPERATIONS_DOUBLES_REVERSE(-); +ARITHMETIC_OPERATIONS_DOUBLES_REVERSE(*); +ARITHMETIC_OPERATIONS_DOUBLES_REVERSE(/); +ARITHMETIC_OPERATIONS_DOUBLES_ASSIGNMENT(+=); +ARITHMETIC_OPERATIONS_DOUBLES_ASSIGNMENT(-=); +ARITHMETIC_OPERATIONS_DOUBLES_ASSIGNMENT(*=); +ARITHMETIC_OPERATIONS_DOUBLES_ASSIGNMENT(/=); + +#define ARITHMETIC_OPERATIONS_SCALAR_OPS(op) \ + scalar_operator scalar_operator::operator op(scalar_operator other) { \ + /* Create an operator that we will store the result in and return to the \ + * user. */ \ + scalar_operator returnOperator; \ + /* Store the previous generator functions in the new operator. This is \ + * needed as the old generator functions would effectively be lost once we \ + * leave this function scope. */ \ + returnOperator._operators_to_compose.push_back(*this); \ + returnOperator._operators_to_compose.push_back(other); \ + auto newGenerator = \ + [&](std::map> parameters) { \ + return returnOperator._operators_to_compose[0] \ + .evaluate(parameters) op returnOperator._operators_to_compose[1] \ + .evaluate(parameters); \ + }; \ + returnOperator.generator = ScalarCallbackFunction(newGenerator); \ + return returnOperator; \ + } + +/// FIXME: Broken implementation +#define ARITHMETIC_OPERATIONS_SCALAR_OPS_ASSIGNMENT(op) \ + void operator op(scalar_operator &self, scalar_operator other) { \ + /* Need to move the existing generating function to a new operator so \ + * that we can modify the generator in `self` in-place. */ \ + scalar_operator selfCopy(self); \ + /* Store the previous generator functions in the new operator. This is \ + * needed as the old generator functions would effectively be lost once we \ + * leave this function scope. */ \ + self._operators_to_compose.push_back(selfCopy); \ + self._operators_to_compose.push_back(other); \ + auto newGenerator = \ + [&](std::map> parameters) { \ + return self._operators_to_compose[0] \ + .evaluate(parameters) op self._operators_to_compose[1] \ + .evaluate(parameters); \ + }; \ + self.generator = ScalarCallbackFunction(newGenerator); \ + } + +ARITHMETIC_OPERATIONS_SCALAR_OPS(+); +ARITHMETIC_OPERATIONS_SCALAR_OPS(-); +ARITHMETIC_OPERATIONS_SCALAR_OPS(*); +ARITHMETIC_OPERATIONS_SCALAR_OPS(/); +ARITHMETIC_OPERATIONS_SCALAR_OPS_ASSIGNMENT(+=); +ARITHMETIC_OPERATIONS_SCALAR_OPS_ASSIGNMENT(-=); +ARITHMETIC_OPERATIONS_SCALAR_OPS_ASSIGNMENT(*=); +ARITHMETIC_OPERATIONS_SCALAR_OPS_ASSIGNMENT(/=); + +operator_sum scalar_operator::operator+(elementary_operator other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + return operator_sum({product_operator({*this}), product_operator({other})}); +} + +operator_sum scalar_operator::operator-(elementary_operator other) { + // Operator sum is composed of product operators, so we must convert + // both underlying types to `product_operators` to perform the arithmetic. + return operator_sum( + {product_operator({*this}), product_operator({-1. * other})}); +} + +product_operator scalar_operator::operator*(elementary_operator other) { + return product_operator({*this, other}); +} + +/// FIXME: division on elementary op needed +// product_operator scalar_operator::operator/(elementary_operator other) { +// return product_operator({*this, (1./other)}); +// } + +} // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/matrix.h b/runtime/cudaq/matrix.h index ff6c883e46..45610fbc88 100644 --- a/runtime/cudaq/matrix.h +++ b/runtime/cudaq/matrix.h @@ -48,6 +48,13 @@ class complex_matrix { complex_matrix(value_type *rawData, const std::size_t rows, const std::size_t cols); + complex_matrix(); + complex_matrix(const complex_matrix &other); + complex_matrix &operator=(const complex_matrix &other); + + static complex_matrix identity(const std::size_t rows, + const std::size_t cols); + /// @brief Return the internal data representation value_type *data() const; @@ -57,6 +64,11 @@ class complex_matrix { /// @brief Return the number of columns std::size_t cols() const { return nCols; } + /// @brief Return the total number of elements in the matrix + std::size_t size() const { return nRows * nCols; } + + complex_matrix operator+(const complex_matrix &other) const; + /// @brief Multiply this matrix with the provided other matrix. /// This does not modify this matrix but instead returns a /// new matrix value. @@ -69,6 +81,12 @@ class complex_matrix { /// @brief Return the element at the row `i` and column `j`. value_type &operator()(std::size_t i, std::size_t j) const; + /// @brief Check for equality between two complex matrices. + bool operator==(complex_matrix &other); + + /// @brief Return the exponential of the matrix. + complex_matrix exp() const; + /// @brief Return the minimal eigenvalue for this matrix. value_type minimal_eigenvalue() const; @@ -82,10 +100,24 @@ class complex_matrix { /// @brief Set all elements in this matrix to 0.0 void set_zero(); + /// @brief Set this matrix equal to the identity. + void set_identity(); + /// @brief Print this matrix to the given output stream void dump(std::ostream &os); /// @brief Print this matrix to standard out void dump(); + + /// @brief Kronecker product between two matrices. + complex_matrix kronecker(complex_matrix &other); + + std::vector get_row(std::size_t i) const; + + void set_row(std::size_t i, const std::vector &row); }; + +// Needed to perform an accumulation of kronecker products. +complex_matrix kronecker(complex_matrix &self, complex_matrix &other); + } // namespace cudaq \ No newline at end of file diff --git a/runtime/cudaq/spin/matrix.cpp b/runtime/cudaq/spin/matrix.cpp index 6a6642b80b..496e3683cb 100644 --- a/runtime/cudaq/spin/matrix.cpp +++ b/runtime/cudaq/spin/matrix.cpp @@ -12,6 +12,9 @@ #include #include +#include +#include + namespace cudaq { /// @brief Hash function for an Eigen::MatrixXcd @@ -40,6 +43,9 @@ std::unordered_map generalEigenSolvers; +complex_matrix::complex_matrix() + : internalOwnedData(nullptr), internalData(nullptr), nRows(0), nCols(0) {} + complex_matrix::complex_matrix(const std::size_t rows, const std::size_t cols) : internalOwnedData(std::unique_ptr( new complex_matrix::value_type[rows * cols])), @@ -53,6 +59,28 @@ complex_matrix::complex_matrix(complex_matrix::value_type *rawData, internalData = rawData; } +complex_matrix::complex_matrix(const complex_matrix &other) { + internalData = other.data(); + // std::memcpy(internalData, other.data(), sizeof(value_type) * other.size()); + nRows = other.rows(); + nCols = other.cols(); +} + +complex_matrix &complex_matrix::operator=(const complex_matrix &other) { + internalData = other.data(); + // std::memcpy(internalData, other.data(), sizeof(value_type) * other.size()); + nRows = other.rows(); + nCols = other.cols(); + return *this; +} + +complex_matrix complex_matrix::identity(const std::size_t rows, + const std::size_t cols) { + auto returnMatrix = complex_matrix(rows, cols); + returnMatrix.set_identity(); + return returnMatrix; +} + complex_matrix::value_type *complex_matrix::data() const { return internalData; } @@ -68,6 +96,25 @@ void complex_matrix::set_zero() { map.setZero(); } +void complex_matrix::set_identity() { + Eigen::Map map(internalData, nRows, nCols); + map.setIdentity(); +} + +complex_matrix complex_matrix::operator+(const complex_matrix &other) const { + if (nRows != other.nRows || nCols != other.nCols) { + throw std::runtime_error("Matrix dimensions must match for addition."); + } + + Eigen::Map map1(internalData, nRows, nCols); + Eigen::Map map2(other.data(), other.rows(), other.cols()); + Eigen::MatrixXcd result = map1 + map2; + + complex_matrix sum(result.rows(), result.cols()); + std::memcpy(sum.data(), result.data(), sizeof(value_type) * result.size()); + return sum; +} + complex_matrix complex_matrix::operator*(complex_matrix &other) const { Eigen::Map map(internalData, nRows, nCols); Eigen::Map otherMap(other.data(), other.nRows, other.nCols); @@ -99,6 +146,12 @@ complex_matrix::value_type &complex_matrix::operator()(std::size_t i, return map(i, j); } +bool complex_matrix::operator==(complex_matrix &other) { + Eigen::Map _self(internalData, nRows, nCols); + Eigen::Map _other(other.data(), nRows, nCols); + return (_self == _other); +} + std::vector complex_matrix::eigenvalues() const { Eigen::Map map(internalData, nRows, nCols); if (map.isApprox(map.adjoint())) { @@ -156,4 +209,56 @@ complex_matrix::value_type complex_matrix::minimal_eigenvalue() const { return eigenvalues()[0]; } -} // namespace cudaq +complex_matrix complex_matrix::exp() const { + Eigen::Map _self(internalData, nRows, nCols); + Eigen::MatrixXcd ret = _self.exp(); + complex_matrix copy(ret.rows(), ret.cols()); + std::memcpy(copy.data(), ret.data(), sizeof(value_type) * ret.size()); + return copy; +} + +complex_matrix complex_matrix::kronecker(complex_matrix &other) { + Eigen::Map map(internalData, nRows, nCols); + Eigen::Map otherMap(other.data(), other.rows(), + other.cols()); + Eigen::MatrixXcd ret = Eigen::kroneckerProduct(map, otherMap).eval(); + complex_matrix copy(ret.rows(), ret.cols()); + std::memcpy(copy.data(), ret.data(), + sizeof(std::complex) * ret.size()); + return copy; +} + +complex_matrix kronecker(complex_matrix &self, complex_matrix &other) { + Eigen::Map map(self.data(), self.rows(), self.cols()); + Eigen::Map otherMap(other.data(), other.rows(), + other.cols()); + Eigen::MatrixXcd ret = Eigen::kroneckerProduct(map, otherMap).eval(); + complex_matrix copy(ret.rows(), ret.cols()); + std::memcpy(copy.data(), ret.data(), + sizeof(std::complex) * ret.size()); + return copy; +} + +std::vector +complex_matrix::get_row(std::size_t i) const { + if (i >= nRows) + throw std::out_of_range("Row index is out of range"); + + std::vector row(nCols); + for (std::size_t j = 0; j < nCols; j++) { + row[j] = internalData[i * nCols + j]; + } + return row; +} + +void complex_matrix::set_row(std::size_t i, + const std::vector &row) { + if (i >= nRows || row.size() != nCols) + throw std::out_of_range("Row index is out of range or invalid row size"); + + for (std::size_t j = 0; j < nCols; j++) { + internalData[i * nCols + j] = row[j]; + } +} + +} // namespace cudaq \ No newline at end of file