Add support to apply unitary and non-unitary gates directly to StateV…

…ector from Python (#121) * Initial reformat and templated of gate tests * Redefine gate mappings with templates * Record given state * Move all gate classes to become StateVector methods * Fix compile errors for SV * Refactor gates implementation * Enable vector param methods * Enable label to gatename map * Support dispatch directly from StateVector class * Add apply methods to SV class * Add log2 utility function * Remove old file arch * Ensure bindings build support for float and double sizes complex data * Remove outdated tests for deleted modules * Remove old definition headers * Tidy dispatch map * Replace header with correct type * Ensure cpp17 is now used * Remove unneeded files in compilation * Ensure apply args are given in correct order * Allow make test to be run from working dir without cleaning * Fix cpp formatting * Remove old code and fix codefactor complaints * Add cast to enable wheel build on MacOS * Enable MSVC intrinsics if using Windows * Ensure C++17 as a requirement * Fix compile-time ifdef * Replace intrinsic with BSR * Avoid intrinsics for portability * Rename log2 function * Add support for 64bit and 128 complex numbers from C++ backend * Update bindings to allow class instantiation and method use * Enable support for different precision parameters passed to backend * Fix log2 change * Fix binding names * Remove io from statevector * Refix the log2 -- replace with instrinsics later * Update format and remove ununsed warnings * Refactor the testing infrastructure for templated StateVector implementation (#115) * Add preliminary catch2 support * Move private methods to public for testing * Overload applyOperations for param and non param calls * Add testing support for X,Y,Z,H gates * Add S, T gate tests * Add support for RX,RY,RZ gate tests * Add PhaseShift gate tests * Add Rot tests * Add CNOT tests * Add support for CSWAP, Toffoli, CZ, CRot tests * Update testing CI * Fix CodeFactor complaints * Update CI image before running tests * Favour use of reverse iterator over counter in for loops * Fix contructor tests * Fix formatting * Fix narrowing complaints * Fix MSVC math errors * Run black * Remove unneeded ops for real*complex products * Fix test builder * Ensure cmake has a version to avoid warnings * Fix formatting of SV * Fix test reporting * Remove whitespace for CF complaints * Apply static analyser fixes * Add imaginary utils * Refactor gate implementation and utility definitions * Fix RY gate defn * Add compile-time complex multiplication functions * Use gate definition functions in tests and add constexpr where applicable * Remove outdated test utilities and tests * Update changelog * Begin support additions fotr arbitrary unitary application * Enable subproject test builds * Enable standalone CPP library header * Update sizing checks for arbitrary unitary application * Add additional utility functions * Add tests for arbitrary unitary application * Add tests for utils and helper functions * Add testing and safety checks for utilities * Remove LAPACK from build system * Add placeholder functions for gate multiplication * Format TestHelpers * Format utility tests * Add unitary apply to bindings * Remove unneeded functions * Relabel matrix applicator methods * Appease the codefactor * Split tests into param and nonparam * Add support for numpy 2D array gate data without copies * Appease codefactor * Offload dim check to util function * Add docstring for utility functions * Codefactor fixes
PennyLaneAI · Aug 9, 2021 · 800eaf3 · 800eaf3
1 parent e309d44
commit 800eaf3
Show file tree

Hide file tree

Showing 9 changed files with 990 additions and 421 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -4,18 +4,18 @@ project(pennylane_lightning)
 
 # Read and set pennylane_lightning version
 function(set_pennylane_lightning_version VERSION_FILE_PATH)
-	file(STRINGS ${VERSION_FILE_PATH} VERSION_FILE_STR)
-	foreach (LINE IN LISTS VERSION_FILE_STR)
-	if("${LINE}" MATCHES "__version__.*")
-		set(VERSION_LINE_STR "${LINE}")
-	endif()
-	endforeach()
-
-	string(REGEX REPLACE "__version__ = \"(.*)\"" "\\1" VERSION_STRING ${VERSION_LINE_STR})
-	set(VERSION_STRING ${VERSION_STRING} PARENT_SCOPE)
+    file(STRINGS ${VERSION_FILE_PATH} VERSION_FILE_STR)
+    foreach (LINE IN LISTS VERSION_FILE_STR)
+    if("${LINE}" MATCHES "__version__.*")
+        set(VERSION_LINE_STR "${LINE}")
+    endif()
+    endforeach()
+
+    string(REGEX REPLACE "__version__ = \"(.*)\"" "\\1" VERSION_STRING ${VERSION_LINE_STR})
+    set(VERSION_STRING ${VERSION_STRING} PARENT_SCOPE)
 endfunction()
 
-set_pennylane_lightning_version("${CMAKE_SOURCE_DIR}/pennylane_lightning/_version.py")
+set_pennylane_lightning_version(${PROJECT_SOURCE_DIR}/pennylane_lightning/_version.py)
 
 message(STATUS "pennylane_lightning version ${VERSION_STRING}")
 set(PROJECT_VERSION ${VERSION_STRING})
@@ -24,14 +24,12 @@ set(CMAKE_CXX_STANDARD 17) # At least C++17 is required
 find_package(OpenMP REQUIRED) # find OpenMP
 
 if(NOT CMAKE_BUILD_TYPE)
-	set(CMAKE_BUILD_TYPE Release)
+    set(CMAKE_BUILD_TYPE Release)
 endif()
 
-
 option(ENABLE_NATIVE "Enable native CPU build tuning" OFF)
 option(BUILD_TESTS "Build cpp tests" OFF)
 
-
 # Add pybind11
 include(FetchContent)
 FetchContent_Declare(
@@ -46,12 +44,6 @@ target_include_directories(pennylane_lightning INTERFACE "pennylane_lightning/sr
 
 add_library(external_dependency INTERFACE)
 
-if ("$ENV{USE_LAPACK}" OR "${USE_LAPACK}")
-    message(STATUS "Use LAPACKE")
-    target_link_libraries(external_dependency INTERFACE lapacke)
-    target_compile_options(external_dependency INTERFACE "-DLAPACKE=1")
-endif()
-
 if ("$ENV{USE_OPENBLAS}" OR "${USE_OPENBLAS}")
     message(STATUS "Use OPENBLAS")
     target_link_libraries(external_dependency INTERFACE openblas)
@@ -74,6 +66,7 @@ if(ENABLE_NATIVE)
     target_compile_options(lightning_qubit_ops PRIVATE -march=native)
 endif()
 
+
 if (BUILD_TESTS)
     enable_testing()
     add_subdirectory("pennylane_lightning/src/tests" "tests")

diff --git a/pennylane_lightning/src/Bindings.cpp b/pennylane_lightning/src/Bindings.cpp
@@ -11,6 +11,7 @@
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 // See the License for the specific language governing permissions and
 // limitations under the License.
+
 #include "StateVector.hpp"
 #include "pybind11/complex.h"
 #include "pybind11/numpy.h"
@@ -59,17 +60,67 @@ template <class fp_t> class StateVecBinder : public StateVector<fp_t> {
                const vector<vector<fp_t>> &params) {
         this->applyOperations(ops, wires, inverse, params);
     }
+    void apply(const vector<string> &ops, const vector<vector<size_t>> &wires,
+               const vector<bool> &inverse) {
+        this->applyOperations(ops, wires, inverse);
+    }
+    /**
+     * @brief Directly apply a given matrix to the specified wires. Matrix data
+     * in 1D row-major format.
+     *
+     * @param matrix
+     * @param wires
+     */
+    void applyMatrixWires(const std::vector<std::complex<fp_t>> &matrix,
+                          const vector<size_t> &wires, bool inverse = false) {
+        this->applyOperation(matrix, wires, inverse);
+    }
+    /**
+     * @brief Directly apply a given matrix to the specified wires. Data in 1/2D
+     * numpy complex array format.
+     *
+     * @param matrix
+     * @param wires
+     * @param inverse
+     */
+    void applyMatrixWires(const py::array_t<complex<fp_t>> &matrix,
+                          const vector<size_t> &wires, bool inverse = false) {
+        const vector<size_t> internalIndices = this->generateBitPatterns(wires);
+        const vector<size_t> externalWires =
+            this->getIndicesAfterExclusion(wires);
+        const vector<size_t> externalIndices =
+            this->generateBitPatterns(externalWires);
+        this->applyMatrix(static_cast<complex<fp_t> *>(matrix.request().ptr),
+                          internalIndices, externalIndices, inverse);
+    }
 };
 
 template <class PrecisionT> void lightning_class_bindings(py::module &m) {
     // Enable module name to be based on size of complex datatype
-    const std::string bitsize = std::to_string(sizeof(PrecisionT) * 8 * 2);
+    const std::string bitsize =
+        std::to_string(sizeof(std::complex<PrecisionT>) * 8);
     const std::string class_name = "StateVectorC" + bitsize;
     py::class_<StateVecBinder<PrecisionT>>(m, class_name.c_str())
         .def(py::init<
              py::array_t<complex<PrecisionT>,
                          py::array::c_style | py::array::forcecast> &>())
-        .def("apply", &StateVecBinder<PrecisionT>::apply)
+        .def("apply",
+             py::overload_cast<
+                 const vector<string> &, const vector<vector<size_t>> &,
+                 const vector<bool> &, const vector<vector<PrecisionT>> &>(
+                 &StateVecBinder<PrecisionT>::apply))
+        .def("apply", py::overload_cast<const vector<string> &,
+                                        const vector<vector<size_t>> &,
+                                        const vector<bool> &>(
+                          &StateVecBinder<PrecisionT>::apply))
+        .def("applyMatrix",
+             py::overload_cast<const std::vector<std::complex<PrecisionT>> &,
+                               const vector<size_t> &, bool>(
+                 &StateVecBinder<PrecisionT>::applyMatrixWires))
+        .def("applyMatrix",
+             py::overload_cast<const py::array_t<complex<PrecisionT>> &,
+                               const vector<size_t> &, bool>(
+                 &StateVecBinder<PrecisionT>::applyMatrixWires))
         .def("PauliX", &StateVecBinder<PrecisionT>::applyPauliX)
         .def("PauliY", &StateVecBinder<PrecisionT>::applyPauliY)
         .def("PauliZ", &StateVecBinder<PrecisionT>::applyPauliZ)

diff --git a/pennylane_lightning/src/StateVector.hpp b/pennylane_lightning/src/StateVector.hpp
@@ -140,6 +140,34 @@ template <class fp_t = double> class StateVector {
         gate(internalIndices, externalIndices, inverse, params);
     }
 
+    /**
+     * @brief Apply a single gate to the state-vector.
+     *
+     * @param matrix Arbitrary unitary gate to apply.
+     * @param wires Wires to apply gate to.
+     * @param inverse Indicates whether to use inverse of gate.
+     * @param params Optional parameter list for parametric gates.
+     */
+    void applyOperation(const std::vector<CFP_t> &matrix,
+                        const vector<size_t> &wires, bool inverse = false,
+                        [[maybe_unused]] const vector<fp_t> &params = {}) {
+
+        auto dim = Util::dimSize(matrix);
+
+        if (dim != wires.size())
+            throw std::invalid_argument(string("The supplied gate requires ") +
+                                        std::to_string(dim) + " wires, but " +
+                                        std::to_string(wires.size()) +
+                                        " were supplied.");
+
+        const vector<size_t> internalIndices = generateBitPatterns(wires);
+        const vector<size_t> externalWires = getIndicesAfterExclusion(wires);
+        const vector<size_t> externalIndices =
+            generateBitPatterns(externalWires);
+
+        applyMatrix(matrix, internalIndices, externalIndices, inverse);
+    }
+
     /**
      * @brief Apply multiple gates to the state-vector.
      *
@@ -234,11 +262,67 @@ template <class fp_t = double> class StateVector {
         return generateBitPatterns(qubitIndices, num_qubits_);
     }
 
-    // Apply Gates
-    void applyUnitary(const vector<CFP_t> &matrix,
-                      const vector<size_t> &indices,
-                      const vector<size_t> &externalIndices, bool inverse) {
-        if (indices.size() != length_)
+    /**
+     * @brief Apply a given matrix directly to the statevector.
+     *
+     * @param matrix Perfect square matrix in row-major order.
+     * @param indices Internal indices participating in the operation.
+     * @param externalIndices External indices unaffected by the operation.
+     * @param inverse Indicate whether inverse should be taken.
+     */
+    void applyMatrix(const vector<CFP_t> &matrix, const vector<size_t> &indices,
+                     const vector<size_t> &externalIndices, bool inverse) {
+        if (static_cast<size_t>(0b1 << (Util::log2(indices.size()) +
+                                        Util::log2(externalIndices.size()))) !=
+            length_)
+            throw std::out_of_range(
+                "The given indices do not match the state-vector length.");
+
+        vector<CFP_t> v(indices.size());
+        for (const size_t &externalIndex : externalIndices) {
+            CFP_t *shiftedState = arr_ + externalIndex;
+            // Gather
+            size_t pos = 0;
+            for (const size_t &index : indices) {
+                v[pos] = shiftedState[index];
+                pos++;
+            }
+
+            // Apply + scatter
+            for (size_t i = 0; i < indices.size(); i++) {
+                size_t index = indices[i];
+                shiftedState[index] = 0;
+
+                if (inverse == true) {
+                    for (size_t j = 0; j < indices.size(); j++) {
+                        const size_t baseIndex = j * indices.size();
+                        shiftedState[index] +=
+                            conj(matrix[baseIndex + i]) * v[j];
+                    }
+                } else {
+                    const size_t baseIndex = i * indices.size();
+                    for (size_t j = 0; j < indices.size(); j++) {
+                        shiftedState[index] += matrix[baseIndex + j] * v[j];
+                    }
+                }
+            }
+        }
+    }
+
+    /**
+     * @brief Apply a given matrix directly to the statevector read directly
+     * from numpy data. Data can be in 1D or 2D format.
+     *
+     * @param matrix Pointer from numpy data.
+     * @param indices Internal indices participating in the operation.
+     * @param externalIndices External indices unaffected by the operation.
+     * @param inverse Indicate whether inverse should be taken.
+     */
+    void applyMatrix(const CFP_t *matrix, const vector<size_t> &indices,
+                     const vector<size_t> &externalIndices, bool inverse) {
+        if (static_cast<size_t>(0b1 << (Util::log2(indices.size()) +
+                                        Util::log2(externalIndices.size()))) !=
+            length_)
             throw std::out_of_range(
                 "The given indices do not match the state-vector length.");
 
@@ -325,7 +409,6 @@ template <class fp_t = double> class StateVector {
 
     void applyT(const vector<size_t> &indices,
                 const vector<size_t> &externalIndices, bool inverse) {
-
         const CFP_t shift = (inverse == true)
                                 ? std::conj(std::exp(CFP_t(0, M_PI / 4)))
                                 : std::exp(CFP_t(0, M_PI / 4));

diff --git a/pennylane_lightning/src/Util.hpp b/pennylane_lightning/src/Util.hpp
@@ -105,10 +105,10 @@ template <class T> inline static constexpr T INVSQRT2() {
 }
 
 /**
- * Calculates 2^n -1 for some integer n > 0 using bitshifts.
+ * Calculates 2^n for some integer n > 0 using bitshifts.
  *
  * @param n the exponent
- * @return value of 2^n -1
+ * @return value of 2^n
  */
 inline size_t exp2(const size_t &n) { return static_cast<size_t>(1) << n; }
 
@@ -134,15 +134,30 @@ inline size_t maxDecimalForQubit(size_t qubitIndex, size_t qubits) {
     return exp2(qubits - qubitIndex - 1);
 }
 
+/**
+ * @brief Returns the number of wires supported by a given qubit gate.
+ *
+ * @tparam T Floating point precision type.
+ * @param data Gate matrix data.
+ * @return size_t Number of wires.
+ */
+template <class T> inline size_t dimSize(const std::vector<T> &data) {
+    const size_t s = data.size();
+    const size_t s_sqrt = std::sqrt(s);
+
+    if (s < 4)
+        throw std::invalid_argument("The dataset must be at least 2x2.");
+    if (((s == 0) || (s & (s - 1))))
+        throw std::invalid_argument("The dataset must be a power of 2");
+    if (s_sqrt * s_sqrt != s)
+        throw std::invalid_argument("The dataset must be a perfect square");
+
+    return log2(sqrt(data.size()));
+}
+
 } // namespace Util
 } // namespace Pennylane
 
-// Helper similar to std::make_unique from c++14
-template <typename T, typename... Args>
-std::unique_ptr<T> make_unique(Args &&...args) {
-    return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
-}
-
 // Exception for functions that aren't implemented
 class NotImplementedException : public std::logic_error {
   public:

diff --git a/pennylane_lightning/src/tests/CMakeLists.txt b/pennylane_lightning/src/tests/CMakeLists.txt
@@ -1,4 +1,5 @@
 cmake_minimum_required(VERSION 3.14)
+
 project(pennylane_lightning_tests)
 
 set(CMAKE_CXX_STANDARD 17)
@@ -31,8 +32,11 @@ include(Catch)
 add_executable(runner runner_main.cpp)
 target_link_libraries(runner pennylane_lightning Catch2::Catch2)
 
-target_sources(runner PRIVATE Test_StateVector.cpp
-                              Test_Bindings.cpp)
+target_sources(runner PRIVATE   Test_Bindings.cpp
+                                Test_StateVector_Nonparam.cpp 
+                                Test_StateVector_Param.cpp 
+                                Test_Util.cpp
+)
 
 target_compile_options(runner PRIVATE "$<$<CONFIG:DEBUG>:-Wall>")
 

diff --git a/pennylane_lightning/src/tests/TestHelpers.hpp b/pennylane_lightning/src/tests/TestHelpers.hpp
@@ -0,0 +1,66 @@
+#include <complex>
+#include <vector>
+
+/**
+ * @brief Utility function to compare complex statevector data.
+ *
+ * @tparam Data_t Floating point data-type.
+ * @param data1 StateVector data 1.
+ * @param data2 StateVector data 2.
+ * @return true Data are approximately equal.
+ * @return false Data are not approximately equal.
+ */
+template <class Data_t>
+inline bool isApproxEqual(
+    const std::vector<Data_t> &data1, const std::vector<Data_t> &data2,
+    const typename Data_t::value_type eps =
+        std::numeric_limits<typename Data_t::value_type>::epsilon() * 100) {
+    if (data1.size() != data2.size())
+        return false;
+
+    for (size_t i = 0; i < data1.size(); i++) {
+        if (data1[i].real() != Approx(data2[i].real()).epsilon(eps) ||
+            data1[i].imag() != Approx(data2[i].imag()).epsilon(eps)) {
+            return false;
+        }
+    }
+    return true;
+}
+
+/**
+ * @brief Utility function to compare complex statevector data.
+ *
+ * @tparam Data_t Floating point data-type.
+ * @param data1 StateVector data 1.
+ * @param data2 StateVector data 2.
+ * @return true Data are approximately equal.
+ * @return false Data are not approximately equal.
+ */
+template <class Data_t>
+inline bool
+isApproxEqual(const Data_t &data1, const Data_t &data2,
+              const typename Data_t::value_type eps =
+                  std::numeric_limits<typename Data_t::value_type>::epsilon() *
+                  100) {
+    if (data1.real() != Approx(data2.real()).epsilon(eps) ||
+        data1.imag() != Approx(data2.imag()).epsilon(eps)) {
+        return false;
+    }
+    return true;
+}
+
+/**
+ * @brief Multiplies every value in a dataset by a given complex scalar value.
+ *
+ * @tparam Data_t Precision of complex data type. Supports float and double
+ * data.
+ * @param data Data to be scaled.
+ * @param scalar Scalar value.
+ */
+template <class Data_t>
+void scaleVector(std::vector<std::complex<Data_t>> &data,
+                 std::complex<Data_t> scalar) {
+    std::transform(
+        data.begin(), data.end(), data.begin(),
+        [scalar](const std::complex<Data_t> &c) { return c * scalar; });
+}