From dd67eacc0e24c233b525850ef9bda56e96018ddd Mon Sep 17 00:00:00 2001 From: Paul Druce <30435358+pauldruce@users.noreply.github.com> Date: Tue, 21 Nov 2023 23:50:55 +0000 Subject: [PATCH] Improving simulation class (#6) * Adding ability to retrieve a ref to dirac operator from simulation class. This change is the start of making this library actually useful for calculations. * Add an eigenvalue recorder for DiracOperator A new method was added to DiracOperator to compute its eigenvalues. Additionally, a Recorder class was created in the Type13Metropolis example to log these eigenvalues into the chosen output file. This enhancement will aid in debugging, performance monitoring, and further scientific analyses during the simulation process. * Add HDF5 support and save eigenvalues to file Added Armadillo's HDF5 support in both 'CMakeLists.txt' and 'Type13Metropolis.cpp'. 'Type13Metropolis.cpp' now includes a 'Recorder' class which can save simulation data to file in HDF5 format. This change aids data storage and analysis as HDF5 allows for easy recall of larger datasets. Also updated CMake files to find and link with HDF5 libraries. * Added getEigenvalues method and updated dependencies In the DiracOperator.cpp and DiracOperator.hpp files, the getEigenvalues method was extracted from the header file and moved to the cpp file for better code organization and encapsulation. This function returns the eigenvalues of the Dirac Operator. Also, the library 'libhdf5-serial-dev' was added to the list of dependencies in .github/workflows/compatability_tests.yml and .github/workflows/build_and_test.yml to ensure successful compilation and testing across different platforms. This library provides the necessary tools for working with hdf5 data files, which are used in this application for data storage and manipulation. The build and test stages were also renamed in the compatibility tests workflow for clarity and consistency. This renaming helps to avoid confusion when viewing the results in GitHub Actions. * Add EigenvalueRecorder and refactor HDF5-related options Added the EigenvalueRecorder class to record and save eigenvalues. This class takes a "DiracOperator" object, a coupling constant "g2", and a "simulationId" string to record the eigenvalues from the Dirac operator into a HDF5 file labelled by the simulation id. It was created to provide a more effective and organized way to handle and record eigenvalues for different simulations. Furthermore, the definition of ARMA_USE_HDF5 macro was moved from the global CMakeLists.txt file to the EigenvalueRecorder.hpp so that it's only defined where actually needed. This was done to avoid potential conflicts with other parts of the code where HDF5 is not required. Lastly, the EigenvalueRecorder files were added into the CMakeLists.txt file to ensure they are included in the build. * Refactor Type13Metropolis.cpp, update build dependencies, and correct a comment typo Refactored Type13Metropolis.cpp file to simplify the code. A large portion of hand-written HDF5 I/O code was replaced with the more streamlined EigenvalueRecorder class, which reduced the length of the file significantly. In the build_and_test.yml file, a long line of dependency installation commands was broken down into multiple lines for better readability and easier maintenance. A minor typo in Simulation.hpp header file comment was corrected for improved code documentation clarity. These changes were made to make the codebase more maintainable and the comments more accurate. * Add documentation to EigenvalueRecorder class Added detailed documentation to EigenvalueRecorder class to enhance code readability and understanding. The comments provide a clear understanding of the purpose and functionality of EigenvalueRecorder class and its recordEigenvalues method. These improvements aim to ease the process of code maintenance and modification. * Add eigenvalue test and update DiracOperator comments Added a unit test for DiracOperator to verify that the eigenvalues are being calculated correctly. * "Refactor Simulation testing and update DiracOperator and Hamiltonian classes The Simulation unit test has been updated to create a temporary variable `dirac` of a Dirac operator, and validate not only its type but also if the Dirac matrix is correct. The `Hamiltonian` class was improved by replacing the description, adding more comments to empty `TODO` areas, and replacing dynamic allocation of `double` arrays with `std::vectors`. This change provides safer memory management and prevents memory leaks. The `DiracOperator` class got a copy constructor. This constructor makes a deep copy of all aspects of the operator, ensuring accurate copying of state details. Memory allocation for storing previous configuration was improved by replacing a dynamic array `mat_bk` with a vector, leading to cleaner and memory safe code. Furthermore, the `Hamiltonian` class got updated with safer memory handling by using `std::vector` instead of raw pointers. Raw pointers were deleted later in the code causing potential memory leaks. The Clifford-type was moved to private data members of DiracOperator for encapsulation. Added the task in TODO file for making an interface for DiracOperator" * Tidying github actions files --- .github/workflows/build_and_test.yml | 92 ++++---- .github/workflows/combination_selection.py | 19 +- .github/workflows/compatability_tests.yml | 32 ++- .github/workflows/linter.yml | 2 - .gitignore | 1 + CMakeLists.txt | 4 +- README.md | 13 +- RFL_source/README.md | 29 ++- .../new_source/BarrettGlaser/Action.cpp | 22 +- .../new_source/BarrettGlaser/Action.hpp | 8 +- .../new_source/BarrettGlaser/Hamiltonian.cpp | 190 +++++++-------- .../new_source/BarrettGlaser/Hamiltonian.hpp | 83 ++++--- .../new_source/BarrettGlaser/Metropolis.cpp | 90 ++++---- .../new_source/BarrettGlaser/Metropolis.hpp | 38 +-- RFL_source/new_source/CMakeLists.txt | 5 +- RFL_source/new_source/Clifford.cpp | 53 ++--- RFL_source/new_source/Clifford.hpp | 2 +- RFL_source/new_source/DiracOperator.cpp | 151 +++++++----- RFL_source/new_source/DiracOperator.hpp | 80 ++++--- RFL_source/new_source/GslRng.hpp | 8 +- RFL_source/new_source/IAction.hpp | 4 +- RFL_source/new_source/IAlgorithm.hpp | 5 +- RFL_source/new_source/IDiracOperator.hpp | 112 +++++++++ .../new_source/IDiracOperatorDerivatives.cpp | 217 ++++++++++++++++++ .../new_source/IDiracOperatorDerivatives.hpp | 15 ++ RFL_source/new_source/Simulation.hpp | 14 +- .../tests/performance/tBenchmark.cpp | 16 +- RFL_source/new_source/tests/tAction.cpp | 29 +-- RFL_source/new_source/tests/tClifford.cpp | 68 +++--- .../new_source/tests/tDiracOperator.cpp | 137 ++++++++++- RFL_source/new_source/tests/tGslRng.cpp | 10 +- RFL_source/new_source/tests/tHamiltonian.cpp | 25 +- RFL_source/new_source/tests/tMetropolis.cpp | 9 +- RFL_source/new_source/tests/tSimulation.cpp | 22 +- .../original_source/source/CMakeLists.txt | 2 +- RFL_source/original_source/source/HMC.cpp | 12 +- RFL_source/original_source/source/cliff.cpp | 41 ++-- TODO.md | 7 +- examples/Type13/CMakeLists.txt | 8 +- examples/Type13/EigenvalueRecorder.cpp | 107 +++++++++ examples/Type13/EigenvalueRecorder.hpp | 45 ++++ examples/Type13/Type13Metropolis.cpp | 27 ++- examples/mauro_hmc_example.cpp | 2 +- examples/mauro_thesis_example.cpp | 10 +- playground/playground.cpp | 8 +- 45 files changed, 1307 insertions(+), 567 deletions(-) create mode 100644 RFL_source/new_source/IDiracOperator.hpp create mode 100644 RFL_source/new_source/IDiracOperatorDerivatives.cpp create mode 100644 RFL_source/new_source/IDiracOperatorDerivatives.hpp create mode 100644 examples/Type13/EigenvalueRecorder.cpp create mode 100644 examples/Type13/EigenvalueRecorder.hpp diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 5ff8f68..8530f6f 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -4,59 +4,22 @@ on: workflow_dispatch: push: branches: - - "**" - pull_request: - branches: - - main + - "*" jobs: - generate_matrix: - name: Generate test matrix - runs-on: ubuntu-latest - outputs: - matrix_includes: ${{steps.set-matrix.outputs.matrix_includes}} - steps: - - uses: actions/checkout@v3 - - id: set-matrix - shell: bash - run: | - echo "Seed for random selection is = ${{github.run_number}}" - echo "${{github.event_name}}" - case "${{github.event_name}}" in - "pull_request") - includes=$(python3 -u .github/workflows/combination_selection.py "all") - echo "matrix_includes=${includes}" >> $GITHUB_OUTPUT - ;; - "push") - includes=$(python3 -u .github/workflows/combination_selection.py ${{github.run_number}} 3) - echo "matrix_includes=${includes}" >> $GITHUB_OUTPUT - ;; - esac - - check-matrix: - runs-on: ubuntu-latest - needs: generate_matrix - steps: - - name: Check matrix definition - run: | - includes='${{ needs.generate_matrix.outputs.matrix_includes }}' - echo "Json received is = $includes" - echo "Check that the json can be converted using GitHubs fromJSON() without error" - matrix_json='${{fromJson(needs.generate_matrix.outputs.matrix_includes)}}' - build_and_test: - needs: [generate_matrix, check-matrix] strategy: fail-fast: false matrix: - BUILD_TYPE: [Release] - OS: [ubuntu-22.04] - ARMA_VERSION: [11.2.3] - include: ${{fromJSON(needs.generate_matrix.outputs.matrix_includes)}} + BUILD_TYPE: [ Release ] + OS: [ ubuntu-latest, macos-latest ] + ARMA_VERSION: [ 12.6.6 ] - # The CMake configure and build commands are platform-agnostic and should work equally well on Windows or Mac. - # You can convert this to a matrix build if you need cross-platform coverage. - # See: https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix + + # The CMake configure and build commands are platform-agnostic and should + # work equally well on Windows or Mac. You can convert this to a matrix + # build if you need cross-platform coverage. See: + # https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix name: "${{matrix.BUILD_TYPE}}, ${{matrix.OS}}, Armadillo v${{matrix.ARMA_VERSION}}." runs-on: ${{matrix.OS}} @@ -70,17 +33,40 @@ jobs: if: startsWith(matrix.OS,'ubuntu') run: | sudo apt-get update -y - sudo apt-get install -y cmake libopenblas-dev liblapack-dev libarpack2-dev libsuperlu-dev wget libgsl-dev doxygen graphviz + sudo apt-get install -y \ + cmake \ + libopenblas-dev \ + liblapack-dev \ + libarpack2-dev \ + libsuperlu-dev \ + wget \ + libgsl-dev \ + doxygen \ + graphviz \ + libhdf5-serial-dev - name: Install macOS dependencies if: startsWith(matrix.OS,'macos') run: | brew update - brew install cmake doxygen graphviz wget - brew install arpack hdf5 libaec openblas superlu pkg-config gsl + brew install \ + cmake \ + doxygen \ + graphviz \ + wget \ + arpack \ + hdf5 \ + libaec \ + openblas \ + superlu \ + pkg-config \ + gsl || true - name: Download Armadillo from SourceForge - run: mkdir armadillo && cd armadillo && wget -O armadillo.tar.xz https://sourceforge.net/projects/arma/files/armadillo-${{matrix.ARMA_VERSION}}.tar.xz/download + run: | + mkdir armadillo + cd armadillo + wget -O armadillo.tar.xz https://sourceforge.net/projects/arma/files/armadillo-${{matrix.ARMA_VERSION}}.tar.xz/download - name: Unzip Armadillo tarball and create build directory working-directory: armadillo @@ -107,8 +93,10 @@ jobs: sudo cmake --install . - name: Configure CMake - # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. - # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type + # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only + # required if you are using a single-configuration generator such as + # make. See: + # https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{matrix.BUILD_TYPE}} - name: Build the entire project diff --git a/.github/workflows/combination_selection.py b/.github/workflows/combination_selection.py index f008e2f..af4167e 100644 --- a/.github/workflows/combination_selection.py +++ b/.github/workflows/combination_selection.py @@ -23,20 +23,20 @@ def create_new_matrix(subset): def main(random_seed, num_selected, select_all): - build_type = ["Release", "Debug"] + build_type = ["Release"] os_versions = [ - "ubuntu-22.04", + "ubuntu-latest", "ubuntu-20.04", - "macos-12", - "macos-11" + "macos-13", + "macos-latest" ] - armadillo_version = ["10.8.2", "11.2.3"] + armadillo_version = ["10.8.2", "11.4.4"] list_of_lists = [build_type, os_versions, armadillo_version] combinations = [p for p in itertools.product(*list_of_lists)] - if select_all == False: + if not select_all: subset = create_combination_selection(combinations, random_seed, num_selected) else: subset = combinations @@ -67,13 +67,8 @@ def main(random_seed, num_selected, select_all): else: random_seed = int(command_line_arguments[1]) - if (len(command_line_arguments) == 3): + if len(command_line_arguments) == 3: random_seed = int(command_line_arguments[1]) num_selected = int(command_line_arguments[2]) sys.exit(main(random_seed, num_selected, select_all)) - -# Proving that we can recreate the combinations by setting the random.seed -# for i in range(10): -# random.seed(1) -# print(random.sample(combinations, 3)) diff --git a/.github/workflows/compatability_tests.yml b/.github/workflows/compatability_tests.yml index 83612da..d55bf05 100644 --- a/.github/workflows/compatability_tests.yml +++ b/.github/workflows/compatability_tests.yml @@ -2,9 +2,7 @@ name: Compatability Tests on: workflow_dispatch: - # pull_request: - # branches: - # - main + pull_request: jobs: generate_matrix: @@ -32,7 +30,7 @@ jobs: echo "Check that the json can be converted using GitHubs fromJSON() without error" matrix_json='${{fromJson(needs.generate_matrix.outputs.matrix_includes)}}' - build_and_test: + compatability_build_and_test: needs: [generate_matrix, check-matrix] strategy: fail-fast: false @@ -52,14 +50,34 @@ jobs: if: startsWith(matrix.OS,'ubuntu') run: | sudo apt-get update -y - sudo apt-get install -y cmake libopenblas-dev liblapack-dev libarpack2-dev libsuperlu-dev wget libgsl-dev doxygen graphviz + sudo apt-get install -y \ + cmake \ + libopenblas-dev \ + liblapack-dev \ + libarpack2-dev \ + libsuperlu-dev \ + wget \ + libgsl-dev \ + doxygen \ + graphviz \ + libhdf5-serial-dev - name: Install macOS dependencies if: startsWith(matrix.OS,'macos') run: | brew update - brew install cmake doxygen graphviz wget - brew install arpack hdf5 libaec openblas superlu pkg-config gsl + brew install --overwrite\ + cmake \ + doxygen \ + graphviz \ + wget \ + arpack \ + hdf5 \ + libaec \ + openblas \ + superlu \ + pkg-config \ + gsl || true - name: Download Armadillo from SourceForge run: mkdir armadillo && cd armadillo && wget -O armadillo.tar.xz https://sourceforge.net/projects/arma/files/armadillo-${{matrix.ARMA_VERSION}}.tar.xz/download diff --git a/.github/workflows/linter.yml b/.github/workflows/linter.yml index 4771b52..48636b7 100644 --- a/.github/workflows/linter.yml +++ b/.github/workflows/linter.yml @@ -2,8 +2,6 @@ name: C++ Linter on: push: - branches: - - "**" jobs: clang-format: diff --git a/.gitignore b/.gitignore index e058cb1..6dd9fde 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,4 @@ lcov* .DS_Store /venv/ /build-xcode/ +.vscode \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 5297e23..179bc69 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,11 +17,11 @@ if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) - # Let's nicely support folders in IDEs + # Nicely supports folders in IDEs set_property(GLOBAL PROPERTY USE_FOLDERS ON) enable_testing() - add_compile_options(-Wall -Wextra -Wpedantic ) + add_compile_options(-Wall -Wextra -Wpedantic) add_compile_options(-fsanitize=address ) add_link_options(-fsanitize=address ) endif() diff --git a/README.md b/README.md index f60a847..a40004c 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,17 @@ Here are my recommended methods of installing CMake per platform: Once you have installed the required dependencies, the RFL library can be built using CMake. +tl;dr +```bash + git clone https://github.com/pauldruce/RFL.git + cd RFL + mkdir build + cd build + cmake .. + cmake --build . --target all -j 4 + ctest -j 4 +``` + Most C/C++ IDEs will have CMake capabilities, and CMake offers a GUI application to make this process easier. However, to build this project via a terminal, you use the following commands: @@ -57,8 +68,6 @@ However, to build this project via a terminal, you use the following commands: The command should generally be `cmake /path/to/source/files` from within an empty build directory. We assumed the source files (notably, the CMakeLists.txt) are in the parent directory where we are running this command. * The command: `cmake -B ./build .` from the root directory of this project - skipping step 1. This is a different but handy way to create the build files. This command will create a folder `./build` and generate the build files within it. - - 3. __Build the project__ This can be done by manually calling `make`, or it is equivalent in the `build` directory. Or CMake has a handy command which is platform-independent: ```bash cmake --build . --target all diff --git a/RFL_source/README.md b/RFL_source/README.md index 6c46079..138e67b 100644 --- a/RFL_source/README.md +++ b/RFL_source/README.md @@ -1,29 +1,40 @@ # RFL Documentation -This directory contains the source code for the original RFL library created by Mauro D'arcangelo and its modern implementation. + +This directory contains the source code for the original RFL library created by Mauro D'arcangelo and its modern +implementation. The aim is to deprecate Mauro's implementation (Original RFL) for version 1.0 of this project. The code will be preserved under a protected branch on this repository. ## Original RFL + The source code for the original library is in `./original_source`. -The library defines only two classes: `Cliff` and `Geom24`, defined in `/include/Cliff.hpp` and `/include/Geom24.hpp`, respectively. +The library defines only two classes: `Cliff` and `Geom24`, defined in `/include/Cliff.hpp` and `/include/Geom24.hpp`, +respectively. - `Cliff` is responsible for creating the 'gamma matrices' for a specific Clifford module. -The general way to specify a Clifford module is by setting two positive integers $p$ and $q$. + The general way to specify a Clifford module is by setting two positive integers $p$ and $q$. - `Geom24` is the main Class for this library. It is responsible for setting up and running the simulation. The original RFL code inherently uses an action of the form: $$S(D) = g_2* Tr(D^2) + g_4*Tr(D^4) $$ -where $g_2$ and $g_4$ are real numbers. This action is the origin of the name of the class `Geom24` as the action contains the quadratic (D^2) and quartic (D^4) traces of the Dirac operator. - +where $g_2$ and $g_4$ are real numbers. This action is the origin of the name of the class `Geom24` as the action +contains the quadratic (D^2) and quartic (D^4) traces of the Dirac operator. ## New RFL -The new implementation of RFL has taken the original RFL source code and made it more flexible without compromising performance. A significant part of the refactoring was to break down the two behemoth classes into smaller classes, implement a shallow layer of inheritance of abstract classes and migrate memory management to modern C++ techniques (read smart pointers with C++17). + +The new implementation of RFL has taken the original RFL source code and made it more flexible without compromising +performance. A significant part of the refactoring was to break down the two behemoth classes into smaller classes, +implement a shallow layer of inheritance of abstract classes and migrate memory management to modern C++ techniques ( +read smart pointers with C++17). The new source code is found in `./new_source`. This refactoring allows for better testing and more flexibility for library users. -It is constantly being developed, and any improvements or enhancements are welcome. Please leave an issue on the GitHub repo. +It is constantly being developed, and any improvements or enhancements are welcome. Please leave an issue on the GitHub +repo. The high performance that Mauro's original implementation strived to develop has been preserved. -To ensure that this performance does not degrade, a unit test that captures the difference between the original and new implementations. -This test requires that the new implementation is no more than 5% slower to run the same simulation in each implementation. +To ensure that this performance does not degrade, a unit test that captures the difference between the original and new +implementations. +This test requires that the new implementation is no more than 5% slower to run the same simulation in each +implementation. diff --git a/RFL_source/new_source/BarrettGlaser/Action.cpp b/RFL_source/new_source/BarrettGlaser/Action.cpp index 0d831bf..4809ea3 100644 --- a/RFL_source/new_source/BarrettGlaser/Action.cpp +++ b/RFL_source/new_source/BarrettGlaser/Action.cpp @@ -6,25 +6,25 @@ using namespace std; using namespace arma; -Action::Action(double g_2) : m_g_2(g_2), m_g_4(1.0) {} +Action::Action(const double g_2) : m_g_2(g_2), m_g_4(1.0) {} -Action::Action(double g_2, double g_4) : m_g_2(g_2), m_g_4(g_4) {} +Action::Action(const double g_2, const double g_4) : m_g_2(g_2), m_g_4(g_4) {} -double Action::calculateSFromDirac(const DiracOperator& dirac) const { - cx_mat dirac_mat = dirac.getDiracMatrix(); - cx_mat dirac_squared = dirac_mat * dirac_mat; - double trace_dirac_squared = trace(dirac_squared).real(); - double trace_dirac_4 = trace(dirac_squared * dirac_squared).real(); +double Action::calculateSFromDirac(const IDiracOperator& dirac) const { + const cx_mat dirac_mat = dirac.getDiracMatrix(); + const cx_mat dirac_squared = dirac_mat * dirac_mat; + const double trace_dirac_squared = trace(dirac_squared).real(); + const double trace_dirac_4 = trace(dirac_squared * dirac_squared).real(); return m_g_2 * trace_dirac_squared + m_g_4 * trace_dirac_4; } -double Action::calculateS(const DiracOperator& dirac) const { +double Action::calculateS(const IDiracOperator& dirac) const { return m_g_2 * dirac.traceOfDiracSquared() + m_g_4 * dirac.traceOfDirac4(); } -void Action::setParams(double g_2, double g_4) { +void Action::setParams(const double g_2, const double g_4) { this->m_g_2 = g_2; this->m_g_4 = g_4; } -void Action::setG4(double value) { this->m_g_4 = value; } -void Action::setG2(double value) { this->m_g_2 = value; } +void Action::setG4(const double value) { this->m_g_4 = value; } +void Action::setG2(const double value) { this->m_g_2 = value; } diff --git a/RFL_source/new_source/BarrettGlaser/Action.hpp b/RFL_source/new_source/BarrettGlaser/Action.hpp index ce1dfcc..d6bbde0 100644 --- a/RFL_source/new_source/BarrettGlaser/Action.hpp +++ b/RFL_source/new_source/BarrettGlaser/Action.hpp @@ -4,10 +4,8 @@ #ifndef RFL_ACTION_HPP #define RFL_ACTION_HPP -#include "DiracOperator.hpp" #include "IAction.hpp" -#include -#include +#include "IDiracOperator.hpp" /** * Action implements the \f$ S(D) = g_2Tr(D^2) + g_4 Tr(D^4) \f$ action as described in the papers @@ -76,7 +74,7 @@ class Action : public IAction { * an optimised algorithm for the B-G action. * @param dirac A reference to the Dirac operator you want to calculate the Barrett-Glaser action of. */ - double calculateS(const DiracOperator& dirac) const override; + double calculateS(const IDiracOperator& dirac) const override; /** * This calculates the Barrett-Glaser action for the supplied Dirac operator by computing the full @@ -84,7 +82,7 @@ class Action : public IAction { * optimum calculation. Please make use of calculateS(...) above for faster results. * @param dirac A reference to the Dirac operator you want to calculate the Barrett-Glaser action of. */ - double calculateSFromDirac(const DiracOperator& dirac) const; + double calculateSFromDirac(const IDiracOperator& dirac) const; private: double m_g_2, m_g_4; diff --git a/RFL_source/new_source/BarrettGlaser/Hamiltonian.cpp b/RFL_source/new_source/BarrettGlaser/Hamiltonian.cpp index a709000..10a8411 100644 --- a/RFL_source/new_source/BarrettGlaser/Hamiltonian.cpp +++ b/RFL_source/new_source/BarrettGlaser/Hamiltonian.cpp @@ -3,16 +3,17 @@ // #include "Hamiltonian.hpp" +#include "../IDiracOperatorDerivatives.hpp" #include "Action.hpp" #include using namespace std; using namespace arma; -Hamiltonian::Hamiltonian(std::unique_ptr&& action, Integrator integrator, double step_size, std::unique_ptr&& rng) +Hamiltonian::Hamiltonian(std::unique_ptr&& action, const Integrator integrator, const double step_size, std::unique_ptr&& rng) : m_action(std::move(action)), m_integrator(integrator), m_dt(step_size), m_rng(std::move(rng)) { } -double Hamiltonian::updateDirac(const DiracOperator& dirac) const { +double Hamiltonian::updateDirac(const IDiracOperator& dirac) const { const double acceptance_val_per_iter = this->run( dirac, 10, @@ -20,22 +21,20 @@ double Hamiltonian::updateDirac(const DiracOperator& dirac) const { return acceptance_val_per_iter; } -void Hamiltonian::sampleMoments(const DiracOperator& dirac) const { +void Hamiltonian::sampleMoments(const IDiracOperator& dirac) const { auto& mom = dirac.getMomenta(); - auto num_matrices = dirac.getNumMatrices(); - auto mat_dim = dirac.getMatrixDimension(); + const auto num_matrices = dirac.getNumMatrices(); + const auto mat_dim = dirac.getMatrixDimension(); for (int i = 0; i < num_matrices; ++i) { // loop on indices for (int j = 0; j < mat_dim; ++j) { - double x; - x = m_rng->getGaussian(1.0); + const double x = m_rng->getGaussian(1.0); mom[i](j, j) = cx_double(x, 0.); for (int k = j + 1; k < mat_dim; ++k) { - double a, b; - a = m_rng->getGaussian(1.0); - b = m_rng->getGaussian(1.0); + const double a = m_rng->getGaussian(1.0); + const double b = m_rng->getGaussian(1.0); mom[i](j, k) = cx_double(a, b) / sqrt(2.); mom[i](k, j) = cx_double(a, -b) / sqrt(2.); } @@ -43,10 +42,10 @@ void Hamiltonian::sampleMoments(const DiracOperator& dirac) const { } } -double Hamiltonian::calculateK(const DiracOperator& dirac) const { +double Hamiltonian::calculateK(const IDiracOperator& dirac) { double res = 0; - auto& mom = dirac.getMomenta(); - auto num_matrices = dirac.getNumMatrices(); + const auto& mom = dirac.getMomenta(); + const auto num_matrices = dirac.getNumMatrices(); for (int i = 0; i < num_matrices; ++i) res += trace(mom[i] * mom[i]).real(); @@ -54,74 +53,74 @@ double Hamiltonian::calculateK(const DiracOperator& dirac) const { return res / 2; } -double Hamiltonian::calculateH(const DiracOperator& dirac) const { +double Hamiltonian::calculateH(const IDiracOperator& dirac) const { return m_action->calculateS(dirac) + calculateK(dirac); } -void Hamiltonian::leapfrog(const DiracOperator& dirac, +void Hamiltonian::leapfrog(const IDiracOperator& dirac, const int& nt, - double g_2) const { + const double g_2) const { auto& mat = dirac.getMatrices(); auto& mom = dirac.getMomenta(); - auto num_matrices = dirac.getNumMatrices(); + const auto num_matrices = dirac.getNumMatrices(); for (int i = 0; i < num_matrices; ++i) { - mat[i] += (m_dt / 2.) * mom[i]; + mat[i] += m_dt / 2. * mom[i]; for (int j = 0; j < nt - 1; ++j) { - mom[i] += -m_dt * dirac.derDirac24(i, true, g_2); + mom[i] += -m_dt * derDirac24(dirac, i, true, g_2); mat[i] += m_dt * mom[i]; } - mom[i] += -m_dt * dirac.derDirac24(i, true, g_2); - mat[i] += (m_dt / 2.) * mom[i]; + mom[i] += -m_dt * derDirac24(dirac, i, true, g_2); + mat[i] += m_dt / 2. * mom[i]; } } -void Hamiltonian::omelyan(const DiracOperator& dirac, +void Hamiltonian::omelyan(const IDiracOperator& dirac, const int& nt, - double g_2) const { - double xi = 0.1931833; + const double g_2) const { auto& mat = dirac.getMatrices(); auto& mom = dirac.getMomenta(); - auto num_matrices = dirac.getNumMatrices(); + const auto num_matrices = dirac.getNumMatrices(); for (int i = 0; i < num_matrices; ++i) { + constexpr double xi = 0.1931833; mat[i] += xi * m_dt * mom[i]; for (int j = 0; j < nt - 1; ++j) { - mom[i] += -(m_dt / 2.) * dirac.derDirac24(i, true, g_2); + mom[i] += -(m_dt / 2.) * derDirac24(dirac, i, true, g_2); mat[i] += (1 - 2 * xi) * m_dt * mom[i]; - mom[i] += -(m_dt / 2.) * dirac.derDirac24(i, true, g_2); + mom[i] += -(m_dt / 2.) * derDirac24(dirac, i, true, g_2); mat[i] += 2 * xi * m_dt * mom[i]; } - mom[i] += -(m_dt / 2.) * dirac.derDirac24(i, true, g_2); + mom[i] += -(m_dt / 2.) * derDirac24(dirac, i, true, g_2); mat[i] += (1 - 2 * xi) * m_dt * mom[i]; - mom[i] += -(m_dt / 2.) * dirac.derDirac24(i, true, g_2); + mom[i] += -(m_dt / 2.) * derDirac24(dirac, i, true, g_2); mat[i] += xi * m_dt * mom[i]; } } -void Hamiltonian::runDualAverage(const DiracOperator& dirac, +void Hamiltonian::runDualAverage(const IDiracOperator& dirac, const int& nt, const int& iter, const double& target) { // initial (_i) and final (_f) potential2, potential4, kinetic, hamiltonian - auto* en_i = new double[4]; - auto* en_f = new double[4]; + auto en_i = vector(4); + auto en_f = vector(4); // dual averaging variables for dt - const double shr = 0.05; - const double kappa = 0.75; - const int i_0 = 10; double stat = 0; - double mu = log(10 * m_dt); + const double mu = log(10 * m_dt); double log_dt_avg = log(m_dt); // iter repetitions of leapfrog for (int i = 0; i < iter; ++i) { + constexpr int i_0 = 10; + constexpr double kappa = 0.75; + constexpr double shr = 0.05; // if it's not the first iteration set potential to // previous final value, otherwise compute it if (i) { @@ -136,25 +135,22 @@ void Hamiltonian::runDualAverage(const DiracOperator& dirac, stat += target - runDualAveragingCore(dirac, nt, en_i, en_f); // perform dual averaging on dt - double log_dt = mu - stat * sqrt(i + 1) / (shr * (i + 1 + i_0)); + const double log_dt = mu - stat * sqrt(i + 1) / (shr * (i + 1 + i_0)); m_dt = exp(log_dt); - double eta = pow(i + 1, -kappa); + const double eta = pow(i + 1, -kappa); log_dt_avg = eta * log_dt + (1 - eta) * log_dt_avg; } // set dt on its final dual averaged value m_dt = exp(log_dt_avg); - - delete[] en_i; - delete[] en_f; } -double Hamiltonian::run(const DiracOperator& dirac, +double Hamiltonian::run(const IDiracOperator& dirac, const int& num_iterations, const int& iter) const { // initial (_i) and final (_f) potential2, potential4, kinetic, hamiltonian - auto* en_i = new double[4]; - auto* en_f = new double[4]; + auto en_i = vector(4); + auto en_f = vector(4); // return statistic double stat = 0; @@ -175,20 +171,17 @@ double Hamiltonian::run(const DiracOperator& dirac, stat += runCore(dirac, num_iterations, en_i, en_f); } - delete[] en_i; - delete[] en_f; - - return (stat / iter); + return stat / iter; } -double Hamiltonian::run(const DiracOperator& dirac, +double Hamiltonian::run(const IDiracOperator& dirac, const int& nt, const double& dt_min, const double& dt_max, const int& iter) { // initial (_i) and final (_f) potential2, potential4, kinetic, hamiltonian - auto* en_i = new double[4]; - auto* en_f = new double[4]; + auto en_i = vector(4); + auto en_f = vector(4); // return statistic double stat = 0; @@ -209,16 +202,13 @@ double Hamiltonian::run(const DiracOperator& dirac, stat += runCore(dirac, nt, dt_min, dt_max, en_i, en_f); } - delete[] en_i; - delete[] en_f; - - return (stat / iter); + return stat / iter; } -double Hamiltonian::runDualAveragingCore(const DiracOperator& dirac, +double Hamiltonian::runDualAveragingCore(const IDiracOperator& dirac, const int& nt, - double* en_i, - double* en_f) const { + vector& en_i, + vector& en_f) const { // acceptance probability (return value) double e = 1; @@ -226,9 +216,8 @@ double Hamiltonian::runDualAveragingCore(const DiracOperator& dirac, sampleMoments(dirac); // store previous configuration - auto num_matrices = dirac.getNumMatrices(); - // TODO: Convert to smart pointer or vector/array - auto* mat_bk = new cx_mat[num_matrices]; + const auto num_matrices = dirac.getNumMatrices(); + auto mat_bk = vector(num_matrices); auto& mat = dirac.getMatrices(); for (int j = 0; j < num_matrices; j++) { mat_bk[j] = mat[j]; @@ -239,9 +228,9 @@ double Hamiltonian::runDualAveragingCore(const DiracOperator& dirac, en_i[3] = m_action->getG2() * en_i[0] + en_i[1] + en_i[2]; // integration - if (m_integrator == Integrator::LEAPFROG) { + if (m_integrator == LEAPFROG) { leapfrog(dirac, nt, m_action->getG2()); - } else if (m_integrator == Integrator::OMELYAN) { + } else if (m_integrator == OMELYAN) { omelyan(dirac, nt, m_action->getG2()); } @@ -268,7 +257,7 @@ double Hamiltonian::runDualAveragingCore(const DiracOperator& dirac, // now do the standard metropolis test else if (en_f[3] > en_i[3]) { // double r = gsl_rng_uniform(m_engine); - double r = m_rng->getUniform(); + const double r = m_rng->getUniform(); e = exp(en_i[3] - en_f[3]); if (r > e) { @@ -282,15 +271,13 @@ double Hamiltonian::runDualAveragingCore(const DiracOperator& dirac, } } - delete[] mat_bk; - return e; } -double Hamiltonian::runCore(const DiracOperator& dirac, +double Hamiltonian::runCore(const IDiracOperator& dirac, const int& nt, - double* en_i, - double* en_f) const { + vector& en_i, + vector& en_f) const { // acceptance probability (return value) double e = 1; @@ -298,8 +285,8 @@ double Hamiltonian::runCore(const DiracOperator& dirac, sampleMoments(dirac); // store previous configuration - auto num_matrices = dirac.getNumMatrices(); - auto* mat_bk = new cx_mat[num_matrices]; + const auto num_matrices = dirac.getNumMatrices(); + auto mat_bk = vector(num_matrices); auto& mat = dirac.getMatrices(); for (int j = 0; j < num_matrices; j++) { mat_bk[j] = mat[j]; @@ -310,9 +297,9 @@ double Hamiltonian::runCore(const DiracOperator& dirac, en_i[3] = m_action->getG2() * en_i[0] + en_i[1] + en_i[2]; // integration - if (m_integrator == Integrator::LEAPFROG) { + if (m_integrator == LEAPFROG) { leapfrog(dirac, nt, m_action->getG2()); - } else if (m_integrator == Integrator::OMELYAN) { + } else if (m_integrator == OMELYAN) { omelyan(dirac, nt, m_action->getG2()); } @@ -324,7 +311,7 @@ double Hamiltonian::runCore(const DiracOperator& dirac, // metropolis test if (en_f[3] > en_i[3]) { - double r = m_rng->getUniform(); + const double r = m_rng->getUniform(); e = exp(en_i[3] - en_f[3]); if (r > e) { @@ -338,68 +325,61 @@ double Hamiltonian::runCore(const DiracOperator& dirac, } } - delete[] mat_bk; - return e; } -double Hamiltonian::runCoreDebug(const DiracOperator& dirac, +double Hamiltonian::runCoreDebug(const IDiracOperator& dirac, const int& nt) const { // exp(-dH) (return value) - double e; // resample momentum sampleMoments(dirac); // store previous configuration - auto num_matrices = dirac.getNumMatrices(); - auto* mat_bk = new cx_mat[num_matrices]; + const auto num_matrices = dirac.getNumMatrices(); + auto mat_bk = vector(num_matrices); auto& mat = dirac.getMatrices(); for (int j = 0; j < num_matrices; j++) { mat_bk[j] = mat[j]; } // calculate initial hamiltonian - double initial_S = m_action->calculateS(dirac); - double initial_K = calculateK(dirac); - double initial_hamiltonian = initial_S + initial_K; + const double initial_S = m_action->calculateS(dirac); + const double initial_K = calculateK(dirac); + const double initial_hamiltonian = initial_S + initial_K; // integration - if (m_integrator == Integrator::LEAPFROG) { + if (m_integrator == LEAPFROG) { leapfrog(dirac, nt, m_action->getG2()); - } else if (m_integrator == Integrator::OMELYAN) { + } else if (m_integrator == OMELYAN) { omelyan(dirac, nt, m_action->getG2()); } // calculate final hamiltonian - double final_S = m_action->calculateS(dirac); - double final_K = calculateK(dirac); - double final_hamiltonian = final_S + final_K; + const double final_S = m_action->calculateS(dirac); + const double final_K = calculateK(dirac); + const double final_hamiltonian = final_S + final_K; - e = exp(initial_hamiltonian - final_hamiltonian); + const double e = exp(initial_hamiltonian - final_hamiltonian); // metropolis test if (final_hamiltonian > initial_hamiltonian) { - double r = m_rng->getUniform(); - - if (r > e) { + if (const double r = m_rng->getUniform(); r > e) { // restore old configuration for (int j = 0; j < num_matrices; ++j) mat[j] = mat_bk[j]; } } - delete[] mat_bk; - return e; } -double Hamiltonian::runCore(const DiracOperator& dirac, +double Hamiltonian::runCore(const IDiracOperator& dirac, const int& nt, const double& dt_min, const double& dt_max, - double* en_i, - double* en_f) { + vector& en_i, + vector& en_f) { // acceptance probability (return value) double e = 1; @@ -410,8 +390,8 @@ double Hamiltonian::runCore(const DiracOperator& dirac, this->m_dt = dt_min + (dt_max - dt_min) * m_rng->getUniform(); // store previous configuration - auto num_matrices = dirac.getNumMatrices(); - auto* mat_bk = new cx_mat[num_matrices]; + const auto num_matrices = dirac.getNumMatrices(); + auto mat_bk = vector(num_matrices); auto& mat = dirac.getMatrices(); for (int j = 0; j < num_matrices; j++) { mat_bk[j] = mat[j]; @@ -422,9 +402,9 @@ double Hamiltonian::runCore(const DiracOperator& dirac, en_i[3] = m_action->getG2() * en_i[0] + en_i[1] + en_i[2]; // integration - if (m_integrator == Integrator::LEAPFROG) { + if (m_integrator == LEAPFROG) { leapfrog(dirac, nt, m_action->getG2()); - } else if (m_integrator == Integrator::OMELYAN) { + } else if (m_integrator == OMELYAN) { omelyan(dirac, nt, m_action->getG2()); } @@ -436,7 +416,7 @@ double Hamiltonian::runCore(const DiracOperator& dirac, // metropolis test if (en_f[3] > en_i[3]) { - double r = m_rng->getUniform(); + const double r = m_rng->getUniform(); e = exp(en_i[3] - en_f[3]); if (r > e) { @@ -450,12 +430,10 @@ double Hamiltonian::runCore(const DiracOperator& dirac, } } - delete[] mat_bk; - return e; } -void Hamiltonian::setStepSize(double dt) { +void Hamiltonian::setStepSize(const double dt) { this->m_dt = dt; } -void Hamiltonian::setIntegrator(Integrator integrator) { this->m_integrator = integrator; } +void Hamiltonian::setIntegrator(const Integrator integrator) { this->m_integrator = integrator; } diff --git a/RFL_source/new_source/BarrettGlaser/Hamiltonian.hpp b/RFL_source/new_source/BarrettGlaser/Hamiltonian.hpp index 53cc142..6349945 100644 --- a/RFL_source/new_source/BarrettGlaser/Hamiltonian.hpp +++ b/RFL_source/new_source/BarrettGlaser/Hamiltonian.hpp @@ -4,11 +4,11 @@ #ifndef RFL_HAMILTONIAN_HPP #define RFL_HAMILTONIAN_HPP +#include + #include "Action.hpp" #include "IAlgorithm.hpp" #include "IRng.hpp" -#include -#include enum Integrator { LEAPFROG, @@ -16,10 +16,14 @@ enum Integrator { }; /** - * A class to represent the Hamiltonian Monte Carlo algorithm for random non-commutative - * geometries involved in a simulation governed by the Barrett-Glaser action. + * @class Hamiltonian + * @brief Represents a Hamiltonian system for performing simulation. + * + * The Hamiltonian class encapsulates the necessary components for performing simulation + * using Hamiltonian dynamics. It combines an action, an integrator, a step size, and a random + * number generator into a coherent simulation algorithm. */ -class Hamiltonian : public IAlgorithm { +class Hamiltonian final : public IAlgorithm { public: /** * The default constructor for this class has been disabled, please use another constructor. @@ -42,7 +46,7 @@ class Hamiltonian : public IAlgorithm { */ Hamiltonian(std::unique_ptr&& action, Integrator integrator, double step_size, std::unique_ptr&& rng); - double updateDirac(const DiracOperator& dirac) const override; + double updateDirac(const IDiracOperator& dirac) const override; /** * @brief Sets the integrator for the system. @@ -57,68 +61,89 @@ class Hamiltonian : public IAlgorithm { */ void setIntegrator(Integrator integrator); - // TODO: Document + /** + * @brief Get the integrator value. + * + * This function returns the current value of the integrator. + * + * @return The current value of the integrator. + */ Integrator getIntegrator() const { return this->m_integrator; }; - // TODO: Document + /** + * @brief Sets the step size for the simulation. + * + * This function sets the step size, represented by the parameter `dt`, for the simulation. + * The step size determines how different the next DiracOperator in the chain is from the + * starting point. + * + * @param dt The step size for the simulation. + */ void setStepSize(double dt); - // TODO: Document + /** + * @brief Retrieves the step size value. + * + * This function returns the step size value associated with the current + * object. + * + * @return The step size value. + */ double getStepSize() const { return this->m_dt; }; private: std::unique_ptr m_action; - Integrator m_integrator = Integrator::LEAPFROG; + Integrator m_integrator = LEAPFROG; double m_dt; std::unique_ptr m_rng; - // This method seems to be the initialiser for the mom variables in DiracOperator - void sampleMoments(const DiracOperator& dirac) const; - double calculateK(const DiracOperator& dirac) const; - double calculateH(const DiracOperator& dirac) const; + // This method seems to be the initialiser for the mom variables in DiracOperator. + void sampleMoments(const IDiracOperator& dirac) const; + static double calculateK(const IDiracOperator& dirac); + double calculateH(const IDiracOperator& dirac) const; - double run(const DiracOperator& dirac, + double run(const IDiracOperator& dirac, const int& num_iterations, const int& iter) const; - double runDualAveragingCore(const DiracOperator& dirac, + double runDualAveragingCore(const IDiracOperator& dirac, const int& nt, - double* en_i, - double* en_f) const; + std::vector& en_i, + std::vector& en_f) const; - double runCore(const DiracOperator& dirac, + double runCore(const IDiracOperator& dirac, const int& nt, - double* en_i, - double* en_f) const; + std::vector& en_i, + std::vector& en_f) const; - double runCoreDebug(const DiracOperator& dirac, + double runCoreDebug(const IDiracOperator& dirac, const int& nt) const; // The methods below modify the step size "this->dt". - void runDualAverage(const DiracOperator& dirac, + void runDualAverage(const IDiracOperator& dirac, const int& nt, const int& iter, const double& target); - double run(const DiracOperator& dirac, + double run(const IDiracOperator& dirac, const int& nt, const double& dt_min, const double& dt_max, const int& iter); - double runCore(const DiracOperator& dirac, + double runCore(const IDiracOperator& dirac, const int& nt, const double& dt_min, const double& dt_max, - double* en_i, - double* en_f); + std::vector& en_i, + std::vector& en_f); // INTEGRATORS - void leapfrog(const DiracOperator& dirac, + void leapfrog(const IDiracOperator& dirac, const int& nt, double g_2) const; - void omelyan(const DiracOperator& dirac, + void omelyan(const IDiracOperator& dirac, const int& nt, double g_2) const; }; diff --git a/RFL_source/new_source/BarrettGlaser/Metropolis.cpp b/RFL_source/new_source/BarrettGlaser/Metropolis.cpp index f89a1b8..9a06fcd 100644 --- a/RFL_source/new_source/BarrettGlaser/Metropolis.cpp +++ b/RFL_source/new_source/BarrettGlaser/Metropolis.cpp @@ -7,7 +7,7 @@ using namespace std; using namespace arma; -double Metropolis::delta24(const DiracOperator& dirac, +double Metropolis::delta24(const IDiracOperator& dirac, const int& x, const int& row_index, const int& column_index, @@ -15,29 +15,29 @@ double Metropolis::delta24(const DiracOperator& dirac, return m_action->getG2() * delta2(dirac, x, row_index, column_index, z) + delta4(dirac, x, row_index, column_index, z); } -double Metropolis::delta2(const DiracOperator& dirac, +double Metropolis::delta2(const IDiracOperator& dirac, const int& x, const int& row_index, const int& column_index, - const cx_double& z) const { - auto& mat = dirac.getMatrices(); - auto& eps = dirac.getEpsilons(); - auto mat_dim = dirac.getMatrixDimension(); - auto gamma_dim = dirac.getGammaDimension(); + const cx_double& z) { + const auto& mat = dirac.getMatrices(); + const auto& eps = dirac.getEpsilons(); + const auto mat_dim = dirac.getMatrixDimension(); + const auto gamma_dim = dirac.getGammaDimension(); if (row_index != column_index) { return 4. * gamma_dim * mat_dim * (2. * (z * mat[x](column_index, row_index)).real() + norm(z)); } else { - double tr_m = trace(mat[x]).real(); + const double tr_m = trace(mat[x]).real(); return 8. * gamma_dim * z.real() * (mat_dim * (mat[x](row_index, row_index).real() + z.real()) + eps[x] * (tr_m + z.real())); } } -double Metropolis::delta4(const DiracOperator& dirac, +double Metropolis::delta4(const IDiracOperator& dirac, const int& x, const int& row_index, const int& column_index, - const cx_double& z) const { + const cx_double& z) { double res = 0.; auto& omega_table_4 = dirac.getOmegaTable4(); @@ -271,26 +271,26 @@ double Metropolis::delta4(const DiracOperator& dirac, return res; } -double Metropolis::runDualAverage(const DiracOperator& dirac, +double Metropolis::runDualAverage(const IDiracOperator& dirac, const double target) { // initial (_i) and final (_f) action2 and action4 auto* s_i = new double[2]; auto* s_f = new double[2]; - auto mat_dim = dirac.getMatrixDimension(); + const auto mat_dim = dirac.getMatrixDimension(); // calculate length of a sweep in terms of dofs - int nsw = dirac.getNumMatrices() * mat_dim * mat_dim - dirac.getNumAntiHermitianMatrices(); + const int nsw = dirac.getNumMatrices() * mat_dim * mat_dim - dirac.getNumAntiHermitianMatrices(); // dual averaging variables - const double shr = 0.05; - const double kappa = 0.75; - const int i_0 = 10; double stat = 0; - double mu = log(10 * m_scale); + const double mu = log(10 * m_scale); double log_scale_avg = log(m_scale); // iter sweeps of metropolis for (int i = 0; i < m_num_steps; ++i) { for (int j = 0; j < nsw; ++j) { + constexpr int i_0 = 10; + constexpr double kappa = 0.75; + constexpr double shr = 0.05; // set action to previous final value, // unless it's the first iteration if (j) { @@ -304,9 +304,9 @@ double Metropolis::runDualAverage(const DiracOperator& dirac, stat += target - runDualAverageCore(dirac, s_i, s_f); // perform dual averaging - double log_scale = mu - stat * sqrt(i + 1) / (shr * (i + 1 + i_0)); + const double log_scale = mu - stat * sqrt(i + 1) / (shr * (i + 1 + i_0)); m_scale = exp(log_scale); - double eta = pow(i + 1, -kappa); + const double eta = pow(i + 1, -kappa); log_scale_avg = eta * log_scale + (1 - eta) * log_scale_avg; } } @@ -319,13 +319,13 @@ double Metropolis::runDualAverage(const DiracOperator& dirac, return (stat / (m_num_steps * nsw)); } -double Metropolis::run(const DiracOperator& dirac) const { +double Metropolis::run(const IDiracOperator& dirac) const { // initial (_i) and final (_f) action2 and action4 auto* s_i = new double[2]; auto* s_f = new double[2]; - auto mat_dim = dirac.getMatrixDimension(); + const auto mat_dim = dirac.getMatrixDimension(); // calculate length of a sweep in terms of dofs - int nsw = dirac.getNumMatrices() * mat_dim * mat_dim - dirac.getNumAntiHermitianMatrices(); + const int nsw = dirac.getNumMatrices() * mat_dim * mat_dim - dirac.getNumAntiHermitianMatrices(); // return statistic double stat = 0; @@ -353,23 +353,23 @@ double Metropolis::run(const DiracOperator& dirac) const { return (stat / (m_num_steps * nsw)); } -double Metropolis::runDualAverageCore(const DiracOperator& dirac, +double Metropolis::runDualAverageCore(const IDiracOperator& dirac, const double* s_i, double* s_f) const { // acceptance probability double e; - auto num_matrices = dirac.getNumMatrices(); - auto mat_dim = dirac.getMatrixDimension(); + const auto num_matrices = dirac.getNumMatrices(); + const auto mat_dim = dirac.getMatrixDimension(); // metropolis - int x = (int)(num_matrices * m_rng->getUniform()); - int row_index = (int)(mat_dim * m_rng->getUniform()); - int column_index = (int)(mat_dim * m_rng->getUniform()); + const int x = (int)(num_matrices * m_rng->getUniform()); + const int row_index = (int)(mat_dim * m_rng->getUniform()); + const int column_index = (int)(mat_dim * m_rng->getUniform()); double re = 0; - double im = 0; cx_double z; if (row_index != column_index) { + double im = 0; re = m_scale * (-1. + 2. * m_rng->getUniform()); im = m_scale * (-1. + 2. * m_rng->getUniform()); z = cx_double(re, im); @@ -378,9 +378,9 @@ double Metropolis::runDualAverageCore(const DiracOperator& dirac, z = cx_double(re, 0); } - double delta_2 = delta2(dirac, x, row_index, column_index, z); - double delta_4 = delta4(dirac, x, row_index, column_index, z); - double action_delta = delta24(dirac, x, row_index, column_index, z); + const double delta_2 = delta2(dirac, x, row_index, column_index, z); + const double delta_4 = delta4(dirac, x, row_index, column_index, z); + const double action_delta = delta24(dirac, x, row_index, column_index, z); auto& mat = dirac.getMatrices(); // metropolis test @@ -401,7 +401,7 @@ double Metropolis::runDualAverageCore(const DiracOperator& dirac, e = 1; } else { e = exp(-action_delta); - double p = m_rng->getUniform(); + const double p = m_rng->getUniform(); if (e > p) { // update matrix element @@ -424,23 +424,23 @@ double Metropolis::runDualAverageCore(const DiracOperator& dirac, return e; } -double Metropolis::runCore(const DiracOperator& dirac, +double Metropolis::runCore(const IDiracOperator& dirac, const double* s_i, double* s_f) const { // acceptance probability double ret = 0; - auto num_matrices = dirac.getNumMatrices(); - auto mat_dim = dirac.getMatrixDimension(); + const auto num_matrices = dirac.getNumMatrices(); + const auto mat_dim = dirac.getMatrixDimension(); // metropolis - int x = (int)(num_matrices * m_rng->getUniform()); - int row_index = (int)(mat_dim * m_rng->getUniform()); - int column_index = (int)(mat_dim * m_rng->getUniform()); + const int x = (int)(num_matrices * m_rng->getUniform()); + const int row_index = (int)(mat_dim * m_rng->getUniform()); + const int column_index = (int)(mat_dim * m_rng->getUniform()); double re = 0; - double im = 0; cx_double z; if (row_index != column_index) { + double im = 0; re = m_scale * (-1. + 2. * m_rng->getUniform()); im = m_scale * (-1. + 2. * m_rng->getUniform()); z = cx_double(re, im); @@ -449,9 +449,9 @@ double Metropolis::runCore(const DiracOperator& dirac, z = cx_double(re, 0); } - double delta_2 = delta2(dirac, x, row_index, column_index, z); - double delta_4 = delta4(dirac, x, row_index, column_index, z); - double action_delta = m_action->getG2() * delta_2 + delta_4; + const double delta_2 = delta2(dirac, x, row_index, column_index, z); + const double delta_4 = delta4(dirac, x, row_index, column_index, z); + const double action_delta = m_action->getG2() * delta_2 + delta_4; auto& mat = dirac.getMatrices(); // metropolis test @@ -471,8 +471,8 @@ double Metropolis::runCore(const DiracOperator& dirac, // move accepted ret = 1; } else { - double e = exp(-action_delta); - double p = m_rng->getUniform(); + const double e = exp(-action_delta); + const double p = m_rng->getUniform(); if (e > p) { // update matrix element diff --git a/RFL_source/new_source/BarrettGlaser/Metropolis.hpp b/RFL_source/new_source/BarrettGlaser/Metropolis.hpp index 09af216..8ad274c 100644 --- a/RFL_source/new_source/BarrettGlaser/Metropolis.hpp +++ b/RFL_source/new_source/BarrettGlaser/Metropolis.hpp @@ -16,7 +16,7 @@ * algorithm for random non-commutative geometries that are following a simulation * governed by the Barrett-Glaser action. */ -class Metropolis : public IAlgorithm { +class Metropolis final : public IAlgorithm { public: /** * The default constructor has been disabled, please use another constructor. @@ -29,20 +29,20 @@ class Metropolis : public IAlgorithm { * @param num_steps is the number of steps to take for each call to the updateDirac method below. * @param rng is a unique_ptr to an implementation of the abstract rng class IRng. */ - Metropolis(std::unique_ptr&& action, double scale, int num_steps, std::unique_ptr&& rng) + Metropolis(std::unique_ptr&& action, const double scale, const int num_steps, std::unique_ptr&& rng) : m_action(std::move(action)), m_scale(scale), m_num_steps(num_steps), m_rng(std::move(rng)){}; /** * setParams updates the scale, number_of_steps and rng parameters that are passed in as part of * the constructor. */ - void setParams(double scale, int number_of_steps, std::unique_ptr&& rng) { + void setParams(const double scale, const int number_of_steps, std::unique_ptr&& rng) { this->m_scale = scale; this->m_num_steps = number_of_steps; this->m_rng = std::move(rng); } - double updateDirac(const DiracOperator& dirac) const override { + double updateDirac(const IDiracOperator& dirac) const override { return this->run(dirac); } @@ -53,37 +53,37 @@ class Metropolis : public IAlgorithm { std::unique_ptr m_rng; // MMC routine that doesn't perform dual averaging - double run(const DiracOperator& dirac) const; + double run(const IDiracOperator& dirac) const; // MMC routine that performs dual averaging - double runDualAverage(const DiracOperator& dirac, + double runDualAverage(const IDiracOperator& dirac, double target); - double delta24(const DiracOperator& dirac, + double delta24(const IDiracOperator& dirac, const int& x, const int& row_index, const int& column_index, const arma::cx_double& z) const; // TODO: These can be made static or members of DiracOperator - double delta2(const DiracOperator& dirac, - const int& x, - const int& row_index, - const int& column_index, - const arma::cx_double& z) const; + static double delta2(const IDiracOperator& dirac, + const int& x, + const int& row_index, + const int& column_index, + const arma::cx_double& z); // TODO: These can be made static or members of DiracOperator - double delta4(const DiracOperator& dirac, - const int& x, - const int& row_index, - const int& column_index, - const arma::cx_double& z) const; + static double delta4(const IDiracOperator& dirac, + const int& x, + const int& row_index, + const int& column_index, + const arma::cx_double& z); - double runDualAverageCore(const DiracOperator& dirac, + double runDualAverageCore(const IDiracOperator& dirac, const double* s_i, double* s_f) const; - double runCore(const DiracOperator& dirac, + double runCore(const IDiracOperator& dirac, const double* s_i, double* s_f) const; }; diff --git a/RFL_source/new_source/CMakeLists.txt b/RFL_source/new_source/CMakeLists.txt index 96e85c6..cb1030a 100644 --- a/RFL_source/new_source/CMakeLists.txt +++ b/RFL_source/new_source/CMakeLists.txt @@ -1,4 +1,7 @@ -add_library(new_RFL) +add_library(new_RFL + IDiracOperator.hpp + IDiracOperatorDerivatives.cpp + IDiracOperatorDerivatives.hpp) set(new_RFL_HEADERS BarrettGlaser/Action.hpp diff --git a/RFL_source/new_source/Clifford.cpp b/RFL_source/new_source/Clifford.cpp index 4987b26..c4ff5ce 100644 --- a/RFL_source/new_source/Clifford.cpp +++ b/RFL_source/new_source/Clifford.cpp @@ -103,23 +103,23 @@ Clifford::Clifford(int p, int q) : m_p(p), m_q(q), m_dim_gamma(0) { //(1,0) if (m_p == 1 && m_q == 0) { - (*this) = Clifford(3); + *this = Clifford(3); } //(0,1) if (m_p == 0 && m_q == 1) { - (*this) = Clifford(4); + *this = Clifford(4); } //(2,0) if (m_p == 2 && m_q == 0) { - (*this) = Clifford(0); + *this = Clifford(0); } //(1,1) if (m_p == 1 && m_q == 1) { - (*this) = Clifford(2); + *this = Clifford(2); } //(0,2) if (m_p == 0 && m_q == 2) { - (*this) = Clifford(1); + *this = Clifford(1); //any other case } else { initGammas(); @@ -166,7 +166,7 @@ Clifford& Clifford::operator=(const Clifford& clifford_to_copy) { * @param q * @param dec */ -static void decomp(int p, int q, int* dec) { +static void decomp(const int p, const int q, int* dec) { if (p) { if (!(p % 2)) { dec[0] = p / 2; @@ -217,33 +217,31 @@ void Clifford::initGammas() { } } - auto begin = vec.begin(); - auto end = vec.end(); - Clifford c_1 = (*begin); + const auto begin = vec.begin(); + const auto end = vec.end(); + Clifford c_1 = *begin; for (auto iter = begin + 1; iter != end; ++iter) { - c_1 *= (*iter); + c_1 *= *iter; } - (*this) = c_1; + *this = c_1; } // TODO: Go through this and make sure it agrees with Lawson+Michelson Clifford& Clifford::operator*=(const Clifford& clifford_2) { // store C2 frequently used variables - int p_2, q_2, dim_2; - p_2 = clifford_2.getP(); - q_2 = clifford_2.getQ(); - dim_2 = clifford_2.getGammaDimension(); + const int p_2 = clifford_2.getP(); + const int q_2 = clifford_2.getQ(); + const int dim_2 = clifford_2.getGammaDimension(); // temporary variables to avoid overwriting on (*this) - int p, q, dim_gamma; vector gamma; - p = m_p + p_2; - q = m_q + q_2; - dim_gamma = m_dim_gamma * dim_2; + const int p = m_p + p_2; + const int q = m_q + q_2; + const int dim_gamma = m_dim_gamma * dim_2; // start computing product - cx_mat id_2(dim_2, dim_2, fill::eye); + const cx_mat id_2(dim_2, dim_2, fill::eye); gamma.reserve(m_p + m_q); for (int i = 0; i < m_p + m_q; ++i) @@ -252,12 +250,11 @@ Clifford& Clifford::operator*=(const Clifford& clifford_2) { gamma.emplace_back(kron(m_chiral, clifford_2.getGammaAtIndex(i))); // compute chirality - int s_2 = (q_2 - p_2 + 8 * p_2) % 8;// +8*p2 is necessary becase % does not mean modulo for negative numbers - bool s_2_even = (s_2 % 8) % 2 == 0; - if (!s_2_even) { - int s = (m_q - m_p + 8 * m_p) % 8; - cx_mat id_1(m_dim_gamma, m_dim_gamma, fill::eye); - if ((s == 2) || (s == 6)) { + const int s_2 = (q_2 - p_2 + 8 * p_2) % 8;// +8*p2 is necessary becase % does not mean modulo for negative numbers + if (const bool s_2_even = s_2 % 8 % 2 == 0; !s_2_even) { + const int s = (m_q - m_p + 8 * m_p) % 8; + const cx_mat id_1(m_dim_gamma, m_dim_gamma, fill::eye); + if (s == 2 || s == 6) { m_chiral = kron(-1 * id_1, clifford_2.getChiral()); } else { m_chiral = kron(id_1, clifford_2.getChiral()); @@ -275,7 +272,7 @@ Clifford& Clifford::operator*=(const Clifford& clifford_2) { for (const auto& v : gamma) m_gammas.push_back(v); - return (*this); + return *this; } ostream& operator<<(ostream& out, const Clifford& clifford) { @@ -285,7 +282,7 @@ ostream& operator<<(ostream& out, const Clifford& clifford) { } bool areHermitian(const cx_mat& m_1, const cx_mat& m_2) { - return (!(m_2.is_hermitian()) && m_1.is_hermitian()); + return !m_2.is_hermitian() && m_1.is_hermitian(); } void Clifford::sortGammas() { diff --git a/RFL_source/new_source/Clifford.hpp b/RFL_source/new_source/Clifford.hpp index cbe9b95..1c6d590 100644 --- a/RFL_source/new_source/Clifford.hpp +++ b/RFL_source/new_source/Clifford.hpp @@ -40,7 +40,7 @@ class Clifford { * getGammaAtIndex returns a copy of a gamma matrix located in the array at index * i. */ - arma::cx_mat getGammaAtIndex(int i) const { return m_gammas.at(i); } + arma::cx_mat getGammaAtIndex(const int i) const { return m_gammas.at(i); } /** * getChiral returns a copy of the chirality matrix generated from the gamma diff --git a/RFL_source/new_source/DiracOperator.cpp b/RFL_source/new_source/DiracOperator.cpp index 1dd2baf..6b36f03 100644 --- a/RFL_source/new_source/DiracOperator.cpp +++ b/RFL_source/new_source/DiracOperator.cpp @@ -4,18 +4,17 @@ #include "DiracOperator.hpp" #include "Clifford.hpp" +#include using namespace std; using namespace arma; DiracOperator::DiracOperator(int p, int q, int dim) - : m_dim(dim) { + : m_clifford(Clifford(p, q)), m_dim(dim) { int n = p + q; - // create a type (p, q) clifford module - Clifford clifford(p, q); - vector gamma = clifford.getGammaMatrices(); - this->m_gamma_dim = clifford.getGammaDimension(); + vector gamma = m_clifford.getGammaMatrices(); + this->m_gamma_dim = m_clifford.getGammaDimension(); vector herm; herm.reserve(p); @@ -56,7 +55,7 @@ DiracOperator::DiracOperator(int p, int q, int dim) // } bool first = true; for (const auto& v : vec) { - if (first) {// skipp first entry + if (first) {// skip first entry first = false; continue; } @@ -100,7 +99,19 @@ DiracOperator::DiracOperator(int p, int q, int dim) } } -DiracOperator::~DiracOperator() { +DiracOperator::DiracOperator(const DiracOperator& original) + : m_clifford(original.m_clifford) { + // Make a copy of all state of original DiracOperator. + this->m_dim = original.m_dim; + this->m_num_matrices = original.m_num_herm; + this->m_num_herm = original.m_num_herm; + this->m_num_antiherm = original.m_num_antiherm; + this->m_gamma_dim = original.m_gamma_dim; + this->m_matrices = std::make_unique>(*original.m_matrices); + this->m_momenta = std::make_unique>(*original.m_momenta); + this->m_omegas = std::make_unique>(*original.m_omegas); + this->m_epsilons = std::make_unique>(*original.m_epsilons); + this->m_omega_table_4 = std::make_unique>(*original.m_omega_table_4); } // TODO: figure out what this function is for. @@ -130,10 +141,10 @@ static vector baseConversion(int dec, const int& base, const int& max) { cx_mat DiracOperator::getDiracMatrix() const { // initialize dirac op to zero - int dim_dirac = m_dim * m_dim * m_gamma_dim; + const int dim_dirac = m_dim * m_dim * m_gamma_dim; cx_mat dirac(dim_dirac, dim_dirac, fill::zeros); - cx_mat id(m_dim, m_dim, fill::eye); + const cx_mat id(m_dim, m_dim, fill::eye); for (int i = 0; i < m_num_matrices; ++i) { cx_mat bracket = kron((*m_matrices)[i], id) + (*m_epsilons)[i] * kron(id, (*m_matrices)[i].st()); dirac += kron((*m_omegas)[i], bracket); @@ -142,6 +153,33 @@ cx_mat DiracOperator::getDiracMatrix() const { return dirac; } +vec DiracOperator::getEigenvalues() const { + const auto diracMatrix = getDiracMatrix(); + assert(diracMatrix.is_hermitian(1e-16)); + auto eigen_vals = arma::eig_sym(diracMatrix); + return eigen_vals; +} + +vector DiracOperator::getHermitianMatrices() const { + std::vector herm_matrices; + for (std::size_t i = 0; i < m_matrices->size(); i++) { + if (m_epsilons->at(i) == 1) { + herm_matrices.push_back(m_matrices->at(i)); + } + } + return herm_matrices; +} + +vector DiracOperator::getAntiHermitianMatrices() const { + vector anti_herm_matrices; + for (size_t i = 0; i < m_matrices->size(); i++) { + if (m_epsilons->at(i) == -1) { + anti_herm_matrices.emplace_back(m_matrices->at(i) * cx_double(0, 1)); + } + } + return anti_herm_matrices; +} + void DiracOperator::printOmegaTable4() const { const int n = (int)pow(m_num_matrices, 4); @@ -247,11 +285,14 @@ cx_mat DiracOperator::derDirac4(const int& k, const bool& herm) const { // all indices equal res += computeB(k); + cx_mat final_result; if (herm) { - return 2 * (res + res.t()); + final_result = 2 * (res + res.t()); } else { - return 4 * res; + final_result = 4 * res; } + + return final_result; } cx_mat DiracOperator::computeB4(const int& k, @@ -356,17 +397,17 @@ cx_mat DiracOperator::computeB2(const int& k, const int& i) const { } cx_mat DiracOperator::computeB(const int& k) const { - auto& matrices = this->getMatrices(); - auto& epsilons = this->getEpsilons(); + const auto& matrices = this->getMatrices(); + const auto& epsilons = this->getEpsilons(); // base matrix products - cx_mat m_2 = matrices[k] * matrices[k]; - cx_mat m_3 = matrices[k] * m_2; + const cx_mat m_2 = matrices[k] * matrices[k]; + const cx_mat m_3 = matrices[k] * m_2; // traces - double tr_3 = trace(m_3).real(); - double tr_2 = trace(m_2).real(); - double tr_1 = trace(matrices[k]).real(); + const double tr_3 = trace(m_3).real(); + const double tr_2 = trace(m_2).real(); + const double tr_1 = trace(matrices[k]).real(); cx_mat res(m_dim, m_dim, fill::eye); res *= epsilons[k] * tr_3; @@ -385,14 +426,12 @@ void DiracOperator::randomiseMatrices(const IRng& rng) const { for (int i = 0; i < m_num_matrices; ++i) { // loop on indices for (int j = 0; j < m_dim; ++j) { - double x; - x = rng.getGaussian(1.0); + const double x = rng.getGaussian(1.0); matrices[i](j, j) = cx_double(x, 0.); for (int k = j + 1; k < m_dim; ++k) { - double a, b; - a = rng.getGaussian(1.0); - b = rng.getGaussian(1.0); + const double a = rng.getGaussian(1.0); + const double b = rng.getGaussian(1.0); matrices[i](j, k) = cx_double(a, b) / sqrt(2.); matrices[i](k, j) = cx_double(a, -b) / sqrt(2.); } @@ -402,14 +441,12 @@ void DiracOperator::randomiseMatrices(const IRng& rng) const { for (int i = 0; i < m_num_matrices; ++i) { // loop on indices for (int j = 0; j < m_dim; ++j) { - double x; - x = rng.getGaussian(1.0); + const double x = rng.getGaussian(1.0); momenta[i](j, j) = cx_double(x, 0.); for (int k = j + 1; k < m_dim; ++k) { - double a, b; - a = rng.getGaussian(1.0); - b = rng.getGaussian(1.0); + const double a = rng.getGaussian(1.0); + const double b = rng.getGaussian(1.0); momenta[i](j, k) = cx_double(a, b) / sqrt(2.); momenta[i](k, j) = cx_double(a, -b) / sqrt(2.); } @@ -418,13 +455,13 @@ void DiracOperator::randomiseMatrices(const IRng& rng) const { } double DiracOperator::traceOfDiracSquared() const { - auto& matrices = this->getMatrices(); - auto& epsilons = this->getEpsilons(); + const auto& matrices = this->getMatrices(); + const auto& epsilons = this->getEpsilons(); double res = 0.; for (int i = 0; i < m_num_matrices; ++i) { - double trace_squared = trace(matrices[i] * matrices[i]).real(); - double trace_mat = trace(matrices[i]).real(); + const double trace_squared = trace(matrices[i] * matrices[i]).real(); + const double trace_mat = trace(matrices[i]).real(); res += (m_dim * trace_squared + epsilons[i] * trace_mat * trace_mat); } @@ -435,7 +472,7 @@ double DiracOperator::traceOfDiracSquared() const { double DiracOperator::traceOfDirac4() const { double res = 0.; - auto num_matrices = this->getNumMatrices(); + const auto num_matrices = this->getNumMatrices(); // four distinct indices for (int i = 0; i < num_matrices; ++i) { @@ -556,28 +593,28 @@ double DiracOperator::computeA4(const int& i_1, const int& i_2, const int& i_3, } double DiracOperator::computeA2(const int& i_1, const int& i_2) const { - auto& omega_table_4 = this->getOmegaTable4(); - auto& matrices = this->getMatrices(); - auto& epsilons = this->getEpsilons(); + const auto& omega_table_4 = this->getOmegaTable4(); + const auto& matrices = this->getMatrices(); + const auto& epsilons = this->getEpsilons(); // clifford product - double cliff = omega_table_4[i_2 + m_num_matrices * (i_1 + m_num_matrices * (i_2 + m_num_matrices * i_1))].real(); + const double cliff = omega_table_4[i_2 + m_num_matrices * (i_1 + m_num_matrices * (i_2 + m_num_matrices * i_1))].real(); // base matrix products - cx_mat m_1_m_1 = matrices[i_1] * matrices[i_1]; - cx_mat m_2_m_2 = matrices[i_2] * matrices[i_2]; - cx_mat m_1_m_2 = matrices[i_1] * matrices[i_2]; + const cx_mat m_1_m_1 = matrices[i_1] * matrices[i_1]; + const cx_mat m_2_m_2 = matrices[i_2] * matrices[i_2]; + const cx_mat m_1_m_2 = matrices[i_1] * matrices[i_2]; // traces - double tr_1122 = trace(m_1_m_1 * m_2_m_2).real(); - double tr_1212 = trace(m_1_m_2 * m_1_m_2).real(); - double tr_122 = trace(m_1_m_2 * matrices[i_2]).real(); - double tr_112 = trace(m_1_m_1 * matrices[i_2]).real(); - double tr_11 = trace(m_1_m_1).real(); - double tr_22 = trace(m_2_m_2).real(); - double tr_12 = trace(m_1_m_2).real(); - double tr_1 = trace(matrices[i_1]).real(); - double tr_2 = trace(matrices[i_2]).real(); + const double tr_1122 = trace(m_1_m_1 * m_2_m_2).real(); + const double tr_1212 = trace(m_1_m_2 * m_1_m_2).real(); + const double tr_122 = trace(m_1_m_2 * matrices[i_2]).real(); + const double tr_112 = trace(m_1_m_1 * matrices[i_2]).real(); + const double tr_11 = trace(m_1_m_1).real(); + const double tr_22 = trace(m_2_m_2).real(); + const double tr_12 = trace(m_1_m_2).real(); + const double tr_1 = trace(matrices[i_1]).real(); + const double tr_2 = trace(matrices[i_2]).real(); if (cliff < 0) { // compute sum @@ -601,18 +638,18 @@ double DiracOperator::computeA2(const int& i_1, const int& i_2) const { } double DiracOperator::computeA(const int& i) const { - auto& matrices = this->getMatrices(); - auto& epsilons = this->getEpsilons(); + const auto& matrices = this->getMatrices(); + const auto& epsilons = this->getEpsilons(); // base matrix products - cx_mat m_2 = matrices[i] * matrices[i]; - cx_mat m_3 = m_2 * matrices[i]; + const cx_mat m_2 = matrices[i] * matrices[i]; + const cx_mat m_3 = m_2 * matrices[i]; // traces - double tr_1 = trace(matrices[i]).real(); - double tr_2 = trace(m_2).real(); - double tr_3 = trace(m_3).real(); - double tr_4 = trace(m_3 * matrices[i]).real(); + const double tr_1 = trace(matrices[i]).real(); + const double tr_2 = trace(m_2).real(); + const double tr_3 = trace(m_3).real(); + const double tr_4 = trace(m_3 * matrices[i]).real(); double res = m_dim * tr_4; res += 4 * epsilons[i] * tr_1 * tr_3; diff --git a/RFL_source/new_source/DiracOperator.hpp b/RFL_source/new_source/DiracOperator.hpp index b7102d8..281d5ac 100644 --- a/RFL_source/new_source/DiracOperator.hpp +++ b/RFL_source/new_source/DiracOperator.hpp @@ -6,26 +6,22 @@ #define RFL_DIRACOPERATOR_HPP #include "Clifford.hpp" +#include "IDiracOperator.hpp" #include "IRng.hpp" #include #include /** - * A class that represents the Dirac operator as described in the paper: - - * John W. Barrett, "Matrix geometries and fuzzy spaces as finite spectral triples", - * J. Math. Phys. 56, 082301 (2015) https://doi.org/10.1063/1.4927224 + * @class DiracOperator + * + * @brief The DiracOperator class represents a Dirac operator with specified + * Clifford type and dimension of matrices. * + * The DiracOperator class provides methods to access various properties and + * perform computations related to the Dirac operator. */ -class DiracOperator { +class DiracOperator final : public IDiracOperator { public: - /** - * The default constructor for this class has been deleted, see the other options - * that are available - */ - DiracOperator() = delete; - ~DiracOperator(); - /** * The constructor for this class. The parameters p and q represent the * Clifford type for the Dirac operator \f$(p,q)\f$. The parameter dim is the @@ -33,37 +29,50 @@ class DiracOperator { */ DiracOperator(int p, int q, int dim); + /** + * A copy constructor to duplicate a DiracOperator + */ + DiracOperator(const DiracOperator& original); + + /** + * Returns the Clifford type of the Dirac operator - encoded as a pair of integers + * representing the (p,q) values of the underlying Clifford module. + * + * @return std::pair of values (p,q) + */ + std::pair getType() const override { return std::pair{m_clifford.getP(), m_clifford.getQ()}; } + /** * getMatrixDimension returns the dimension of the H and L matrices of the * Dirac operator */ - int getMatrixDimension() const { return m_dim; }; + int getMatrixDimension() const override { return m_dim; }; /** * getGammaDimension returns the dimension of the gamma matrices of the Dirac * operator */ - int getGammaDimension() const { return m_gamma_dim; }; + int getGammaDimension() const override { return m_gamma_dim; }; /** * getNumMatrices returns the total number of H and L matrices */ - int getNumMatrices() const { return m_num_matrices; }; + int getNumMatrices() const override { return m_num_matrices; }; /** * getNumHermitianMatrices returns the number of H matrices */ - int getNumHermitianMatrices() const { return m_num_herm; }; + int getNumHermitianMatrices() const override { return m_num_herm; }; /** * getNumHermitianMatrices returns the number of L matrices */ - int getNumAntiHermitianMatrices() const { return m_num_antiherm; }; + int getNumAntiHermitianMatrices() const override { return m_num_antiherm; }; /** * This method returns a reference to the vector of H and L matrices of the Dirac operator */ - std::vector& getMatrices() const { return *m_matrices; } + std::vector& getMatrices() const override { return *m_matrices; } /** * getEpsilons returns a reference to a vector of +/-1. The sign of the entry relates @@ -71,30 +80,48 @@ class DiracOperator { * the associated matrix is Hermitian. If the value is -1, then the associated matrix * is anti-Hermitian. */ - std::vector& getEpsilons() const { return *m_epsilons; } + std::vector& getEpsilons() const override { return *m_epsilons; } /** * getMomenta returns a reference to a vector of matrices which correspond to.. */ // TODO: This should probably be moved to the Hamiltonian class, or needs good documentation about it's use - std::vector& getMomenta() const { return *m_momenta; } + std::vector& getMomenta() const override { return *m_momenta; } /** * Returns the value of \f$\text{Tr}(D^2)\f$ */ - double traceOfDiracSquared() const; + double traceOfDiracSquared() const override; /** * Returns the value of \f$\text{Tr}(D^4)\f$ * @return */ - double traceOfDirac4() const; + double traceOfDirac4() const override; /** * getDiracMatrix returns a fully constructed matrix representation of the * Dirac operator. */ - arma::cx_mat getDiracMatrix() const; + arma::cx_mat getDiracMatrix() const override; + + /** + * getEigenvalues returns an armadillo vector of the eigenvalues of the Dirac Operator + * when expressed in it's fully assembled matrix form. + */ + arma::vec getEigenvalues() const override; + + /** + * Returns a vector of the Hermitian matrices, H_i, in the decomposition of + * the Dirac Operator. + */ + std::vector getHermitianMatrices() const override; + + /** + * Returns a vector of the Hermitian matrices, H_i, in the decomposition of + * the Dirac Operator. + */ + std::vector getAntiHermitianMatrices() const override; // TODO: Document arma::cx_mat derDirac24(const int& k, const bool& herm, double g_2) const; @@ -103,13 +130,15 @@ class DiracOperator { // TODO: Document arma::cx_mat derDirac4(const int& k, const bool& herm) const; // TODO: Document - std::vector& getOmegaTable4() const { return *m_omega_table_4; } + std::vector& getOmegaTable4() const override { return *m_omega_table_4; } // TODO: Document void printOmegaTable4() const; // TODO: Document - void randomiseMatrices(const IRng& rng_engine) const; + void randomiseMatrices(const IRng& rng) const override; private: + // The clifford module that makes up part of the Dirac operator. + Clifford m_clifford; // CONSTANTS // The dimension of the H and L matrices. int m_dim; @@ -117,7 +146,6 @@ class DiracOperator { int m_num_matrices, m_num_herm, m_num_antiherm; // size of gamma matrices int m_gamma_dim; - // MATRICES // H and L matrices (all hermitian) std::unique_ptr> m_matrices; diff --git a/RFL_source/new_source/GslRng.hpp b/RFL_source/new_source/GslRng.hpp index 6e16326..c0393c4 100644 --- a/RFL_source/new_source/GslRng.hpp +++ b/RFL_source/new_source/GslRng.hpp @@ -19,13 +19,13 @@ * * To do so, a new class that implements the interface defined by IRng needs to be constructed. */ -class GslRng : public IRng { +class GslRng final : public IRng { public: /** * This default constructor creates a GSL random number generator that is * seeded by the current time. */ - GslRng() { + GslRng() : m_rng(nullptr) { gsl_rng_env_setup(); m_rng = gsl_rng_alloc(gsl_rng_ranlxd1); gsl_rng_set(m_rng, time(nullptr)); @@ -35,7 +35,7 @@ class GslRng : public IRng { * This constructor creates a GSL random number generator that is seeded by * the provided parameter 'seed'. */ - explicit GslRng(unsigned long seed) { + explicit GslRng(const unsigned long seed) : m_rng(nullptr) { m_rng = gsl_rng_alloc(gsl_rng_ranlxd1); gsl_rng_set(m_rng, seed); } @@ -46,7 +46,7 @@ class GslRng : public IRng { } } - double getGaussian(double sigma) const override { + double getGaussian(const double sigma) const override { return gsl_ran_gaussian(m_rng, sigma); } double getUniform() const override { diff --git a/RFL_source/new_source/IAction.hpp b/RFL_source/new_source/IAction.hpp index be2d17a..2023ce1 100644 --- a/RFL_source/new_source/IAction.hpp +++ b/RFL_source/new_source/IAction.hpp @@ -5,7 +5,7 @@ #ifndef RFL_IACTION_HPP #define RFL_IACTION_HPP -#include "DiracOperator.hpp" +#include "IDiracOperator.hpp" /** * IAction is an abstract class that provides an interface for any action you may want to use in a @@ -22,7 +22,7 @@ class IAction { * @param dirac is the Dirac operator that the action should be calculated for. * @return the value of the action for the input Dirac operator. */ - virtual double calculateS(const DiracOperator& dirac) const = 0; + virtual double calculateS(const IDiracOperator& dirac) const = 0; virtual ~IAction() = default; }; diff --git a/RFL_source/new_source/IAlgorithm.hpp b/RFL_source/new_source/IAlgorithm.hpp index cacabe8..b4b1356 100644 --- a/RFL_source/new_source/IAlgorithm.hpp +++ b/RFL_source/new_source/IAlgorithm.hpp @@ -5,8 +5,7 @@ #ifndef RFL_IALGORITHM_HPP #define RFL_IALGORITHM_HPP -#include "DiracOperator.hpp" -#include "IAction.hpp" +#include "IDiracOperator.hpp" /** * This is an Abstract Class/Interface for the various Monte Carlo @@ -21,7 +20,7 @@ class IAlgorithm { * @param dirac is the DiracOperator the algorithm will operate on. * @return the acceptance rate of the process. */ - virtual double updateDirac(const DiracOperator& dirac) const = 0; + virtual double updateDirac(const IDiracOperator& dirac) const = 0; virtual ~IAlgorithm() = default; }; diff --git a/RFL_source/new_source/IDiracOperator.hpp b/RFL_source/new_source/IDiracOperator.hpp new file mode 100644 index 0000000..40ce599 --- /dev/null +++ b/RFL_source/new_source/IDiracOperator.hpp @@ -0,0 +1,112 @@ +// +// Created by Paul Druce on 16/11/2023. +// + +#ifndef IDIRACOPERATOR_HPP +#define IDIRACOPERATOR_HPP + +#include "IRng.hpp" +#include + +class IDiracOperator { +public: + virtual ~IDiracOperator() = default; + + // TODO: Document + virtual void randomiseMatrices(const IRng& rng_engine) const = 0; + + /** + * Returns the Clifford type of the Dirac operator - encoded as a pair of integers + * representing the (p,q) values of the underlying Clifford module. + * + * @return std::pair of values (p,q) + */ + virtual std::pair getType() const = 0; + + /** + * getMatrixDimension returns the dimension of the H and L matrices of the + * Dirac operator + */ + virtual int getMatrixDimension() const = 0; + + /** + * getGammaDimension returns the dimension of the gamma matrices of the Dirac + * operator + */ + virtual int getGammaDimension() const = 0; + + /** + * getNumMatrices returns the total number of H and L matrices + */ + virtual int getNumMatrices() const = 0; + + /** + * getNumHermitianMatrices returns the number of H matrices + */ + virtual int getNumHermitianMatrices() const = 0; + + /** + * getNumHermitianMatrices returns the number of L matrices + */ + virtual int getNumAntiHermitianMatrices() const = 0; + + /** + * This method returns a reference to the vector of H and L matrices of the Dirac operator + */ + virtual std::vector& getMatrices() const = 0; + + /** + * getEpsilons returns a reference to a vector of +/-1. The sign of the entry relates + * to the matrix with the same index returned from getMatrices. If the value is +1 then + * the associated matrix is Hermitian. If the value is -1, then the associated matrix + * is anti-Hermitian. + */ + virtual std::vector& getEpsilons() const = 0; + + /** + * getMomenta returns a reference to a vector of matrices which correspond to.. + */ + // TODO: This should probably be moved to the Hamiltonian class, or needs good documentation about it's use + virtual std::vector& getMomenta() const = 0; + + /** + * Returns the value of \f$\text{Tr}(D^2)\f$ + */ + virtual double traceOfDiracSquared() const = 0; + + /** + * Returns the value of \f$\text{Tr}(D^4)\f$ + * @return + */ + virtual double traceOfDirac4() const = 0; + + /** + * getDiracMatrix returns a fully constructed matrix representation of the + * Dirac operator. + */ + virtual arma::cx_mat getDiracMatrix() const = 0; + + /** + * getEigenvalues returns an armadillo vector of the eigenvalues of the Dirac Operator + * when expressed in it's fully assembled matrix form. + */ + virtual arma::vec getEigenvalues() const = 0; + + /** + * Returns a vector of the Hermitian matrices, H_i, in the decomposition of + * the Dirac Operator. + */ + virtual std::vector getHermitianMatrices() const = 0; + + /** + * Returns a vector of the Hermitian matrices, H_i, in the decomposition of + * the Dirac Operator. + */ + virtual std::vector getAntiHermitianMatrices() const = 0; + + // TODO: Document what this should do. + // This feels like an implementation detail, but some other methods require an access to this + // for the moment. + virtual std::vector& getOmegaTable4() const = 0; +}; +#endif//IDIRACOPERATOR_HPP diff --git a/RFL_source/new_source/IDiracOperatorDerivatives.cpp b/RFL_source/new_source/IDiracOperatorDerivatives.cpp new file mode 100644 index 0000000..ae201ba --- /dev/null +++ b/RFL_source/new_source/IDiracOperatorDerivatives.cpp @@ -0,0 +1,217 @@ +// +// Created by Paul Druce on 17/11/2023. +// + +#include "IDiracOperatorDerivatives.hpp" + +using namespace arma; + +static cx_mat computeB4(const IDiracOperator& dirac, const int& k, + const int& i_2, + const int& i_3, + const int& i_4, + const double& cliff, + const bool& neg) { + auto& matrices = dirac.getMatrices(); + auto& epsilons = dirac.getEpsilons(); + const auto dim = dirac.getMatrixDimension(); + + // base matrix products + cx_mat m_2_m_3 = matrices[i_2] * matrices[i_3]; + cx_mat m_2_m_4 = matrices[i_2] * matrices[i_4]; + cx_mat m_3_m_4 = matrices[i_3] * matrices[i_4]; + cx_mat m_2_m_3_m_4 = m_2_m_3 * matrices[i_4]; + + // return value + cx_mat res(dim, dim, fill::eye); + + if (neg) { + + // traces + double tr_234 = trace(m_2_m_3_m_4).imag(); + double tr_2 = trace(matrices[i_2]).real(); + double tr_3 = trace(matrices[i_3]).real(); + double tr_4 = trace(matrices[i_4]).real(); + + // compute sum + res *= -2 * epsilons[k] * tr_234; + res += cx_double(0., dim) * (m_2_m_3_m_4 - m_2_m_3_m_4.t()); + res += cx_double(0., epsilons[i_2] * tr_2) * (m_3_m_4 - m_3_m_4.t()); + res += cx_double(0., epsilons[i_3] * tr_3) * (m_2_m_4 - m_2_m_4.t()); + res += cx_double(0., epsilons[i_4] * tr_4) * (m_2_m_3 - m_2_m_3.t()); + } else { + // traces + double tr_234 = trace(m_2_m_3_m_4).real(); + double tr_23 = trace(m_2_m_3).real(); + double tr_24 = trace(m_2_m_4).real(); + double tr_34 = trace(m_3_m_4).real(); + double tr_2 = trace(matrices[i_2]).real(); + double tr_3 = trace(matrices[i_3]).real(); + double tr_4 = trace(matrices[i_4]).real(); + + // compute sum + res *= 2 * epsilons[k] * tr_234; + res += dim * (m_2_m_3_m_4 + m_2_m_3_m_4.t()); + res += epsilons[i_2] * tr_2 * (m_3_m_4 + m_3_m_4.t()); + res += epsilons[i_3] * tr_3 * (m_2_m_4 + m_2_m_4.t()); + res += epsilons[i_4] * tr_4 * (m_2_m_3 + m_2_m_3.t()); + res += 2 * epsilons[k] * epsilons[i_2] * tr_34 * matrices[i_2]; + res += 2 * epsilons[k] * epsilons[i_3] * tr_24 * matrices[i_3]; + res += 2 * epsilons[k] * epsilons[i_4] * tr_23 * matrices[i_4]; + } + + return cliff * res; +} + +static cx_mat computeB2(const IDiracOperator& dirac, const int& k, const int& i) { + const auto num_matrices = dirac.getNumMatrices(); + const auto dim = dirac.getMatrixDimension(); + const auto gamma_dim = dirac.getGammaDimension(); + + auto& omega_table_4 = dirac.getOmegaTable4(); + auto& matrices = dirac.getMatrices(); + auto& epsilons = dirac.getEpsilons(); + + // clifford product + double cliff = omega_table_4[i + num_matrices * (k + num_matrices * (i + num_matrices * k))].real(); + + // base matrix products + cx_mat mi_mk = matrices[i] * matrices[k]; + cx_mat mi_mi = matrices[i] * matrices[i]; + cx_mat mi_mi_mk = matrices[i] * mi_mk; + cx_mat mi_mk_mi = mi_mk * matrices[i]; + + // traces + double triki = trace(mi_mk_mi).real(); + double trik = trace(mi_mk).real(); + double trii = trace(mi_mi).real(); + double tri = trace(matrices[i]).real(); + double trk = trace(matrices[k]).real(); + + // return value + cx_mat res(dim, dim, fill::eye); + + if (cliff < 0) { + // compute sum + res *= epsilons[k] * triki; + res += dim * (mi_mi_mk + mi_mi_mk.t() - mi_mk_mi); + res += epsilons[i] * tri * (mi_mk + mi_mk.t()); + res += 2 * epsilons[k] * epsilons[i] * trik * matrices[i]; + res += epsilons[k] * trk * mi_mi; + res += trii * matrices[k]; + } else { + // compute sum + res *= 3 * epsilons[k] * triki; + res += dim * (mi_mi_mk + mi_mi_mk.t() + mi_mk_mi); + res += 3 * epsilons[i] * tri * (mi_mk + mi_mk.t()); + res += 6 * epsilons[k] * epsilons[i] * trik * matrices[i]; + res += 3 * epsilons[k] * trk * mi_mi; + res += 3 * trii * matrices[k]; + } + + return 2 * gamma_dim * res; +} + +static cx_mat computeB(const IDiracOperator& dirac, const int& k) { + const auto m_dim = dirac.getMatrixDimension(); + const auto m_gamma_dim = dirac.getGammaDimension(); + const auto& matrices = dirac.getMatrices(); + const auto& epsilons = dirac.getEpsilons(); + + // base matrix products + const cx_mat m_2 = matrices[k] * matrices[k]; + const cx_mat m_3 = matrices[k] * m_2; + + // traces + const double tr_3 = trace(m_3).real(); + const double tr_2 = trace(m_2).real(); + const double tr_1 = trace(matrices[k]).real(); + + cx_mat res(m_dim, m_dim, fill::eye); + res *= epsilons[k] * tr_3; + res += m_dim * m_3; + res += 3 * tr_2 * matrices[k]; + res += 3 * epsilons[k] * tr_1 * m_2; + + return 2 * m_gamma_dim * res; +} + +cx_mat derDirac24(const IDiracOperator& dirac, const int& k, const bool& herm, const double g_2) { + return g_2 * derDirac2(dirac, k) + derDirac4(dirac, k, herm); +} + +arma::cx_mat derDirac2(const IDiracOperator& dirac, const int& k) { + const auto dim = dirac.getMatrixDimension(); + const auto gamma_dim = dirac.getGammaDimension(); + const auto& epsilons = dirac.getEpsilons(); + const auto& matrices = dirac.getMatrices(); + arma::cx_mat res(dim, dim, fill::eye); + res *= epsilons[k] * trace(matrices[k]).real(); + res += dim * matrices[k]; + return 4 * gamma_dim * res; +} + +cx_mat derDirac4(const IDiracOperator& dirac, const int& k, const bool& herm) { + const auto dim = dirac.getMatrixDimension(); + const auto num_matrices = dirac.getNumMatrices(); + const auto& epsilons = dirac.getEpsilons(); + const auto& omega_table_4 = dirac.getOmegaTable4(); + + cx_mat res(dim, dim, fill::zeros); + + // four distinct indices + for (int i_1 = 0; i_1 < num_matrices; ++i_1) { + if (i_1 != k) { + for (int i_2 = i_1 + 1; i_2 < num_matrices; ++i_2) { + if (i_2 != k) { + for (int i_3 = i_2 + 1; i_3 < num_matrices; ++i_3) { + if (i_3 != k) { + // epsilon factor + + if (double e = epsilons[k] * epsilons[i_1] * epsilons[i_2] * epsilons[i_3]; e < 0) { + // clifford product + double cliff_1 = omega_table_4[i_3 + num_matrices * (i_2 + num_matrices * (i_1 + num_matrices * k))].imag(); + double cliff_2 = omega_table_4[i_2 + num_matrices * (i_3 + num_matrices * (i_1 + num_matrices * k))].imag(); + double cliff_3 = omega_table_4[i_3 + num_matrices * (i_1 + num_matrices * (i_2 + num_matrices * k))].imag(); + + if (fabs(cliff_1) > 1e-10) { + res += computeB4(dirac, k, i_1, i_2, i_3, cliff_1, true); + res += computeB4(dirac, k, i_1, i_3, i_2, cliff_2, true); + res += computeB4(dirac, k, i_2, i_1, i_3, cliff_3, true); + } + } else { + // clifford product + double cliff_1 = omega_table_4[i_3 + num_matrices * (i_2 + num_matrices * (i_1 + num_matrices * k))].real(); + double cliff_2 = omega_table_4[i_2 + num_matrices * (i_3 + num_matrices * (i_1 + num_matrices * k))].real(); + double cliff_3 = omega_table_4[i_3 + num_matrices * (i_1 + num_matrices * (i_2 + num_matrices * k))].real(); + + if (fabs(cliff_1) > 1e-10) { + res += computeB4(dirac, k, i_1, i_2, i_3, cliff_1, false); + res += computeB4(dirac, k, i_1, i_3, i_2, cliff_2, false); + res += computeB4(dirac, k, i_2, i_1, i_3, cliff_3, false); + } + } + } + } + } + } + } + } + res = res + res.t(); + + // two distinct pairs of equal indices + for (int i = 0; i < num_matrices; ++i) { + if (i != k) { + res += computeB2(dirac, k, i); + } + } + + // all indices equal + res += computeB(dirac, k); + + if (herm) { + return 2 * (res + res.t()); + } else { + return 4 * res; + } +} diff --git a/RFL_source/new_source/IDiracOperatorDerivatives.hpp b/RFL_source/new_source/IDiracOperatorDerivatives.hpp new file mode 100644 index 0000000..fca45e3 --- /dev/null +++ b/RFL_source/new_source/IDiracOperatorDerivatives.hpp @@ -0,0 +1,15 @@ +// +// Created by Paul Druce on 17/11/2023. +// + +#ifndef IDIRACOPERATORDERIVATIVES_HPP +#define IDIRACOPERATORDERIVATIVES_HPP +#include "IDiracOperator.hpp" + +arma::cx_mat derDirac24(const IDiracOperator& dirac, const int& k, const bool& herm, double g_2); + +arma::cx_mat derDirac2(const IDiracOperator& dirac, const int& k); + +arma::cx_mat derDirac4(const IDiracOperator& dirac, const int& k, const bool& herm); + +#endif//IDIRACOPERATORDERIVATIVES_HPP diff --git a/RFL_source/new_source/Simulation.hpp b/RFL_source/new_source/Simulation.hpp index ecef752..a35fbfc 100644 --- a/RFL_source/new_source/Simulation.hpp +++ b/RFL_source/new_source/Simulation.hpp @@ -6,9 +6,7 @@ #define RFL_RFL_SOURCE_NEW_SOURCE_SIMULATION_HPP_ #include "DiracOperator.hpp" -#include "IAction.hpp" #include "IAlgorithm.hpp" -#include /** * A class to control the simulation. This class is responsible for the lifespan of @@ -33,8 +31,16 @@ class Simulation { /** * This method starts the simulation and will return when it is complete. */ - void run() { - this->m_algorithm->updateDirac(*m_dirac); + double run() const { + return this->m_algorithm->updateDirac(*m_dirac); + } + + /** + * This method returns a reference to DiracOperator in use in this class. + * @return const reference to DiracOperator object + */ + const DiracOperator& getDiracOperator() const { + return *m_dirac; } private: diff --git a/RFL_source/new_source/tests/performance/tBenchmark.cpp b/RFL_source/new_source/tests/performance/tBenchmark.cpp index d8f6e9d..8d204e8 100644 --- a/RFL_source/new_source/tests/performance/tBenchmark.cpp +++ b/RFL_source/new_source/tests/performance/tBenchmark.cpp @@ -40,7 +40,7 @@ void NewMethod() { auto metropolis = std::make_unique(std::move(action), SCALE, NUM_STEPS, std::move(rng)); - auto simulation = Simulation( + const auto simulation = Simulation( std::move(dirac), std::move(metropolis)); @@ -70,25 +70,25 @@ TEST(BenchmarkTests, NewImplementationIsNotSignificantlySlower) { using std::chrono::high_resolution_clock; using std::chrono::milliseconds; - auto mauro_start = high_resolution_clock::now(); + const auto mauro_start = high_resolution_clock::now(); MauroMethod(); - auto mauro_stop = high_resolution_clock::now(); + const auto mauro_stop = high_resolution_clock::now(); /* Getting number of milliseconds as an integer. */ - auto mauro_ms_int = duration_cast(mauro_stop - mauro_start); + const auto mauro_ms_int = duration_cast(mauro_stop - mauro_start); std::cout << "Mauro's implementation:\n"; std::cout << mauro_ms_int.count() << "ms\n"; - auto new_start = high_resolution_clock::now(); + const auto new_start = high_resolution_clock::now(); NewMethod(); - auto new_stop = high_resolution_clock::now(); + const auto new_stop = high_resolution_clock::now(); - auto new_ms_int = duration_cast(new_stop - new_start); + const auto new_ms_int = duration_cast(new_stop - new_start); std::cout << "New implementation:\n"; std::cout << new_ms_int.count() << "ms\n"; - auto winner = new_ms_int.count() < mauro_ms_int.count() ? "the new implementation" : "Mauro's implementation"; + const auto winner = new_ms_int.count() < mauro_ms_int.count() ? "the new implementation" : "Mauro's implementation"; std::cout << "The winner is: " << winner << std::endl; ASSERT_TRUE(new_ms_int.count() < mauro_ms_int.count() * MAX_FRACTION_DIFF) << "New implementation is slower by more than 5%"; diff --git a/RFL_source/new_source/tests/tAction.cpp b/RFL_source/new_source/tests/tAction.cpp index 4157ef1..412ea3f 100644 --- a/RFL_source/new_source/tests/tAction.cpp +++ b/RFL_source/new_source/tests/tAction.cpp @@ -2,18 +2,19 @@ // Created by Paul Druce on 27/12/2022. // #include "BarrettGlaser/Action.hpp" +#include "DiracOperator.hpp" #include "GslRng.hpp" #include -static void CompareActions(int p, int q, int dim, double g_2) { +static void CompareActions(const int p, const int q, const int dim, const double g_2) { constexpr int num_of_test_repeats = 100; - Action action(g_2); - DiracOperator dirac(p, q, dim); - GslRng rng; + const Action action(g_2); + const DiracOperator dirac(p, q, dim); for (int i = 0; i < num_of_test_repeats; ++i) { + const GslRng rng; dirac.randomiseMatrices(rng); - double d_2 = dirac.getMatrixDimension() * dirac.getMatrixDimension(); + const double d_2 = dirac.getMatrixDimension() * dirac.getMatrixDimension(); auto s_1 = action.calculateS(dirac) / d_2; @@ -41,39 +42,39 @@ TEST(ActionTests, ActionMethodsDontDiffer) { {2, 1, 4, -0.1}, {0, 5, 4, -2.8}}; - for (auto& d : testing_params) { - CompareActions(d.p, d.q, d.dim, d.g_2); + for (const auto& [p, q, dim, g_2] : testing_params) { + CompareActions(p, q, dim, g_2); } } TEST(ActionTests, Set_g2) { - const double initial_g_2 = 1.1; + constexpr double initial_g_2 = 1.1; Action action(initial_g_2); EXPECT_EQ(action.getG2(), initial_g_2); - const double new_g_2 = -3.4; + constexpr double new_g_2 = -3.4; action.setG2(new_g_2); EXPECT_EQ(action.getG2(), new_g_2); } TEST(ActionTests, Set_g4) { - const double initial_g_4 = 1.1; + constexpr double initial_g_4 = 1.1; Action action(1.0, initial_g_4); EXPECT_EQ(action.getG4(), initial_g_4); - const double new_g_4 = -3.4; + constexpr double new_g_4 = -3.4; action.setG4(new_g_4); EXPECT_EQ(action.getG4(), new_g_4); } TEST(ActionTests, SetParameters) { - const double initial_g_2 = 1.1, initial_g_4 = 2.2; + constexpr double initial_g_2 = 1.1, initial_g_4 = 2.2; Action action(initial_g_2, initial_g_4); EXPECT_EQ(action.getG2(), initial_g_2); EXPECT_EQ(action.getG4(), initial_g_4); - const double new_g_2 = -5.0, new_g_4 = 7.2; + constexpr double new_g_2 = -5.0, new_g_4 = 7.2; action.setParams(new_g_2, new_g_4); EXPECT_EQ(action.getG2(), new_g_2); @@ -82,7 +83,7 @@ TEST(ActionTests, SetParameters) { TEST(ActionTests, CreateWithNoParams) { EXPECT_NO_THROW( - Action action; + const Action action; EXPECT_EQ(action.getG2(), 0); EXPECT_EQ(action.getG4(), 0);); } \ No newline at end of file diff --git a/RFL_source/new_source/tests/tClifford.cpp b/RFL_source/new_source/tests/tClifford.cpp index b185b38..2cabf7c 100644 --- a/RFL_source/new_source/tests/tClifford.cpp +++ b/RFL_source/new_source/tests/tClifford.cpp @@ -14,7 +14,7 @@ typedef struct { } CliffordData; TEST(CliffordTests, NoErrorWhenConstructing) { - std::vector cliff_data = { + const std::vector cliff_data = { {1, 0}, {0, 1}, {1, 1}, @@ -26,9 +26,9 @@ TEST(CliffordTests, NoErrorWhenConstructing) { {2, 3}, {3, 3}}; - for (auto& d : cliff_data) { + for (const auto& [p, q] : cliff_data) { EXPECT_NO_THROW( - Clifford cliff(d.p, d.q);); + Clifford cliff(p, q);); } } @@ -54,7 +54,7 @@ TEST(CliffordTests, GetDimGamma) { int dim_gamma; } DimGammaData; - std::vector data = { + const std::vector data = { {1, 0, 1}, {0, 1, 1}, {1, 1, 2}, @@ -63,9 +63,9 @@ TEST(CliffordTests, GetDimGamma) { {1, 3, 4}, {3, 1, 4}, {3, 3, 8}}; - for (auto& d : data) { - Clifford cliff(d.p, d.q); - EXPECT_EQ(d.dim_gamma, cliff.getGammaDimension()) << "(p,q) = (" << d.p << "," << d.q << ")"; + for (const auto& [p, q, dim_gamma] : data) { + Clifford cliff(p, q); + EXPECT_EQ(dim_gamma, cliff.getGammaDimension()) << "(p,q) = (" << p << "," << q << ")"; } } @@ -76,7 +76,7 @@ TEST(CliffordTests, GetGammas) { int num_gammas; } GetGammasData; - std::vector data = { + const std::vector data = { {1, 0, 1}, {0, 1, 1}, {1, 1, 2}, @@ -85,9 +85,9 @@ TEST(CliffordTests, GetGammas) { {1, 3, 4}, {3, 3, 6}}; - for (auto& d : data) { - Clifford cliff(d.p, d.q); - EXPECT_EQ(d.num_gammas, cliff.getGammaMatrices().size()) << "(p,q) = (" << d.p << "," << d.q << ")"; + for (const auto& [p, q, num_gammas] : data) { + Clifford cliff(p, q); + EXPECT_EQ(num_gammas, cliff.getGammaMatrices().size()) << "(p,q) = (" << p << "," << q << ")"; } } @@ -101,21 +101,21 @@ TEST(CliffordTests, GammasHaveCorrectHermiticity) { } } - for (const auto& d : data) { - Clifford C(d.p, d.q); + for (const auto& [p, q] : data) { + Clifford C(p, q); auto gammas = C.getGammaMatrices(); - std::vector herm_gammas(gammas.begin(), gammas.begin() + d.p); - std::vector anti_herm_gammas(gammas.begin() + d.p, gammas.end()); + std::vector herm_gammas(gammas.begin(), gammas.begin() + p); + std::vector anti_herm_gammas(gammas.begin() + p, gammas.end()); - EXPECT_EQ(d.p, herm_gammas.size()); + EXPECT_EQ(p, herm_gammas.size()); for (auto& hg : herm_gammas) { - EXPECT_TRUE(hg.is_hermitian()) << "(p,q) = (" << d.p << "," << d.q << ")\n "; + EXPECT_TRUE(hg.is_hermitian()) << "(p,q) = (" << p << "," << q << ")\n "; } - EXPECT_EQ(d.q, anti_herm_gammas.size()); + EXPECT_EQ(q, anti_herm_gammas.size()); for (auto& ahg : anti_herm_gammas) { - EXPECT_TRUE(!ahg.is_hermitian()) << "(p,q) = (" << d.p << "," << d.q << ")"; + EXPECT_TRUE(!ahg.is_hermitian()) << "(p,q) = (" << p << "," << q << ")"; } } } @@ -131,11 +131,11 @@ TEST(CliffordTests, GammasHaveCorrectDims) { } } - for (const auto& d : data) { - Clifford C(d.p, d.q); - int n = d.p + d.q; - int exponent = (n % 2 == 0) ? (int)n / 2 : (int)((n - 1) / 2); - int expect_dim = 1 << exponent; + for (const auto& [p, q] : data) { + Clifford C(p, q); + const int n = p + q; + const int exponent = (n % 2 == 0) ? (int)n / 2 : (int)((n - 1) / 2); + const int expect_dim = 1 << exponent; // Easy check EXPECT_EQ(expect_dim, C.getGammaDimension()); @@ -144,8 +144,8 @@ TEST(CliffordTests, GammasHaveCorrectDims) { auto gammas = C.getGammaMatrices(); for (auto& g : gammas) { ASSERT_EQ(g.n_rows, g.n_cols) - << "Gamma matrices not squarefailed for (p,q) = (" << d.p << "," << d.q << ")"; - EXPECT_EQ(expect_dim, g.n_rows) << "Gamma dim not correct for (p,q) = (" << d.p << "," << d.q << ")"; + << "Gamma matrices not squarefailed for (p,q) = (" << p << "," << q << ")"; + EXPECT_EQ(expect_dim, g.n_rows) << "Gamma dim not correct for (p,q) = (" << p << "," << q << ")"; } } } @@ -154,20 +154,20 @@ TEST(CliffordTests, GammasHaveCorrectDims) { TEST(CliffordTests, ChiralityIsCorrect) { // TODO: Fix source code, because this test fails if max_p, max_q is >= 6 constexpr int max_p = 5; - constexpr int max_q = 5; vector data = {}; data.push_back({1, 0}); data.push_back({0, 1}); for (int q = 1; q < max_p; q++) { + constexpr int max_q = 5; for (int p = 1; p < max_q; p++) { data.push_back({p, q}); } } - for (auto& d : data) { - Clifford C(d.p, d.q); - int s = (d.q - d.p + 8 * d.p) % 8;// Need to add 8*p because % does behave for negative values. + for (const auto& [p, q] : data) { + Clifford C(p, q); + const int s = (q - p + 8 * p) % 8;// Need to add 8*p because % does behave for negative values. auto gammas = C.getGammaMatrices(); auto chirality = C.getChiral(); @@ -191,7 +191,7 @@ TEST(CliffordTests, ChiralityIsCorrect) { TEST(CliffordTests, CliffordModuleInlineMultiplication) { Clifford C1(1, 2); - Clifford C2(2, 1); + const Clifford C2(2, 1); C1 *= C2; @@ -200,10 +200,10 @@ TEST(CliffordTests, CliffordModuleInlineMultiplication) { } TEST(CliffordTests, CliffordModuleMultiplication) { - Clifford C1(1, 2); - Clifford C2(2, 1); + const Clifford C1(1, 2); + const Clifford C2(2, 1); - auto C3 = C1 * C2; + const auto C3 = C1 * C2; // New Clifford module has correct dimensions EXPECT_EQ(C3.getP(), 3); diff --git a/RFL_source/new_source/tests/tDiracOperator.cpp b/RFL_source/new_source/tests/tDiracOperator.cpp index 63f9965..f44efba 100644 --- a/RFL_source/new_source/tests/tDiracOperator.cpp +++ b/RFL_source/new_source/tests/tDiracOperator.cpp @@ -25,15 +25,15 @@ TEST(DiracOperatorTests, NumHermitianMatricesIsCorrectlySet) { int num_h_mat; } NumHMatData; - std::vector data = { + const std::vector data = { {1, 1, 1}, {1, 2, 1}, {2, 1, 3}, {3, 3, 16}}; - for (auto& d : data) { - DiracOperator dirac(d.p, d.q, 5); - EXPECT_EQ(d.num_h_mat, dirac.getNumHermitianMatrices()) << "(p,q) = (" << d.p << "," << d.q << ")"; + for (const auto& [p, q, num_h_mat] : data) { + DiracOperator dirac(p, q, 5); + EXPECT_EQ(num_h_mat, dirac.getNumHermitianMatrices()) << "(p,q) = (" << p << "," << q << ")"; } } @@ -44,15 +44,15 @@ TEST(DiracOperatorTests, NumAntiHermitianMatricesIsCorrectlySet) { int num_h_mat; } NumLMatData; - std::vector data = { + const std::vector data = { {1, 1, 1}, {1, 2, 3}, {2, 1, 1}, {3, 3, 16}}; - for (auto& d : data) { - DiracOperator dirac(d.p, d.q, 5); - EXPECT_EQ(d.num_h_mat, dirac.getNumAntiHermitianMatrices()) << "(p,q) = (" << d.p << "," << d.q << ")"; + for (const auto& [p, q, num_h_mat] : data) { + DiracOperator dirac(p, q, 5); + EXPECT_EQ(num_h_mat, dirac.getNumAntiHermitianMatrices()) << "(p,q) = (" << p << "," << q << ")"; } } @@ -63,14 +63,127 @@ TEST(DiracOperatorTests, TotalNumMatricesIsSetCorrectly) { int num_hl_mat; } NumTotalMatricesData; - std::vector data = { + const std::vector data = { {1, 1, 2}, {1, 2, 4}, {2, 1, 4}, {3, 3, 32}}; - for (auto& d : data) { - DiracOperator dirac(d.p, d.q, 5); - EXPECT_EQ(d.num_hl_mat, dirac.getNumMatrices()) << "(p,q) = (" << d.p << "," << d.q << ")"; + for (const auto& [p, q, num_hl_mat] : data) { + DiracOperator dirac(p, q, 5); + EXPECT_EQ(num_hl_mat, dirac.getNumMatrices()) << "(p,q) = (" << p << "," << q << ")"; + } +} + +TEST(DiracOperatorTests, GetType) { + const std::vector> data = { + {1, 1}, + {1, 2}, + {2, 1}, + {3, 3}}; + + for (const auto& [p, q] : data) { + DiracOperator dirac(p, q, 5); + EXPECT_EQ(p, dirac.getType().first); + EXPECT_EQ(q, dirac.getType().second); + } +} + +TEST(DiracOperatorTests, GetEigenvalues) { + const std::vector> data = { + {1, 1}, + {1, 2}, + {2, 1}, + {3, 3}}; + + for (const auto& [p, q] : data) { + DiracOperator dirac(p, q, 5); + + auto eigenvalues = dirac.getEigenvalues(); + + auto dirac_matrix = dirac.getDiracMatrix(); + auto expected_eigenvalues = arma::eig_sym(dirac_matrix); + + ASSERT_EQ(eigenvalues.n_elem, expected_eigenvalues.n_elem); + + eigenvalues = arma::sort(eigenvalues); + expected_eigenvalues = arma::sort(expected_eigenvalues); + + for (unsigned i = 0; i < (unsigned)eigenvalues.n_elem; i++) { + EXPECT_FLOAT_EQ(eigenvalues[i], expected_eigenvalues[i]); + } + } +} + +TEST(DiracOperatorTests, CopyConstructorWorks) { + const DiracOperator dirac_1(1, 3, 5); + const DiracOperator dirac_2 = dirac_1;// NOLINT(*-unnecessary-copy-initialization) + + ASSERT_TRUE(arma::approx_equal( + dirac_1.getDiracMatrix(), + dirac_2.getDiracMatrix(), + "absdiff", + 1e-10)); +} + +TEST(DiracOperatorTests, GetHermitianMatricesReturnsHermitianMatrices) { + typedef struct { + int p; + int q; + int dim; + int n_herm_matrices; + } DiracParams; + + const std::vector params{ + {2, 1, 5, 3}, + {1, 2, 6, 1}, + {3, 3, 7, 16}}; + + for (const auto& [p, q, dim, n_herm_matrices] : params) { + DiracOperator dirac(p, q, dim); + + auto hermitian_matrices = dirac.getHermitianMatrices(); + EXPECT_EQ(n_herm_matrices, hermitian_matrices.size()) + << "For dirac params" + << "(p,q)=(" << p << "," << q << ")" + << " dim = " << dim; + + for (auto& matrix : hermitian_matrices) { + EXPECT_TRUE(matrix.is_hermitian()) + << "For dirac params" + << "(p,q)=(" << p << "," << q << ")" + << " dim = " << dim; + } + } +} + +TEST(DiracOperatorTests, GetAntiHermitianMatricesReturnsAntiHermitianMatrices) { + typedef struct { + int p; + int q; + int dim; + int n_anti_herm_matrices; + } DiracParams; + + const std::vector params{ + {2, 1, 5, 1}, + {1, 2, 6, 3}, + {3, 3, 7, 16}}; + + for (const auto& [p, q, dim, n_anti_herm_matrices] : params) { + DiracOperator dirac(p, q, dim); + + auto hermitian_matrices = dirac.getAntiHermitianMatrices(); + EXPECT_EQ(n_anti_herm_matrices, hermitian_matrices.size()) + << "For dirac params" + << "(p,q)=(" << p << "," << q << ")" + << " dim = " << dim; + + for (auto& matrix : hermitian_matrices) { + EXPECT_TRUE(!matrix.is_hermitian()) + << "For dirac params" + << "(p,q)=(" << p << "," << q << ")" + << " dim = " << dim; + } } } \ No newline at end of file diff --git a/RFL_source/new_source/tests/tGslRng.cpp b/RFL_source/new_source/tests/tGslRng.cpp index bb052e4..531f89c 100644 --- a/RFL_source/new_source/tests/tGslRng.cpp +++ b/RFL_source/new_source/tests/tGslRng.cpp @@ -6,7 +6,7 @@ #include TEST(GslRngTests, SettingSeedProducesSameValues) { - GslRng rng(1); + const GslRng rng(1); double expected_values[5] = {0.83451879245814453, 0.61670202724383927, @@ -14,22 +14,22 @@ TEST(GslRngTests, SettingSeedProducesSameValues) { 0.038517618778346474, 0.5896974345675261}; std::cout.precision(17); - for (double expected_value : expected_values) { + for (const double expected_value : expected_values) { EXPECT_DOUBLE_EQ(rng.getUniform(), expected_value); } } TEST(GslRngTests, NoSeedProducesDifferentValues) { - GslRng rng; // Get first 5 numbers double first_values[10]; for (double& first_value : first_values) { + const GslRng rng; first_value = rng.getUniform(); } - GslRng new_rng; - for (double orig_val : first_values) { + for (const double orig_val : first_values) { + const GslRng new_rng; EXPECT_LE(std::abs(orig_val - new_rng.getUniform()), 1e-8); } } \ No newline at end of file diff --git a/RFL_source/new_source/tests/tHamiltonian.cpp b/RFL_source/new_source/tests/tHamiltonian.cpp index f9944d5..08bfbca 100644 --- a/RFL_source/new_source/tests/tHamiltonian.cpp +++ b/RFL_source/new_source/tests/tHamiltonian.cpp @@ -3,51 +3,52 @@ // #include "BarrettGlaser/Hamiltonian.hpp" +#include "DiracOperator.hpp" #include "GslRng.hpp" #include TEST(HamiltonianTests, ConstructorDoesNotThrow) { - const Integrator integrator = LEAPFROG; + constexpr Integrator integrator = LEAPFROG; // TODO: How do we test random stuff? -> find out - double step_size = 0.1; + constexpr double step_size = 0.1; auto action = std::make_unique(1.0, 1.0); ASSERT_NO_THROW( const Hamiltonian hamiltonian(std::move(action), integrator, step_size, std::make_unique());); } TEST(HamiltonianTests, CanChangeIntegrator) { - const Integrator integrator = LEAPFROG; - double step_size = 0.1; + constexpr Integrator integrator = LEAPFROG; + constexpr double step_size = 0.1; auto action = std::make_unique(1.0, 1.0); Hamiltonian hamiltonian(std::move(action), integrator, step_size, std::make_unique()); ASSERT_EQ(hamiltonian.getIntegrator(), integrator); - const Integrator new_integrator = OMELYAN; + constexpr Integrator new_integrator = OMELYAN; hamiltonian.setIntegrator(new_integrator); ASSERT_EQ(hamiltonian.getIntegrator(), new_integrator); } TEST(HamiltonianTests, CanChangeStepSize) { - const Integrator integrator = LEAPFROG; - double step_size = 0.1; + constexpr Integrator integrator = LEAPFROG; + constexpr double step_size = 0.1; auto action = std::make_unique(1.0, 1.0); Hamiltonian hamiltonian(std::move(action), integrator, step_size, std::make_unique()); ASSERT_EQ(hamiltonian.getStepSize(), step_size); - const double new_step_size = 0.5; + constexpr double new_step_size = 0.5; hamiltonian.setStepSize(new_step_size); ASSERT_EQ(hamiltonian.getStepSize(), new_step_size); } TEST(HamiltonianTests, UpdateDiracUpdatesTheDirac) { - Hamiltonian hamiltonian(std::make_unique(), Integrator::LEAPFROG, 0.2, std::make_unique()); - auto dirac = DiracOperator(1, 1, 5); - auto old_dirac_matrix = dirac.getDiracMatrix(); + const Hamiltonian hamiltonian(std::make_unique(), Integrator::LEAPFROG, 0.2, std::make_unique()); + const auto dirac = DiracOperator(1, 1, 5); + const auto old_dirac_matrix = dirac.getDiracMatrix(); hamiltonian.updateDirac(dirac); - auto new_dirac_matrix = dirac.getDiracMatrix(); + const auto new_dirac_matrix = dirac.getDiracMatrix(); const auto diracs_are_equal = arma::approx_equal(new_dirac_matrix, old_dirac_matrix, "absdiff", 1e-6); diff --git a/RFL_source/new_source/tests/tMetropolis.cpp b/RFL_source/new_source/tests/tMetropolis.cpp index af57de1..3c14f0c 100644 --- a/RFL_source/new_source/tests/tMetropolis.cpp +++ b/RFL_source/new_source/tests/tMetropolis.cpp @@ -3,6 +3,7 @@ // #include "BarrettGlaser/Metropolis.hpp" +#include "DiracOperator.hpp" #include "GslRng.hpp" #include @@ -23,12 +24,12 @@ TEST(MetropolisTests, UpdateDiracUpdatesTheDirac) { constexpr int num_steps = 20; auto action = std::make_unique(1.0, 1.0); - Metropolis metropolis(std::move(action), scale, num_steps, std::move(rng)); + const Metropolis metropolis(std::move(action), scale, num_steps, std::move(rng)); - auto dirac = DiracOperator(1, 1, 5); - auto old_dirac_matrix = dirac.getDiracMatrix(); + const auto dirac = DiracOperator(1, 1, 5); + const auto old_dirac_matrix = dirac.getDiracMatrix(); metropolis.updateDirac(dirac); - auto new_dirac_matrix = dirac.getDiracMatrix(); + const auto new_dirac_matrix = dirac.getDiracMatrix(); const auto diracs_are_equal = arma::approx_equal(new_dirac_matrix, old_dirac_matrix, "absdiff", 1e-6); diff --git a/RFL_source/new_source/tests/tSimulation.cpp b/RFL_source/new_source/tests/tSimulation.cpp index 93fc3af..c9114c3 100644 --- a/RFL_source/new_source/tests/tSimulation.cpp +++ b/RFL_source/new_source/tests/tSimulation.cpp @@ -8,9 +8,9 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" -class MockAlgorithm : public IAlgorithm { +class MockAlgorithm final : public IAlgorithm { public: - MOCK_METHOD(double, updateDirac, (const DiracOperator& dirac), (const override)); + MOCK_METHOD(double, updateDirac, (const IDiracOperator& dirac), (const override)); }; TEST(SimulationTests, ConstructorDoesNotThrow) { @@ -25,9 +25,25 @@ TEST(SimulationTests, RunCallsUpdateDirac) { EXPECT_CALL(*algo, updateDirac(testing::_)) .Times(1); - Simulation simulation( + const Simulation simulation( std::make_unique(1, 1, 5), std::move(algo)); simulation.run(); } + +TEST(SimulationTests, GetDiracReturnsSameDirac) { + auto dirac = DiracOperator(1, 1, 5); + auto dirac_ptr = std::make_unique(dirac); + const Simulation simulation( + std::move(dirac_ptr), + std::make_unique()); + + const DiracOperator& retrieved_dirac = simulation.getDiracOperator(); + ASSERT_EQ(retrieved_dirac.getType(), std::pair(1, 1)); + ASSERT_TRUE(arma::approx_equal( + retrieved_dirac.getDiracMatrix(), + dirac.getDiracMatrix(), + "absdiff", + 1e-10)); +} diff --git a/RFL_source/original_source/source/CMakeLists.txt b/RFL_source/original_source/source/CMakeLists.txt index 257f226..1218c02 100644 --- a/RFL_source/original_source/source/CMakeLists.txt +++ b/RFL_source/original_source/source/CMakeLists.txt @@ -24,4 +24,4 @@ target_link_libraries(RFL ${ARMADILLO_LIBRARIES} ${GSL_LIBRARIES}) install(TARGETS RFL DESTINATION ${CMAKE_SOURCE_DIR}/lib/RFL/bin/${CMAKE_BUILD_TYPE} PUBLIC_HEADER DESTINATION ${CMAKE_SOURCE_DIR}/lib/RFL/include - ) +) diff --git a/RFL_source/original_source/source/HMC.cpp b/RFL_source/original_source/source/HMC.cpp index 8df68af..7280d79 100644 --- a/RFL_source/original_source/source/HMC.cpp +++ b/RFL_source/original_source/source/HMC.cpp @@ -18,15 +18,15 @@ void Geom24::HMC_duav(const int& Nt, auto* en_f = new double[4]; // dual averaging variables for dt - const double shr = 0.05; - const double kappa = 0.75; - const int i0 = 10; double Stat = 0; - double mu = log(10 * dt); + const double mu = log(10 * dt); double log_dt_avg = log(dt); // iter repetitions of leapfrog for (int i = 0; i < iter; ++i) { + constexpr int i0 = 10; + constexpr double kappa = 0.75; + constexpr double shr = 0.05; // if it's not the first interation set potential to // previous final value, otherwise compute it if (i) { @@ -41,9 +41,9 @@ void Geom24::HMC_duav(const int& Nt, Stat += target - HMC_duav_core(Nt, dt, engine, en_i, en_f, integrator); // perform dual averaging on dt - double log_dt = mu - Stat * sqrt(i + 1) / (shr * (i + 1 + i0)); + const double log_dt = mu - Stat * sqrt(i + 1) / (shr * (i + 1 + i0)); dt = exp(log_dt); - double eta = pow(i + 1, -kappa); + const double eta = pow(i + 1, -kappa); log_dt_avg = eta * log_dt + (1 - eta) * log_dt_avg; } diff --git a/RFL_source/original_source/source/cliff.cpp b/RFL_source/original_source/source/cliff.cpp index f060c7c..fae7eab 100644 --- a/RFL_source/original_source/source/cliff.cpp +++ b/RFL_source/original_source/source/cliff.cpp @@ -97,7 +97,7 @@ Cliff::Cliff(int mode) { } Cliff::Cliff(int p_, int q_) - : p(p_), q(q_) { + : p(p_), q(q_), dim_gamma() { //(1,0) if (p == 1 && q == 0) { (*this) = Cliff(3); @@ -154,7 +154,7 @@ Cliff& Cliff::operator=(const Cliff& C) { return *this; } -static void decomp(int p, int q, int* dec) { +static void decomp(const int p, const int q, int* dec) { if (p) { if (!(p % 2)) { dec[0] = p / 2; @@ -201,34 +201,31 @@ void Cliff::init_gamma() { vec.emplace_back(i); } - auto begin = vec.begin(); - auto end = vec.end(); + const auto begin = vec.begin(); + const auto end = vec.end(); - Cliff C1 = (*begin); + Cliff C1 = *begin; for (auto iter = begin + 1; iter != end; ++iter) - C1 *= (*iter); + C1 *= *iter; - (*this) = C1; + *this = C1; } Cliff& Cliff::operator*=(const Cliff& C2) { // store C2 frequently used variables - int p2, q2, dim2; - p2 = C2.get_p(); - q2 = C2.get_q(); - dim2 = C2.get_dim_gamma(); + const int p2 = C2.get_p(); + const int q2 = C2.get_q(); + const int dim2 = C2.get_dim_gamma(); // temporary variables to avoid overwriting on (*this) - int p_, q_, dim_gamma_; vector gamma_; - - p_ = p + p2; - q_ = q + q2; - dim_gamma_ = dim_gamma * dim2; + const int p_ = p + p2; + const int q_ = q + q2; + const int dim_gamma_ = dim_gamma * dim2; // start computing product - cx_mat id2(dim2, dim2, fill::eye); + const cx_mat id2(dim2, dim2, fill::eye); gamma_.reserve(p + q); for (int i = 0; i < p + q; ++i) @@ -237,11 +234,11 @@ Cliff& Cliff::operator*=(const Cliff& C2) { gamma_.emplace_back(kron(chiral, C2.get_gamma(i))); // compute chirality - int s2 = (q2 - p2 + 8 * p2) % 8;// +8*p2 is necessary becase % does not mean modulo for negative numbers + const int s2 = (q2 - p2 + 8 * p2) % 8;// +8*p2 is necessary becase % does not mean modulo for negative numbers if ((s2 % 8) % 2) { - int s = (q - p + 8 * p) % 8; - cx_mat id1(dim_gamma, dim_gamma, fill::eye); - if ((s == 2) || (s == 6)) { + const int s = (q - p + 8 * p) % 8; + const cx_mat id1(dim_gamma, dim_gamma, fill::eye); + if (s == 2 || s == 6) { chiral = kron(-1 * id1, C2.get_chiral()); } else { chiral = kron(id1, C2.get_chiral()); @@ -260,7 +257,7 @@ Cliff& Cliff::operator*=(const Cliff& C2) { for (const auto& v : gamma_) gamma.push_back(v); - return (*this); + return *this; } ostream& operator<<(ostream& out, const Cliff& C) { diff --git a/TODO.md b/TODO.md index 0846ecd..a8698d8 100644 --- a/TODO.md +++ b/TODO.md @@ -3,14 +3,17 @@ - [ ] Add documentation string for each class and public method - [ ] Add documentation generation to github actions - [ ] Transfer methods that can be made static in Action to more appropriate place -- [ ] Investigate the use of multithreading to improve performance. - - Potentially can just use Armadillo here, rather than do the multithreading myself. +- [ ] Investigate multithreading to improve performance. + - It might be possible to use Armadillo here, rather than do the multithreading myself. - [ ] Replace raw pointer use to use smart pointers and references. - [x] DiracOperator - [x] Action - [ ] Metropolis - [ ] Hamiltonian - [ ] Refactor the github actions to have reusable setup actions for linux, windows, macos. +- [ ] Create github actions to create built libraries for macOS, Windows and linux OS. +- [ ] Make an interface for DiracOperator, rather than using concrete DiracOperator in + Simulation class and Algorithm. ## Completed diff --git a/examples/Type13/CMakeLists.txt b/examples/Type13/CMakeLists.txt index c6dcfde..9ad2e17 100644 --- a/examples/Type13/CMakeLists.txt +++ b/examples/Type13/CMakeLists.txt @@ -1,4 +1,8 @@ -add_executable(Type13Metropolis Type13Metropolis.cpp) +find_package(HDF5 REQUIRED) +add_executable(Type13Metropolis Type13Metropolis.cpp + EigenvalueRecorder.hpp + EigenvalueRecorder.cpp) + +target_link_libraries(Type13Metropolis new_RFL HDF5::HDF5) target_include_directories(Type13Metropolis PUBLIC ${new_RFL_INCUDE_DIR}) -target_link_libraries(Type13Metropolis new_RFL) \ No newline at end of file diff --git a/examples/Type13/EigenvalueRecorder.cpp b/examples/Type13/EigenvalueRecorder.cpp new file mode 100644 index 0000000..d55fa12 --- /dev/null +++ b/examples/Type13/EigenvalueRecorder.cpp @@ -0,0 +1,107 @@ +// +// Created by Paul Druce on 30/12/2023. +// + +#include "EigenvalueRecorder.hpp" +#include +#include + +namespace fs = std::filesystem; + +struct diracData { + int p; + int q; + int matrixSize; +}; + +static std::string create_eigenvalues_dataset(const int diracId) { + std::ostringstream ss; + ss << "/eigenvalues_" << diracId; + return ss.str(); +} + +static std::string create_simulation_group(const std::string& timeString) { + std::ostringstream ss; + ss << "Simulation_" << timeString; + return ss.str(); +} + +static std::string create_eigenvalues_hdf_filename(const diracData& data, const double g2) { + std::ostringstream ss; + + ss << -1 * g2; + std::string g2_string = ss.str(); + ss.str(""); + ss.clear(); + + // Replace decimal with underscore in g2 string. + std::replace(g2_string.begin(), g2_string.end(), '.', '_'); + + ss << "Eigenvalues_" + << data.p << "_" << data.q + << "_N_" << data.matrixSize + << std::fixed << std::setprecision(3) << "_g2_" << g2_string + << ".h5"; + auto eigenvalues_hdf_filename = ss.str(); + return eigenvalues_hdf_filename; +} + +static std::string create_eigenvalues_dir(const diracData& data, const std::string& outputDirectoryPath) { + std::ostringstream ss; + ss << "eigenvalues_" << data.p << "_" << data.q << "_N_" << data.matrixSize << "/"; + auto eigenvalues_folder = outputDirectoryPath + ss.str(); + if (!fs::exists(eigenvalues_folder)) { + fs::create_directory(eigenvalues_folder); + } + return eigenvalues_folder; +} + +static std::string create_output_dir(const diracData& data, const std::string& outputRootPath) { + std::ostringstream ss; + + // Create output dir if needed. + const auto output_path = outputRootPath + "/output/"; + if (!fs::exists(output_path)) { + fs::create_directory(output_path); + } + + // Create type folder if needed. + ss.str(""); + ss.clear(); + ss << output_path << "Type" << data.p << data.q << "/"; + const std::string diracTypeDir = ss.str(); + if (!fs::exists(diracTypeDir)) { + fs::create_directory(diracTypeDir); + } + return ss.str(); +} + +void EigenvalueRecorder::recordEigenvalues(int diracId) const { + auto [p, q] = m_dirac.getType(); + auto matrixSize = m_dirac.getMatrixDimension(); + struct diracData data = {p, q, matrixSize}; + + std::cout << "Generating and saving eigenvalues for " + << "{p,q} = {" << p << ", " << q << "} " + << "and N = " << matrixSize << "\n"; + + auto eigenvalues = m_dirac.getEigenvalues(); + + auto outputRootPathEnv = std::getenv("RFL_OUTPUT_DIR"); + std::string outputRootPath = outputRootPathEnv ? outputRootPathEnv : ""; + if (outputRootPath.empty()) { + outputRootPath = "/tmp/RFL"; + } + if (!fs::exists(outputRootPath)) { + fs::create_directory(outputRootPath); + } + auto outputDirectoryPath = create_output_dir(data, outputRootPath); + auto eigenvaluesOutputPath = create_eigenvalues_dir(data, outputDirectoryPath); + auto hdf5Filename = create_eigenvalues_hdf_filename(data, m_g2); + auto simulation_hdf5_group = create_simulation_group(m_simulationId); + auto datasetName = create_eigenvalues_dataset(diracId); + + std::string hdf_filepath = eigenvaluesOutputPath + hdf5Filename; + std::cout << "Saving eigenvalues to " << hdf_filepath << '\n'; + eigenvalues.save(arma::hdf5_name(hdf_filepath, simulation_hdf5_group + "/" + datasetName, arma::hdf5_opts::append)); +} diff --git a/examples/Type13/EigenvalueRecorder.hpp b/examples/Type13/EigenvalueRecorder.hpp new file mode 100644 index 0000000..308c89c --- /dev/null +++ b/examples/Type13/EigenvalueRecorder.hpp @@ -0,0 +1,45 @@ +// +// Created by Paul Druce on 30/12/2023. +// + +#ifndef EXAMPLES_EIGENVALUERECORDER_HPP +#define EXAMPLES_EIGENVALUERECORDER_HPP + +#define ARMA_USE_HDF5 +#include "DiracOperator.hpp" + +/** + * @class EigenvalueRecorder + * @brief A simple class to record the eigenvalues of a DiracOperator to a file. + * + * This recorder uses a HDF5 file to store the eigenvalues. + * + * The output directory of the files is controlled by the environment variable + * RFL_OUTPUT_DIR. + * If this variable is not set, then all output is placed in /tmp/RFL + * + * Each recorder writes to a dataset group in the HDF5 file named by the constructor + * parameter "simulationId". Each eigenvalue data set is appended to this dataset group. + * + */ +class EigenvalueRecorder { +public: + EigenvalueRecorder(const DiracOperator& dirac, const double g2, const std::string& simulationId) + : m_dirac(dirac), m_g2(g2), m_simulationId(simulationId) { + } + + /** + * @brief Records the eigenvalues of a given Dirac operator. + * + * The recordEigenvalues function accepts a diracId as input and records the eigenvalues of + * the dirac operator to a file. + */ + void recordEigenvalues(int diracId) const; + +private: + const DiracOperator& m_dirac; + const double m_g2; + std::string m_simulationId; +}; + +#endif//EXAMPLES_EIGENVALUERECORDER_HPP \ No newline at end of file diff --git a/examples/Type13/Type13Metropolis.cpp b/examples/Type13/Type13Metropolis.cpp index 86aa955..ce16372 100644 --- a/examples/Type13/Type13Metropolis.cpp +++ b/examples/Type13/Type13Metropolis.cpp @@ -5,28 +5,45 @@ #include "BarrettGlaser/Action.hpp" #include "BarrettGlaser/Metropolis.hpp" #include "DiracOperator.hpp" +#include "EigenvalueRecorder.hpp" #include "GslRng.hpp" #include "Simulation.hpp" +#include + +using namespace arma; + +static std::string getCurrentDateTime() { + const auto t = std::time(nullptr); + const auto tm = *std::localtime(&t); + std::stringstream ss; + ss << std::put_time(&tm, "%Y%m%d%H%M%S"); + return ss.str(); +} int main() { - double metropolis_scale = 0.2; + double metropolisScale = 0.2; int iter = 10; auto rng = std::make_unique(); - auto dirac = std::make_unique(1, 3, 10); - auto action = std::make_unique(-2.7, 1.0); + auto g2 = -2.7; + auto g4 = 1.0; + auto action = std::make_unique(g2, g4); auto metropolis = std::make_unique( std::move(action), - metropolis_scale, + metropolisScale, iter, std::move(rng)); - auto simulation = Simulation(std::move(dirac), std::move(metropolis)); + const auto simulation = Simulation(std::move(dirac), std::move(metropolis)); + const auto timeStamp = getCurrentDateTime(); for (int i = 0; i < 10; i++) { simulation.run(); + const auto& dirac2 = simulation.getDiracOperator(); + EigenvalueRecorder recorder(dirac2, g2, timeStamp); + recorder.recordEigenvalues(i); } return 0; diff --git a/examples/mauro_hmc_example.cpp b/examples/mauro_hmc_example.cpp index e931c1b..2bc5b39 100644 --- a/examples/mauro_hmc_example.cpp +++ b/examples/mauro_hmc_example.cpp @@ -12,7 +12,7 @@ using namespace arma; int main() { // Initialize the random number generator gsl_rng* engine = gsl_rng_alloc(gsl_rng_ranlxd1); - gsl_rng_set(engine, time(NULL)); + gsl_rng_set(engine, time(nullptr)); Geom24 G(2, 0, 10, -2.7); diff --git a/examples/mauro_thesis_example.cpp b/examples/mauro_thesis_example.cpp index 1f13a46..03b6973 100644 --- a/examples/mauro_thesis_example.cpp +++ b/examples/mauro_thesis_example.cpp @@ -12,18 +12,18 @@ int main() { gsl_rng* engine = gsl_rng_alloc(gsl_rng_ranlxd1); gsl_rng_set(engine, time(nullptr)); // Clifford module parameters - int p = 2; - int q = 0; + constexpr int p = 2; + constexpr int q = 0; // Matrix algebra dimension - int n = 32; + constexpr int n = 32; // Coupling constant value - double g2 = -3; + constexpr double g2 = -3; // Metropolis scale factor - double scale = 0.05; // Create the Dirac operator Geom24 G(p, q, n, g2); // Metropolis simulation for (int i = 0; i < 100; ++i) { + constexpr double scale = 0.05; // Metropolis evolution for 100 steps G.MMC(scale, 100, engine); // Print the value of the action diff --git a/playground/playground.cpp b/playground/playground.cpp index ddae083..6dac331 100644 --- a/playground/playground.cpp +++ b/playground/playground.cpp @@ -11,13 +11,13 @@ class Action { virtual std::string calculate() = 0; }; -class quadraticAction : public Action { +class quadraticAction final : public Action { std::string calculate() override { return "calculated by quadratic action"; }; }; -class barrettGlaserAction : public Action { +class barrettGlaserAction final : public Action { std::string calculate() override { return "calculated by Barrett-Glaser action"; } @@ -33,8 +33,8 @@ class ActionManager { m_action = std::move(new_action); } - void Print() { - auto action_val = m_action->calculate(); + void Print() const { + const auto action_val = m_action->calculate(); std::cout << action_val << std::endl; }