diff --git a/.clang-format b/.clang-format
new file mode 100755
index 0000000..eb26853
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,102 @@
+---
+AccessModifierOffset: 0
+AlignAfterOpenBracket: Align
+AlignArrayOfStructures: None
+AlignConsecutiveAssignments:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: true
+ AlignFunctionDeclarations: true
+ AlignFunctionPointers: true
+ PadOperators: true
+AlignConsecutiveBitFields:
+ Enabled: false
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: false
+ AlignFunctionDeclarations: false
+ AlignFunctionPointers: false
+ PadOperators: false
+AlignConsecutiveDeclarations:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: false
+ AlignFunctionDeclarations: true
+ AlignFunctionPointers: false
+ PadOperators: false
+AlignConsecutiveMacros:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: false
+ AlignFunctionDeclarations: false
+ AlignFunctionPointers: false
+ PadOperators: false
+AlignConsecutiveShortCaseStatements:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCaseArrows: false
+ AlignCaseColons: false
+AlignConsecutiveTableGenCondOperatorColons:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: false
+ AlignFunctionDeclarations: false
+ AlignFunctionPointers: false
+ PadOperators: false
+AlignConsecutiveTableGenDefinitionColons:
+ Enabled: true
+ AcrossEmptyLines: false
+ AcrossComments: false
+ AlignCompound: true
+ AlignFunctionDeclarations: false
+ AlignFunctionPointers: false
+ PadOperators: false
+AlignEscapedNewlines: Left
+AlignOperands: Align
+AlignTrailingComments:
+ Kind: Always
+ OverEmptyLines: 0
+AllowAllArgumentsOnNextLine: true
+AllowBreakBeforeNoexceptSpecifier: Never
+AllowShortCaseExpressionOnASingleLine: true
+AllowShortCompoundRequirementOnASingleLine: true
+AllowAllParametersOfDeclarationOnNextLine: false
+AllowShortBlocksOnASingleLine: true
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortEnumsOnASingleLine: true
+AllowShortFunctionsOnASingleLine: true
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLambdasOnASingleLine: Empty
+AllowShortLoopsOnASingleLine: false
+BinPackArguments: false
+BinPackParameters: false
+BraceWrapping:
+ AfterCaseLabel: false
+ AfterClass: false
+ AfterControlStatement: false
+ AfterEnum: false
+ AfterFunction: false
+ AfterNamespace: false
+ AfterObjCDeclaration: false
+ AfterStruct: false
+ AfterUnion: false
+ BeforeCatch: false
+ BeforeElse: false
+ IndentBraces: false
+BreakBeforeBraces: Custom
+ColumnLimit: 120
+IncludeBlocks: Preserve
+IndentWidth: 4
+InsertNewlineAtEOF: true
+KeepEmptyLinesAtTheStartOfBlocks: true
+QualifierAlignment: Left
+CommentPragmas: '^.*A2Lfactory:'
+---
+Language: Cpp
+Standard: c++20
+...
diff --git a/.clang-tidy b/.clang-tidy
new file mode 100644
index 0000000..19d823e
--- /dev/null
+++ b/.clang-tidy
@@ -0,0 +1,15 @@
+Checks: >
+ -*,
+ clang-analyzer-*,
+ cert-*,
+ cppcoreguidelines-*,
+ bugprone-*,
+ google-*,
+ misc-*,
+ performance-*,
+ readability-*,
+ modernize-*,
+ -modernize-use-trailing-return-type,
+ -readability-identifier-length,
+ -misc-non-private-member-variables-in-classes
+
diff --git a/.gitignore b/.gitignore
index 1e17110..c3dc3b5 100755
--- a/.gitignore
+++ b/.gitignore
@@ -1,23 +1,28 @@
.DS_Store
.directory
+toolchain.cmake
+actions-runner/*
+actions-runner/_work
actions-runner
-CMakeLists.txt
assets/flippy.psd
cmake-build*
.idea
.run
-
+compile_commands.json
+build-*
+build/
+cmakerc
+logs/
+.clangd
demo_out
flippy/flippy.hpp-e
discarded
tests/test.info
tests/test_out
-logs
tests/res
banner.old
doxygen_warnings.log
make_doxygen_html.sh
make_single_header_flippy.sh
-benchmarks
-code_overview
/actions-runner/_work
+/demo/planar_membrane_sheet_fluctuation_MC/
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100755
index 0000000..c8fa672
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,51 @@
+cmake_minimum_required(VERSION 3.16)
+project(flippy)
+file(READ "VERSION.json" ver)
+string(JSON ver GET ${ver} VERSION)
+string(REGEX MATCHALL "[0-9]+.[0-9]+.[0-9]+" version ${ver})
+string(REGEX MATCHALL "-.+" pre_release_flag ${ver})
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+if("${pre_release_flag}" STREQUAL "")
+ message("no prerelease flags to set")
+ project(flippy VERSION ${version})
+else()
+ project(
+ flippy
+ VERSION ${version}
+ DESCRIPTION ${pre_release_flag})
+ message("prerelease flag: ${pre_release_flag}")
+endif()
+
+if(${CMAKE_CXX_COMPILER_ID} STREQUAL MSVC)
+ message("building for MSVC")
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O2")
+else()
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
+endif()
+
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -g -fno-omit-frame-pointer")
+
+set(FLIPPY_COMPILE_OPTIONS
+ -Wall
+ -Wextra
+ -Wpedantic
+ -Wshadow
+ -Wconversion
+ -Wswitch-enum
+ -Werror
+ -ffast-math
+ # gToDo: remove this after fixing nlohmann dependency
+ -Wno-nan-infinity-disabled
+ # gToDo: remove this after fixing catch2 dependency
+ -Wno-c2y-extensions)
+
+set(CMAKE_CXX_STANDARD 26)
+set(CXX_STANDARD_REQUIRED ON)
+
+enable_testing()
+
+add_subdirectory(third_party)
+add_subdirectory(tests)
+add_subdirectory(demo)
+add_subdirectory(benchmarks)
diff --git a/README.md b/README.md
index 3a26268..ba82af1 100755
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@
[](mailto:flippy@mailbox.org)
# *flippy*
-
+
c++ package for dynamically triangulated membrane simulations.
@@ -19,11 +19,11 @@ c++ package for dynamically triangulated membrane simulations.
| | |
# Support
-*flippy* is still in active development but has the first stable release. And a full public [API documentation](https://flippy-software-package.github.io/flippy). I am actively working on flippy and welcome pull requests and bug reports.
+*flippy* is still in active development but has the first stable release. And a full public [API documentation](https://flippy-software-package.github.io/flippy). I am actively working on flippy and welcome pull requests and bug reports.
I will also consider feature requests, but in that case, please provide a use case.
### for questions about general usage
please use the support email [](mailto:flippy@mailbox.org).
-### for bugfixes
+### for bugfixes
please create an [issue](https://github.com/flippy-software-package/flippy/issues).
### for feature requests
again the [issues](https://github.com/flippy-software-package/flippy/issues) page can be used, but be aware that new features will be slow to come.
@@ -32,10 +32,8 @@ again the [issues](https://github.com/flippy-software-package/flippy/issues) pag
*flippy* is a header-only library, so all you need to do is to download the `flippy` sub-folder and copy it into your project.
-Or, if you prefer using a single header file, you can download the [flippy.hpp](https://raw.githubusercontent.com/flippy-software-package/flippy/master/single_header_flippy/flippy.hpp) header from the `single_header_flippy` folder.
-
# Documentation
-*flippy*'s public API is fully [documented](https://flippy-software-package.github.io/flippy). However, if you have never used *flippy* before, you might want to start with the demo projects. These are located in sub-folders of the [`demo`](https://github.com/flippy-software-package/flippy/tree/master/demo) folder. Each sub-folder contains a readme like this one, explaining how to set up a project and what to expect.
+*flippy*'s public API is fully [documented](https://flippy-software-package.github.io/flippy). However, if you have never used *flippy* before, you might want to start with the demo projects. These are located in sub-folders of the [`demo`](https://github.com/flippy-software-package/flippy/tree/master/demo) folder. Each sub-folder contains a readme like this one, explaining how to set up a project and what to expect.
Automatically generated code [documentation](https://github.com/flippy-software-package/flippy/wiki/Documentation) over on the [wiki](https://github.com/flippy-software-package/flippy/wiki).
# Citing flippy
@@ -64,86 +62,125 @@ This code can be found in the sub-folder `demo/biconcave_shapes_MC` together wit
```cpp
//demo/biconcave_shapes_MC/main.cpp
-#include // needed for random displacement generation
-#include // need for std::vector
-#include // needed for std::cout
#include "flippy.hpp"
+#include // needed for std::cout
+#include // needed for random displacement generation and node index shuffling
+#include // need for std::vector
+
+fp::Real sphere_vol(fp::Real R) { return fp::Real(4. / 3.) * fp::PI * R * R * R; }
+fp::Real sphere_area(fp::Real R) { return fp::Real(4.) * fp::PI * R * R; }
+
+struct EnergyParameters {
+ fp::Real kappa, K_V, K_A, V_t, A_t;
+};
+
+// This is the energy function that is used by flippy's built-in updater to decide if a move was energetically favorable
+// or not
+fp::Real surface_energy(const fp::Node &node,
+ const fp::Triangulation &trg,
+ const EnergyParameters &p,
+ const std::vector &changed_neighbourhood) {
+ fp::Real V = trg.global_geometry().volume;
+ fp::Real A = trg.global_geometry().area;
+ fp::Real dV = V - p.V_t;
+ fp::Real dA = A - p.A_t;
+ fp::Real E_bending = fp::node_unit_bending_energy(node);
+ for (fp::Index changed_node_id : changed_neighbourhood) {
+ E_bending += fp::node_unit_bending_energy(node);
+ }
-double sphere_vol(double R){return 4./3. * M_PI *R*R*R;}
-double sphere_area(double R){return 4. * M_PI *R*R;}
-
-struct EnergyParameters{double kappa, K_V, K_A, V_t, A_t;};
-
-// This is the energy function that is used by flippy's built-in updater to decide if a move was energetically favorable or not
-double surface_energy([[maybe_unused]]fp::Node const& node,
- fp::Triangulation const& trg,
- EnergyParameters const& prms){
- double V = trg.global_geometry().volume;
- double A = trg.global_geometry().area;
- double dV = V-prms.V_t;
- double dA = A-prms.A_t;
- double energy = prms.kappa*trg.global_geometry().unit_bending_energy +
- prms.K_V*dV*dV/prms.V_t + prms.K_A*dA*dA/prms.A_t;
+ fp::Real energy = (p.kappa * E_bending) + (p.K_V * dV * dV / p.V_t) + (p.K_A * dA * dA / p.A_t);
return energy;
}
-int main(){
- int n_triang = 7; // triangulation iteration number of nodes N_node=12+30*n+20*n*(n-1)/2 where n is the same as n_trng
- double l_min = 2;
- double R = l_min/(2*sin(asin(1./(2*sin(2.*M_PI/5.)))/(n_triang+1.))); // estimate of a typical bond length in the initial triangulation and then create a sphere such that the initial bond length is close to minimal. This formula is derived from the equidistant sub-triangulation of an icosahedron, where geodesic distances are used as a distance measure.
- double l_max = 2.*l_min; // if you make l_max closer to l_min bond_flip acceptance rate will go down
- double r_Verlet = 2*l_max;
- EnergyParameters prms{.kappa=10 /*kBT*/,
- .K_V=100 /*kBT/area*/, .K_A=1000 /*kBT/volume*/,
- .V_t=0.6*sphere_vol(R), .A_t=sphere_area(R)};
- double linear_displ=l_min/8.; // side length of a voxel from which the displacement of the node is drawn
- int max_mc_steps=1e5; // max number of iteration steps (depending on the strength of your CPU, this should take anywhere from a couple of seconds to a couple of minutes
+#include "../benchmarks/external/code_utils.hpp" //ToDo Remove before prod
+
+int main() {
+ cutils::Timer t;
+ int n_triang =
+ 7; // triangulation iteration number of nodes N_node=12+30*n+20*n*(n-1)/2 where n is the same as n_trng
+ fp::Real l_min = 2;
+ fp::Real R =
+ l_min /
+ (2 * sin(asin(1. / (2 * sin(2. * M_PI / 5.))) / (n_triang + 1.))); // estimate of a typical bond length in the
+ // initial triangulation and then create a sphere such that the initial bond length is close to minimal. This
+ // formula is derived from the equidistant sub-triangulation of an icosahedron, where geodesic distances are used as
+ // a distance measure.
+ fp::Real l_max = 1.8 * l_min; // if you make l_max closer to l_min bond_flip acceptance rate will go down
+ fp::Real r_Verlet = 2 * l_max;
+ fp::Real red_vol = 0.6;
+ EnergyParameters prms{.kappa = 10 /*kBT*/,
+ .K_V = 100 /*kBT/area*/,
+ .K_A = 1000 /*kBT/volume*/,
+ .V_t = red_vol * sphere_vol(R),
+ .A_t = sphere_area(R)};
+ fp::Real linear_displ = l_min / 18.; // side length of a voxel from which the displacement of the node is drawn
+ int max_mc_steps = 2e5; // max number of iteration steps (depending on the strength of your CPU, this should take
+ // anywhere from a couple of seconds to a couple of minutes
std::random_device random_number_generator_seed;
- std::mt19937 rng(random_number_generator_seed()); // create a random number generator and seed it with the current time
+ std::mt19937 rng(
+ random_number_generator_seed()); // create a random number generator and seed it with the current time
// All the flippy magic is happening on the following two lines
- fp::Triangulation guv(n_triang, R, r_Verlet);
- fp::MonteCarloUpdater mc_updater(guv, prms, surface_energy, rng, l_min, l_max);
+ fp::Triangulation guv(n_triang, R, r_Verlet);
+ fp::MonteCarloUpdater mc_updater(guv, prms, surface_energy, rng, l_min, l_max);
- fp::vec3 displ{}; // declaring a 3d vector (using flippy's built in vec3 type) for later use as a random direction vector
- std::uniform_real_distribution displ_distr(-linear_displ, linear_displ); //define a distribution from which the small displacements in x y and z directions will be drawn
+ fp::vec3 displ{}; // declaring a 3d vector (using flippy's built in vec3 type) for later use as a random
+ // direction vector
+ std::uniform_real_distribution displ_distr(-linear_displ, linear_displ); // define a distribution from
+ // which the small displacements in x y and z directions will be drawn
- guv.scale_node_coordinates(1, 1, 0.8); // squish the sphere in the z-direction to break the initial symmetry. This speeds up the convergence to a biconcave shape greatly.
+ guv.scale_node_coordinates(1, 1, 0.85); // squish the sphere in the
+ // z-direction to break the initial symmetry. This speeds up the convergence to a biconcave shape greatly.
fp::Json data_init = guv.make_egg_data();
- fp::json_dump("test_run_init", data_init); // ATTENTION!!! this file will be saved in the same folder as the executable
+ fp::json_dump("test_run_init",
+ data_init); // ATTENTION!!! this file will be saved in the same folder as the executable
std::vector shuffled_ids;
shuffled_ids.reserve(guv.size());
- for(auto const& node: guv.nodes()){ shuffled_ids.push_back(node.id);} //create a vector that contains all node ids. We can shuffle this vector in each MC step to iterate randomly through the nodes
+ for (const auto &node : guv.nodes()) {
+ shuffled_ids.push_back(node.id);
+ } // create a vector that contains all node ids. We can shuffle this vector in each MC step to iterate randomly
+ // through the nodes
- for(int mc_step=0; mc_step
+
+
diff --git a/assets/flippy_icon.png b/assets/flippy_icon.png
old mode 100644
new mode 100755
diff --git a/assets/move and flipp updates.png b/assets/move and flipp updates.png
old mode 100644
new mode 100755
diff --git a/assets/structure_of_flippy.png b/assets/structure_of_flippy.png
old mode 100644
new mode 100755
diff --git a/assets/triangulation.png b/assets/triangulation.png
old mode 100644
new mode 100755
diff --git a/assets/triangulation.svg b/assets/triangulation.svg
old mode 100644
new mode 100755
diff --git a/assets/triangulation_four_panels.svg b/assets/triangulation_four_panels.svg
old mode 100644
new mode 100755
diff --git a/banner b/banner
index 2c2649e..5fa220b 100755
--- a/banner
+++ b/banner
@@ -9,7 +9,7 @@
* 888 888 888 888 888 888 888 888 888 dynamically triangulated
* 888 888 888 888 d88P 888 d88P Y88b 888 surfaces
* 888 888 888 88888P" 88888P" "Y88888
- * 888 888 888 version 1.0.1
+ * 888 888 888 version 1.1.0
* 888 888 Y8b d88P
* 888 888 "Y88P"
*
diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
new file mode 100755
index 0000000..c6456d3
--- /dev/null
+++ b/benchmarks/CMakeLists.txt
@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 3.16)
+
+# include_directories(micro_benchmarks)
+include_directories(external/google_benchmark)
+include_directories(external)
+include_directories(benchmark_logger)
+include_directories(../flippy)
+include_directories(../flippy/utilities)
+include_directories(.)
+
+# add_subdirectory(micro_benchmarks/vec3_bench)
+add_subdirectory(rbc_bench)
diff --git a/benchmarks/benchmark-main.zip b/benchmarks/benchmark-main.zip
new file mode 100755
index 0000000..b7d5ff3
Binary files /dev/null and b/benchmarks/benchmark-main.zip differ
diff --git a/benchmarks/benchmark_logger/FlippyBenchmarkLogger.h b/benchmarks/benchmark_logger/FlippyBenchmarkLogger.h
new file mode 100644
index 0000000..901fdba
--- /dev/null
+++ b/benchmarks/benchmark_logger/FlippyBenchmarkLogger.h
@@ -0,0 +1,108 @@
+#ifndef NODES_HPP_FLIPPYBENCHMARKLOGGER_H
+#define NODES_HPP_FLIPPYBENCHMARKLOGGER_H
+#include
+#include
+#include
+
+template
+
+concept Updater = requires(U u,
+ fp::Triangulation triangulation,
+ fp::Node node,
+ fp::Index nn_id,
+ fp::vec3 displacement,
+ unsigned long long z,
+ std::ostream &os) {
+ // { U(triangulation, node) } -> std::same_as;
+
+ { u.move_MC_updater(node, displacement) } -> std::same_as>;
+ { u.flip_MC_updater(node, nn_id) } -> std::same_as;
+
+ { u.move_attempt_count() } -> std::same_as;
+ { u.bond_length_move_rejection_count() } -> std::same_as;
+ { u.move_back_count() } -> std::same_as;
+
+ { u.flip_attempt_count() } -> std::same_as;
+ { u.bond_length_flip_rejection_count() } -> std::same_as;
+ { u.flip_back_count() } -> std::same_as;
+};
+
+template class FlippyBenchmarkLogger {
+ const Updater &mcu_;
+ const fp::Triangulation &triangulation_;
+ cutils::Timer timer{};
+
+ FlippyBenchmarkLogger(const Updater &mcu, const fp::Triangulation &triangulation)
+ : mcu_(mcu), triangulation_(triangulation) {};
+
+ public:
+ static FlippyBenchmarkLogger start_benchmark(const Updater &mcu, const fp::Triangulation &triangulation) {
+ FlippyBenchmarkLogger logger(mcu, triangulation);
+
+ PRINT("each of ", triangulation.size(), " beads will be attempted to be moved: ", "XXX", "(predicted) times");
+ PRINT("leading to a total of", "XXX", "(predicted) MC moves");
+
+ return logger;
+ }
+
+ std::string log(const std::string &log_dir, const fp::Json &config_data) {
+ cutils::HumanReadableTime elapsed_time = timer.stop();
+ // MonteCarloUpdater counts the number of accepted and rejected moves, distinguishing whether a rejection
+ // occurred because of the energy or the bond length constraint. We can use this to print simple statistics
+ // here. For example, this will help us decide if our displacement size is too large.
+ auto attempt = mcu_.move_attempt_count();
+ auto put_back = mcu_.move_back_count();
+ auto bond_length_legal_move = attempt - mcu_.bond_length_move_rejection_count();
+
+ auto attempt_flip = mcu_.flip_attempt_count();
+ auto flip_back = mcu_.flip_back_count();
+
+ auto topologically_successful_flip = mcu_.flip_attempt_count() - flip_back;
+ PRINT("percentage of failed moves: ",
+ static_cast(mcu_.move_back_count() + mcu_.bond_length_move_rejection_count()) /
+ static_cast(mcu_.move_attempt_count()));
+ PRINT("percentage of failed flips: ",
+ static_cast(mcu_.flip_back_count() + mcu_.bond_length_flip_rejection_count()) /
+ static_cast(mcu_.flip_attempt_count()));
+
+ long double time_in_seconds = ((long double)elapsed_time.diff_ns) * 1.e-9l;
+ long double moves_ps = static_cast(attempt) / static_cast(time_in_seconds);
+ long double e_moves_ps =
+ static_cast(bond_length_legal_move) / static_cast(time_in_seconds);
+ long double flips_ps = static_cast(attempt_flip) / static_cast(time_in_seconds);
+ long double e_flips_ps =
+ static_cast(topologically_successful_flip) / static_cast(time_in_seconds);
+
+ auto now = std::chrono::system_clock::now();
+ std::time_t now_time = std::chrono::system_clock::to_time_t(now);
+ std::tm *time = std::localtime(&now_time);
+ std::stringstream date, time_stamp;
+ date << 1900 + time->tm_year << "-" << time->tm_mon + 1 << "-" << time->tm_mday;
+ time_stamp << time->tm_hour << "-" << time->tm_min << "-" << time->tm_sec;
+
+ std::string log_sub_dir = log_dir + date.str() + "/";
+ std::filesystem::create_directories(log_sub_dir);
+
+ std::string log_name = log_sub_dir + time_stamp.str() + ".yml";
+
+ std::ofstream log_file(log_name);
+ log_file << "BENCHMARK VERSION: " << config_data["benchmark_version"].get() << '\n';
+ log_file << "config: " << config_data << '\n';
+ log_file << "all_attempted_bead_moves: " << attempt << '\n';
+ log_file << "expensive_moves: " << bond_length_legal_move << '\n';
+ log_file << "rejected_moves_stochastic: " << put_back << '\n';
+ log_file << "all_attempted_flips: " << attempt_flip << '\n';
+ log_file << "expensive_flips: " << topologically_successful_flip << '\n';
+ log_file << "rejected_flip_stochastic: " << flip_back << '\n';
+ log_file << "total_simulation_time: " << elapsed_time.diff << " " << elapsed_time.unit << " ("
+ << elapsed_time.diff_fine << " " << elapsed_time.unit_fine << ")" << '\n';
+ log_file << "mps: " << moves_ps << '\n';
+ log_file << "emps: " << e_moves_ps << '\n';
+ log_file << "fps: " << flips_ps << '\n';
+ log_file << "efps: " << e_flips_ps << '\n';
+ log_file.close();
+ return log_name;
+ }
+};
+
+#endif // NODES_HPP_FLIPPYBENCHMARKLOGGER_H
diff --git a/benchmarks/experimental/Node_Static_NN_ARRAYS.hpp b/benchmarks/experimental/Node_Static_NN_ARRAYS.hpp
new file mode 100644
index 0000000..84bf4c3
--- /dev/null
+++ b/benchmarks/experimental/Node_Static_NN_ARRAYS.hpp
@@ -0,0 +1,250 @@
+#ifndef FLIPPY_NODES_EXPERIMENTAL_HPP
+#define FLIPPY_NODES_EXPERIMENTAL_HPP
+/**
+ * @file
+ * @brief This file contains the fp::Node and fp::Nodes classes, data structures that represent a single node of the
+ * triangulation and the collection of all nodes of the triangulation, respectively.
+ */
+#include "external/json.hpp"
+#include "vec3.hpp"
+#include
+#include
+
+namespace fp::experimental {
+
+constexpr unsigned int MAX_ALLOWED_NEIGHBORS = 16;
+
+template class NeighborsArray : public std::array {
+ public:
+ using std::array::array;
+ size_t size() { return size_; }
+
+ void emplace(size_t idx, T &&val) {
+ for (size_t i = size_; i > idx; --i) {
+ this->data()[i] = this->data()[i - 1];
+ }
+ this->data()[idx] = val;
+ size_++;
+ }
+ // void pop(size_t idx){
+ // for(size_t i = size_ - 1; )
+ // move_left(idx, this);
+ // }
+
+ private:
+ // void move_left(size_t start, )
+ size_t size_{};
+};
+
+using Json = nlohmann::json;
+//! A data structure containing all geometric and topological information associated with a node.
+/**
+ * This is a DUMB DATA STRUCTURE, meaning that it is not responsible for the coherence of the data it contains.
+ * For performance reasons, methods associated with Node struct will never check if the Node::curvature is the norm of
+ * the Node::curvature_vector or if the Node::nn_ids and Node::nn_distances are in the correct order. It is the
+ * responsibility of higher-order structures like Nodes and Triangulation to check that correct data is stored and
+ * updated correctly. However, it does check the data for consistency. It will match the length of Node::nn_ids and
+ * Node::nn_distances and pop and add both of them together.
+ * @tparam Real @RealStub
+ * @tparam Index @IndexStub
+ */
+template struct Node {
+ //! @NodeIDStub
+ Index id;
+ //! Voronoi area associated with the node.
+ /**
+ * The Voronoi area is the sum of (mixed) Voronoi areas inside the triangles, incident to the node.
+ * Definition follows [Gueguen et al. 2017](https://doi.org/10.1039/C7SM01272A).
+ * \f[
+ * A_{i} = \sum_{j} A'_{ij}.
+ * \f]
+ *
+ * @see Triangulation::mixed_area
+ * See Figure tr1. C in Triangulation.
+ * @see Node::curvature_vec Triangulation::update_bulk_node_geometry(Index)
+ */
+ Real area;
+ //! If the node is part of a closed surface triangulation, then the `volume` contains the volume of the tetrahedron
+ //! connected to each voronoi cell sub-triangle and the center of the lab coordinate system as defined in [Gueguen
+ //! et al. 2017](https://doi.org/10.1039/C7SM01272A).
+ /** * This means that the volume of an individual node does not have a proper physical interpretation.
+ * Only the sum of all node volumes, which is given by the triangulation
+ * is interpretable as a physical volume of an object.
+ * The definition follows [Gueguen et al. 2017](https://doi.org/10.1039/C7SM01272A).
+ * \f[
+ * V_{ij} = A_{ij} \vec{x}_{i}\cdot \frac{\vec{n}_{ij,j+1}}{\| \vec{n}_{ij,j+1} \|}.
+ * \f]
+ * See Figure tr1. D in Triangulation.
+ * @see Node::curvature_vec Triangulation::update_bulk_node_geometry(Index)
+ */
+ Real volume;
+ //! `unit_bending_energy` corresponds to the [Helfrich bending
+ //! energy](https://en.wikipedia.org/wiki/Elasticity_of_cell_membranes) with bending rigidity 1 and gaussian bending
+ //! stiffness 0.
+ /**
+ * \f[
+ * \mathrm{unit\_bending\_energy} = \frac{1}{2} A_{\mathrm{node}} (2 H_{node})^2
+ * \f]
+ * where \f$ H_{node} \f$ is the mean curvature of the node given by:
+ * \f[
+ * H_{node}^2 = \frac{\vec{K}_{node}}{2A_{node}} \cdot \frac{\vec{K}_{node}}{2A_{node}}
+ * \f],
+ * with \f$ \vec{K} \f$ denoting the Node::curvature_vector.
+ * @see Node::curvature_vec Triangulation::update_bulk_node_geometry(Index)
+ */
+ Real unit_bending_energy;
+ //! Position of the node in the lab frame.
+ vec3 pos;
+ //! Curvature vector of the node.
+ /**
+ * The definition of the curvature vector follows [Meyer et al. 2003](https://doi.org/10.1007/978-3-662-05105-4_2).
+ * \f[
+ * \vec{K}_i = \frac{1}{2A_i}\sum_{j(i)} \left( \cot\left(\alpha_{ij}^{j+1}\right) +
+ * \cot\left(\alpha_{ij}^{j-1}\right) \right)\vec{\ell}_{ij}
+ * \f]
+ * See Figure tr1. B in Triangulation.
+ * @see Node::curvature_vec Triangulation::update_bulk_node_geometry(Index)
+ */
+ vec3 curvature_vec;
+ //! A vector containing the global ids of the current node's next neighbors.
+ /**
+ * `nn_ids` contains the ids of nodes that are connected to this node in the triangulation.
+ * The next neighbors that are also mutual neighbors in the triangulation are stored sequentially in the vector.
+ * The last and the first elements are also neighbors, i.e., the nn_ids vector wraps around.
+ * During the calculation, this is facilitated through the use of @ref fp::Neighbors.
+ * @note The order of the next neighbors matters for the proper function of fp::Triangulation but is not guaranteed
+ * by this data structure. See Figure tr1. A, in Triangulation.
+ */
+ NeighborsArray nn_ids;
+ //! Distance vectors pointing from the node to its next neighbors.
+ NeighborsArray> nn_distances;
+ //! The Verlet list contains the ids of nodes that are close to this node.
+ std::vector verlet_list;
+
+ // unit-tested
+ //! Find and deletes the element with the id `to_pop_nn_id` in the `nn_id` vector.
+ void pop_nn(Index to_pop_nn_id) {
+ /**
+ * @param to_pop_nn_id @NNIDStub This id is supposed to be removed from the next neighbor id vector.
+ * @see Node::nn_ids
+ * @note this will lead to resizing of the vector, which can be expensive!
+ * @warning If the provided next neighbor id is not part of the Node::nn_ids, this function will fail silently.
+ * It will not delete anything and won't throw any errors or warnings;
+ */
+ auto pop_pos = find_nns_loc_pointer(to_pop_nn_id);
+ auto dist = pop_pos - nn_ids.begin();
+
+ if (pop_pos != nn_ids.end()) {
+ // I checked that this would work on example code on cppreference https://godbolt.org/z/6qf8c9nTz
+ nn_ids.erase(pop_pos);
+ nn_distances.erase(nn_distances.begin() + dist);
+ }
+ }
+
+ auto find_nns_loc_pointer(Index nn_id) {
+ /**
+ * @brief Given the global id of the next neighbor, this function can be used to locate it in the Node::nn_ids
+ * vector.
+ *
+ * This function is just a convenient wrapper around the
+ * [std::find](https://en.cppreference.com/w/cpp/algorithm/find) function.
+ * ```
+ * std::find(nn_ids.begin(), nn_ids.end(), to_pop_nn_id);
+ * ```
+ * @param nn_id @NNIDStub
+ * @return if `nn_id` is contained in Node::nn_ids then the pointer to the position of that id in the `nn_ids`
+ * vector will be returned. Otherwise `nn_ids.end()`.
+ * @warning This function is not responsible for graceful handling of `nn_id`'s that are not found in the
+ * Node::nn_ids vector. If the `nn_id` is not contained in Node::nn_ids then the `nn_ids.end()` iterator will be
+ * returned. It is up to the user to perform the necessary checks to avoid undefined behavior that might result
+ * from trying to delete uninitiated memory.
+ */
+ return std::find(nn_ids.begin(), nn_ids.end(), nn_id);
+ }
+
+ // unit-tested
+ void emplace_nn_id(Index to_emplace_nn_id, const vec3 &to_emplace_nn_pos, Index loc_idx) {
+ /**
+ * @brief This function can be used to add new next neighbors to this node.
+ *
+ * This function constructs `to_emplace_nn_id` right before `to_emplace_pos`,
+ * i.e. if to_emplace_nn_id is 3, to_emplace_nn_id will be constructed right before the
+ * 3rd element and will become the new 3rd element.
+ * @param to_emplace_nn_id @NNIDStub This id is supposed to be added to the Node::nn_ids vector of this node.
+ * @param to_emplace_nn_pos const reference to the 3 dimensional position vector (type vec3) containing
+ * the position of the new next neighbour. This input is used to calculate the correct distance between this
+ * node and the new next neighbor, which will then be added to the Node::nn_distances vector.
+ * @param loc_idx @LocNNIndexStub
+ * @note This function causes the resizing of two vectors, which can be costly.
+ * @warning Making next neighbors is a symmetric operation. I.e., if node one becomes the next neighbor of node
+ * two, node two also has to become the next neighbor of node one. However, this function is not responsible for
+ * this relationship. It only adds a new next neighbor to this node, and the higher-order structures, like
+ * Triangulation, are responsible for guaranteeing the reciprocal relationship.
+ * @see Triangulation::emplace_before(Index, Index, Index)
+ */
+ if (loc_idx < nn_ids.size()) {
+ nn_ids.emplace(loc_idx, to_emplace_nn_id);
+ nn_distances.emplace(loc_idx, to_emplace_nn_pos - pos);
+ }
+ }
+
+ // unit-tested
+ //! This function can provide the stored distance vector to the next neighbor.
+ const vec3 &get_distance_vector_to(Index nn_id) const {
+ /**
+ * @param nn_id @NNIDStub.
+ * @return returns the distance currently stored in the Node::nn_distances vector for the requested next
+ * neighbor. If the provided `nn_id` can not be found in the Node::nn_ids vector, then the function writes an
+ * error message to standard error output and terminates the program with exit code 12.
+ * @note @TerminationNoteStub
+ */
+ auto id_pos = std::find(nn_ids.begin(), nn_ids.end(), nn_id);
+ if (id_pos != nn_ids.end()) {
+ return nn_distances[static_cast(id_pos - nn_ids.begin())];
+ } else {
+ std::cerr << "nn_id:" << nn_id
+ << " provided to `get_distance_vector_to` is not a next neighbour of the node " << id;
+ exit(12);
+ }
+ }
+
+ // defaulted operators are not explicitly unit-tested
+ /**
+ * @brief Default equality operator.
+ *
+ * @param other_node constant reference to the other Node.
+ * @return True if both nodes are equal.
+ */
+ bool operator==(const Node &other_node) const = default;
+
+ /**
+ * @brief Streaming operator that can print formatted output to standard out with all Node data fields.
+ *
+ * @param os This is intended to be std::cout or any other std::ofstream reference.
+ * @param node The streamed node.
+ * @return Updated stream.
+ */
+ friend std::ostream &operator<<(std::ostream &os, const Node &node) {
+
+ os << "node: " << node.id << '\n'
+ << "area: " << node.area << '\n'
+ << "volume: " << node.volume << '\n'
+ << "unit_bending_energy: " << node.unit_bending_energy << '\n'
+ << "curvature_vec: " << node.curvature_vec << '\n'
+ << "pos: " << node.pos << '\n'
+ << "nn_ids: ";
+ for (const auto &nn_id : node.nn_ids) {
+ os << nn_id << ' ';
+ }
+ os << '\n' << "nn_distances: ";
+ for (const auto &nn_dist : node.nn_distances) {
+ os << nn_dist << '\n';
+ }
+ os << '\n';
+
+ return os;
+ }
+};
+
+} // namespace fp::experimental
+#endif // FLIPPY_NODES_HPP
diff --git a/benchmarks/experimental/vec3_valarray.h b/benchmarks/experimental/vec3_valarray.h
new file mode 100644
index 0000000..0f67dd2
--- /dev/null
+++ b/benchmarks/experimental/vec3_valarray.h
@@ -0,0 +1,335 @@
+#ifndef FLIPPY_VEC3_EXPERIMENTAL_HPP
+#define FLIPPY_VEC3_EXPERIMENTAL_HPP
+/**
+ * @file
+ * @brief Header file containing the definition and implementation a 3 dimensional vector class, with useful
+ * mathematical operations like cross and dot products as member methods.
+ */
+
+#include "custom_concepts.hpp"
+#include
+#include
+#include
+#include
+
+namespace fp::experimental {
+
+/**
+ * \brief Internal implementation of a 3D vector.
+ *
+ * !!! vec3 does not throw !!! This means that if you ask vec3 to divide a vector by 0 or more realistically if you
+ * normalize a zero length vector vec3 will not check for the division by zero and will return a nan result!
+ * Since vec3 is used everywhere in flippy, including in very expensive calculations, I decided to omit the security
+ * check for the sake of speed.
+ *
+ * To keep the external dependencies low, flippy implements it's own 3D vector class with basic functionality like dot
+ * product and cross product
+ *
+ * Example:
+ * ```c++
+ * fp::vec3 v1{1,0,0};
+ * fp::vec3 v2{0,0,1};
+ *
+ * assert(v1.dot(v2)==0);
+ * assert(v1.cross(v2).norm()==1);
+ * assert(((v1-v2)==fp::vec3{1.,0.,-1.}));
+ * ```
+ *
+ * @tparam Real @RealStub
+ */
+
+template class vec3 {
+ std::valarray data_;
+
+ public:
+ vec3() : data_(std::valarray(3)) {}
+ vec3(Real x, Real y, Real z) : data_(std::valarray{x, y, z}) {}
+ vec3(Real &&x, Real &&y, Real &&z) : data_(std::valarray{x, y, z}) {}
+ // vec3(Real const&x, Real const&y, Real const&z): data_(std::valarray{x, y, z}){}
+
+ Real x() { return data_[0]; } //!< The x component of the vector.
+ Real y() { return data_[1]; } //!< The y component of the vector.
+ Real z() { return data_[2]; } //!< The z component of the vector.
+
+ //! In place addition method.
+ /**
+ * Example:
+ * ```c++
+ * fp::vec3 v1{1,0,0};
+ * fp::vec3 v2{0,0,1};
+ * v1.add(v2); // v1 will contain {1, 0, 1}
+ * ```
+ * @param v add this vector elementwise to the vector that is calling the *add* method.
+ */
+ void add(const vec3 &v) { data_ += v.data_; }
+
+ //! In place subtraction method.
+ /**
+ * Example:
+ * ```c++
+ * fp::vec3 v1{2,0,0};
+ * fp::vec3 v2{1,0,1};
+ * v1.subtract(v2); // v1 will contain {1, 0, -1}
+ * ```
+ * @param v subtract this vector elementwise from the vector that is calling the *subtract* method.
+ */
+ void subtract(const vec3 &v) { data_ -= v.data_; }
+ //! Scale the vector by a real number s.
+ /**
+ * This function scales the vector in-place by the provided number `s`.
+ * @param s multiplicative prefactor.
+ */
+ void scale(Real s) { data_ *= s; }
+
+ //! Calculate dot product with another vector.
+ /**
+ * Example:
+ * @code{c++}
+ * fp::vec3 v1{1,0,0};
+ * fp::vec3 v2{2,0,1};
+ * double res = v1.dot(v2); // res will contain 2*1 + 0*0 + 0*1=2
+ * @endcode
+ * @param v the other vec3 vector
+ * @return result of the dot product between the original vector and `v`.
+ */
+ Real dot(const vec3 &v) const { return (data_ * v.data_).sum(); }
+
+ //! Always returns 3.
+ /**
+ * This function always returns 3 since vec3 can only have three elements.
+ * It was implemented for completeness, to make it more easy for vec3 to be used as a drop-in replacement for other
+ * vector types.
+ * @return Size (number of elements) of vec3.
+ */
+ [[nodiscard]] constexpr std::size_t size() const { return data_.size(); }
+
+ //! Calculate cross product between two vectors.
+ /**
+ * A static method to calculate cross product between two vectors.
+ * Example:
+ * @code{c++}
+ * fp::vec3 v1{1,0,0};
+ * fp::vec3 v2{0,1,0};
+ * fp::vec3 v3 = cross(v1, v2); // v3 will contain {0,0,1}
+ * @endcode
+ * @param a first vector of the cross product
+ * @param b second vector of the cross product
+ * @return result of the cross product between the original vector and `v`.
+ */
+ static inline vec3 cross(const vec3 &a, const vec3 &b) {
+ vec3 res;
+ std::valarray as = a.data_.cshift(1);
+ std::valarray bs = b.data_.cshift(-1);
+
+ res.data_ = as * bs;
+
+ as = as.cshift(1);
+ bs = bs.cshift(-1);
+
+ res.data_ -= as * bs;
+
+ return res;
+ }
+
+ //! Calculate cross product with another vector.
+ /**
+ * Example:
+ * @code{c++}
+ * fp::vec3 v1{1,0,0};
+ * fp::vec3 v2{0,1,0};
+ * fp::vec3 v3 = v1.cross(v2); // v3 will contain {0,0,1}
+ * @endcode
+ * @param other the other vec3 vector.
+ * @return result of the cross product between the original vector and `other`.
+ */
+ vec3 cross(const vec3 &other) const { return cross(*this, other); }
+
+ //! Returns the norm of the vector.
+ /**
+ * Example:
+ * @code{c++}
+ * fp::vec3 v{1,0,1};
+ * double res = v.norm(); // res will contain 1,4142135624... i.e. sqrt(2)
+ * @endcode
+ * @return The euclidian norm of the vector.
+ */
+ Real norm() const { return std::sqrt(this->dot(*this)); }
+
+ //! Returns the square of the norm of the vector.
+ /**
+ * Example:
+ * @code{c++}
+ * fp::vec3 v{1,0,1};
+ * double res = v.norm_square(); // res will contain 2
+ * @endcode
+ * @return Square of the euclidian norm of the vector.
+ */
+ Real norm_square() const { return this->dot(*this); }
+
+ //! Normalize the vector in place. And return a reference to the new normalized vector.
+ /**
+ * @warning If you normalize a zero length vector, you effectively
+ * demand to divide by zero! this function will not do a security check
+ * for you and will just return nan!
+ * @return Reference to the normalized vector.
+ */
+ const vec3 &normalize() {
+ *this = *this / this->norm();
+ return *this;
+ }
+
+ //! Streaming operator for easy printing of the vector.
+ friend std::ostream &operator<<(std::ostream &os, const vec3 &obj) {
+ os << "{" << obj.x() << ',' << obj.y() << ',' << obj.z() << '}';
+ return os;
+ }
+
+ //! default equality operator.
+ /**
+ * @param other vec3 on the right hand side of the comparison operator.
+ * @return `true` if all elements of the compared vectors are equal and to `false` otherwise.
+ */
+ bool operator==(const vec3 &other) const = default;
+
+ //! Overloaded operator defined in terms of vec2::add.
+ /**
+ *
+ * @param lhs left hand side of the `+` operator
+ * @param rhs right hand side oif the `+` operator
+ * @return equivalent to a new copy of `lhs.add(rhs)`.
+ */
+ friend vec3 operator+(vec3 lhs, const vec3 &rhs) {
+ lhs += rhs;
+ return lhs;
+ }
+
+ //! Overloaded operator defined in terms of vec3::add.
+ /**
+ * Equivalent to `lhs.add(rhs)`.
+ * @param lhs left hand side of the `+=` operator
+ * @param rhs right hand side oif the `+=` operator
+ */
+ friend void operator+=(vec3 &lhs, const vec3 &rhs) { lhs.add(rhs); }
+
+ //! Overloaded operator defined in terms of vec3::subtract.
+ /**
+ *
+ * @param lhs left hand side of the `-` operator
+ * @param rhs right hand side oif the `-` operator
+ * @return equivalent to a new copy of `lhs.subtract(rhs)`.
+ */
+ friend vec3 operator-(vec3 lhs, const vec3 &rhs) {
+ lhs -= rhs;
+ return lhs;
+ }
+
+ //! Overloaded operator defined in terms of vec3::subtract.
+ /**
+ * Equivalent to `lhs.subtract(rhs)`.
+ * @param lhs left hand side of the `-=` operator
+ * @param rhs right hand side oif the `-=` operator
+ */
+ friend void operator-=(vec3 &lhs, const vec3 &rhs) { lhs.subtract(rhs); }
+
+ //! Overloaded operator defined in terms of vec3::scale.
+ /**
+ * Left multiplication by a scalar `s*v`.
+ * @param lhs left hand side of the `*` operator
+ * @param rhs right hand side oif the `*` operator
+ * @return equivalent to a new copy of `rhs.scale(lhs)`.
+ */
+ friend vec3 operator*(const Real &lhs, vec3 rhs) {
+ rhs.scale(lhs);
+ return rhs;
+ }
+
+ //! Overloaded operator defined in terms of vec3::scale.
+ /**
+ * Right multiplication by a scalar `v*s`.
+ * @param lhs left hand side of the `*` operator
+ * @param rhs right hand side oif the `*` operator
+ * @return equivalent to a new copy of `lhs.scale(rhs)`.
+ */
+ friend vec3 operator*(vec3 lhs, const Real &rhs) {
+ lhs.scale(rhs);
+ return lhs;
+ }
+
+ //! Overloaded operator defined in terms of vec3::scale.
+ /**
+ * In place division by a scalar `v/s`, equivalent to `lhs.scale(1/rhs)`.
+ * @param lhs left hand side of the `/=` operator
+ * @param rhs right hand side oif the `/=` operator
+ * @warning for performance reasons, this function will not check for zero division!
+ */
+ friend void operator/=(vec3 &lhs, const Real &rhs) { lhs.scale((Real)1. / rhs); }
+
+ //! Overloaded operator defined in terms of vec3::scale.
+ /**
+ * Division by a scalar `v/s`.
+ * @param lhs left hand side of the `/` operator
+ * @param rhs right hand side oif the `/` operator
+ * @return equivalent to a new copu of `lhs.scale(1/rhs)`.
+ * @warning for performance reasons, this function will not check for zero division!
+ */
+ friend vec3 operator/(vec3 lhs, const Real &rhs) {
+ lhs /= rhs;
+ return lhs;
+ }
+
+ //! element access operator.
+ /**
+ * @tparam Index automatically deduced type of the index.
+ * @param idx can only be 0 1 or 2. Any other number will cause the program to exit with an error.
+ * @return for a vec3 v: v[1] returns v.x, v[2] returns v.y and v[3] returns v.z.
+ *
+ * @note: The use of the subscription operator might be slower than the direct access of the data member.
+ */
+ template
+ requires std::is_integral_v
+ Real &operator[](Index idx) {
+ assert(idx < 3);
+ return data_[idx];
+ }
+
+ //! element access operator for constant environments.
+ /**
+ * @tparam Index automatically deduced type of the index.
+ * @param idx can only be 0 1 or 2. Any other number will cause the program to exit with an error.
+ * @return for a vec3 v: v[1] returns a constant reference to v.x, v[2] returns a constant reference to v.y and v[3]
+ * returns a constant reference to v.z.
+ *
+ * @note: The use of the subscription operator might be slower than the direct access of the data member.
+ */
+ template
+ requires std::is_integral_v
+ const Real &operator[](Index idx) const {
+ assert(idx < 3);
+ return data_[idx];
+ }
+
+ //! Unary minus operator.
+ /**
+ *
+ * @param v original vector.
+ * @return A copy of -v the vector v itself stays unaffected.
+ */
+ friend vec3 operator-(vec3 v) {
+ v.data_ = -v.data_;
+ return v;
+ }
+
+ //! Unary minus operator for rvalues.
+ /**
+ *
+ * @param v an rvalue vec3 vector.
+ * @return The rvalue vector `v` is moved into the function and `-v` is returned.
+ */
+ friend vec3 operator-(vec3 &&v) {
+ v.data_ = -v.data_;
+ return v;
+ }
+};
+} // namespace fp::experimental
+
+#endif // FLIPPY_VEC3_HPP
diff --git a/benchmarks/external/code_utils.hpp b/benchmarks/external/code_utils.hpp
new file mode 100644
index 0000000..392c5bc
--- /dev/null
+++ b/benchmarks/external/code_utils.hpp
@@ -0,0 +1,279 @@
+#pragma once
+#include
+#include
+#include
+#include
+#include
+
+namespace cutils {
+
+#if LOGGING_ON
+#define LOGN(x) \
+ std::cout << #x << ": "; \
+ cutils::print(x)
+#else
+#define LOGN(x)
+#endif
+
+#if LOGGING_ON
+#define LOG(...) cutils::print(__VA_ARGS__)
+#else
+#define LOG(...) ;
+#endif
+
+#if PRINTING_ON
+#define PRINT(...) cutils::print(__VA_ARGS__)
+#else
+#define PRINT(...)
+#endif
+
+template
+concept Beginable = requires(C c) { std::begin(c); };
+template
+concept Endable = requires(C c) { std::end(c); };
+template
+concept NeqableBeginAndEnd = requires(C c) {
+ { std::begin(c) != std::end(c) } -> std::same_as;
+};
+template
+concept BeginIncrementable = requires(C c) { std::begin(c)++; };
+template
+concept BeginDerefable = requires(C c) { *std::begin(c); };
+template
+concept BeginDerefToVoid = requires(C c) {
+ { *std::begin(c) } -> std::same_as;
+};
+
+template
+concept BeginAndEndCopyConstructibleAndDestructible = requires(C c) {
+ requires std::destructible && std::destructible &&
+ std::copy_constructible && std::copy_constructible;
+};
+
+template
+concept Container = Beginable && Endable && NeqableBeginAndEnd && BeginIncrementable && BeginDerefable &&
+ !BeginDerefToVoid && BeginAndEndCopyConstructibleAndDestructible;
+
+template
+concept StdOutStreamable = requires(T t) { std::cout << t; };
+template
+concept OutStreamable = requires(std::ostream os, T t) { os << t; };
+template
+concept InStreamable = requires(std::istream os, T t) { os >> t; };
+
+template
+concept Printable = StdOutStreamable || Container;
+
+template static void print();
+template static void print(const Container auto &output);
+template static void print(const Printable auto &first, const Printable auto &...rest);
+
+struct HumanReadableTime {
+ std::string unit;
+ std::string unit_fine;
+ long long diff;
+ long long diff_fine;
+ long long diff_ns;
+};
+
+// ------------------------ IMPLEMENTATIONS ---------------------------- //
+
+template [[maybe_unused]] void print() { std::cout << end; }
+
+template [[maybe_unused]] void print(const Container auto &output) {
+ print('{');
+ const auto &last_elem = *(output.end() - 1);
+ for (const auto &elem : output) {
+ if (elem != last_elem) {
+ print(elem, ',');
+ } else {
+ print(elem, '}');
+ }
+ }
+}
+
+/**
+ * Print function that can take an arbitrary number of arguments. It was implemented to work similarly to Pythons print
+ * function.
+ * @tparam Tfirst Type of the first argument (automatically deduced)
+ * @tparam Trest Type of the remaining part of the variadic number of arguments (automatically deduced)
+ * @param first First argument
+ * @param rest rest of the variadic arguments
+ * @returns nothing.
+ * @see print(const C>& output)
+ * @attention can not handle containers of containers.
+ * @warning if a container of zero length is provided the code will crush.
+ */
+// template
+//[[maybe_unused]] void print(const Tfirst& first, const Trest& ... rest)
+template
+[[maybe_unused]] void print(const Printable auto &first, const Printable auto &...rest) {
+ if constexpr (StdOutStreamable) {
+ std::cout << first << sep;
+ } else {
+ print(first);
+ }
+ print(rest...);
+}
+
+/**
+ * given a time difference in the highest clock resolution this struct constructs a human readable string of elapsed
+ * time.
+ */
+static HumanReadableTime human_readable_time(long long diff) {
+ std::string unit;
+ std::string unit_fine;
+ long long diff_ns = diff;
+ long long diff_fine = diff;
+ double diff_d = static_cast(diff);
+ if (diff_d / (1.e3l) < 1.l) {
+ unit = " ns";
+ } else if (diff_d / 1.e6 < 1.) {
+ unit = " µs";
+ diff /= 1000;
+ } else if (diff_d / 1.e9 < 1.) {
+ unit = " ms";
+ diff /= 1000'000ll;
+ unit_fine = " µs";
+ diff_fine /= 1000ll;
+ } else if (diff_d / (1.e9l * 60.l) < 1.l) {
+ unit = " s";
+ diff /= 1000'000'000ll;
+ unit_fine = " ms";
+ diff_fine /= 1000'000ll;
+ } else if (diff_d / (1.e9 * 3600.) < 1.) {
+ unit = " m";
+ diff /= 60'000'000'000ll;
+ unit_fine = " s";
+ diff_fine /= 1000'000'000ll;
+ } else {
+ unit = " h";
+ diff /= 3600'000'000'000ll;
+ unit_fine = " m";
+ diff_fine /= 60'000'000ll;
+ }
+ return {.unit = unit, .unit_fine = unit_fine, .diff = diff, .diff_fine = diff_fine, .diff_ns = diff_ns};
+}
+
+class Timer {
+ /**
+ * this class keeps time form its creation to destruction. I.e. it can
+ * time the duration of a scope.
+ *
+ * */
+ private:
+ using Clock = std::conditional_t;
+ Clock::time_point Start = Clock::now();
+ Clock::time_point Now;
+ bool stopped = false;
+
+ public:
+ [[maybe_unused]] Timer() = default;
+ void restart() { Start = Clock::now(); }
+ HumanReadableTime stop() {
+ Now = Clock::now();
+ long long diff_ = std::chrono::duration_cast(Now - Start).count();
+
+ HumanReadableTime hrt = human_readable_time(diff_);
+
+ auto end = std::chrono::system_clock::now();
+ auto end_time = std::chrono::system_clock::to_time_t(end);
+ std::cout << "finished computation at " << std::ctime(&end_time) << "elapsed time: " << hrt.diff << hrt.unit
+ << " (" << hrt.diff_fine << hrt.unit_fine << ")" << '\n';
+ stopped = true;
+ return hrt;
+ }
+
+ ~Timer() {
+ if (not stopped) { stop(); }
+ }
+};
+
+template
+concept UniformRandomBitGenerator = requires(G g) {
+ typename G::result_type; // Must have result_type
+ { G::min() } -> std::same_as; // min function
+ { G::max() } -> std::same_as; // max function
+ { g() } -> std::same_as; // operator()
+ requires std::is_unsigned_v; // result_type is unsigned integer type
+ requires G::min() < G::max(); // min is strictly less than max
+};
+
+template
+concept RandomNumberEngine = requires(E e,
+ const E x,
+ const E y,
+ typename E::result_type s,
+ unsigned long long z,
+ std::ostream &os,
+ std::istream &is) {
+ typename E::result_type;
+ { E() } -> std::same_as;
+ { E(x) } -> std::same_as;
+ { E(s) } -> std::same_as;
+ { e.seed() } -> std::same_as;
+ { e.seed(s) } -> std::same_as;
+ { e() } -> std::same_as;
+ { e.discard(z) } -> std::same_as;
+ { x == y } -> std::same_as;
+ { x != y } -> std::same_as;
+ { os << x } -> std::same_as;
+ { is >> e } -> std::same_as;
+ { E::min() } -> std::same_as;
+ { E::max() } -> std::same_as;
+};
+
+// source https://en.wikipedia.org/wiki/Xorshift#xoshiro_and_xoroshiro
+// template
+// requires std::is_integral::value
+// requires std::is_same_v
+class xorshift32 {
+ using ResultType = uint32_t;
+ ResultType seed_;
+
+ public:
+ typedef ResultType result_type;
+
+ xorshift32() : seed_(12) {}
+
+ explicit xorshift32(ResultType seed) : seed_(seed) {}
+
+ /* The state word must be initialized to non-zero */
+ ResultType operator()() {
+ /* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */
+ ResultType x = seed_;
+ x ^= x << 13;
+ x ^= x >> 17;
+ x ^= x << 5;
+ return seed_ = x;
+ }
+ void discard(long long z) {
+ for (long long i = 0; i < z; ++i) {
+ (*this)();
+ }
+ }
+ constexpr static ResultType min() { return std::numeric_limits::min(); }
+ constexpr static ResultType max() { return std::numeric_limits::max(); }
+
+ void seed() { *this = xorshift32(); }
+ void seed(ResultType seed_inp) { *this = xorshift32(seed_inp); }
+
+ bool operator==(const xorshift32 &rhs) const { return seed_ == rhs.seed_; }
+ bool operator!=(const xorshift32 &rhs) const { return seed_ != rhs.seed_; }
+
+ friend auto operator<<(std::ostream &os, const xorshift32 &rng) -> std::ostream & {
+ os.flags(std::ostream::dec | std::ostream::skipws);
+ os << rng.seed_;
+ return os;
+ }
+
+ friend auto operator>>(std::istream &is, xorshift32 &rng) -> std::istream & {
+ is.flags(std::istream::dec | std::istream::skipws);
+ is >> rng.seed_;
+ return is;
+ }
+};
+
+} // namespace cutils
diff --git a/benchmarks/make_benchmark_timeline.py b/benchmarks/make_benchmark_timeline.py
new file mode 100755
index 0000000..6b16735
--- /dev/null
+++ b/benchmarks/make_benchmark_timeline.py
@@ -0,0 +1,126 @@
+import os
+import sys
+from typing import List, Dict, NamedTuple
+from collections import namedtuple
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import plotnine as p9
+import yaml
+
+LogPathData = namedtuple('LogData', 'path timestamp')
+
+ParsedData = Dict[str, str]
+
+
+class LogData(NamedTuple):
+ timestamp: str
+ parsed_data: ParsedData
+
+
+def get_subdirectories(directory: str) -> ([str], [str]):
+ sub_dir_paths = []
+ sub_dir_names = []
+ for d in os.listdir(directory):
+ full_path = os.path.join(directory, d)
+ if os.path.isdir(full_path):
+ sub_dir_paths.append(full_path)
+ sub_dir_names.append(d)
+ return sub_dir_paths, sub_dir_names
+
+
+def get_files(directory: str, day: str) -> List[LogPathData]:
+ path_data: List[LogPathData] = []
+ for f in os.listdir(directory):
+ full_path = os.path.join(directory, f)
+ if os.path.isfile(full_path) and f.endswith('.yml'):
+ hour = f.split('.')[0]
+ data = LogPathData(path=full_path, timestamp=f'{day}-{hour}')
+ path_data.append(data)
+ return path_data
+
+
+def _parse_log(log: str)-> Dict[str, str]:
+ return yaml.safe_load(log)
+
+
+def parse_log(log_path_data: LogPathData):
+ log_data: LogData
+ with open(file=log_path_data.path, mode='r') as log_file:
+ raw_data = log_file.read()
+ log_data = LogData(parsed_data=_parse_log(raw_data),
+ timestamp=log_path_data.timestamp)
+ return log_data
+
+
+def get_all_logs(logs_dir: str,
+ benchmark_version: int | None = None) -> List[ParsedData]:
+
+ assert os.path.isdir(logs_dir), f'{logs_dir=} is not a valid path'
+
+ log_sub_dirs, log_dates = get_subdirectories(logs_dir)
+
+ all_log_paths: List[LogPathData] = []
+ for lsdp, day in zip(log_sub_dirs, log_dates):
+ files = get_files(lsdp, day)
+ all_log_paths += files
+
+ all_logs: List[ParsedData] = []
+ for log_path_data in all_log_paths:
+ log = parse_log(log_path_data=log_path_data)
+ if benchmark_version is not None:
+ if ("BENCHMARK VERSION" not in log.parsed_data or
+ log.parsed_data['BENCHMARK VERSION'] != benchmark_version):
+ continue
+ log.parsed_data['timestamp'] = log.timestamp
+ all_logs.append(log.parsed_data)
+ return all_logs
+
+
+def main():
+ logs_dir = sys.argv[1]
+ all_logs: List[ParsedData] = get_all_logs(logs_dir)
+
+ df = pd.DataFrame(all_logs)
+ def version_num_to_category(num: float):
+ if not pd.isna(num):
+ return str(int(num))
+ else:
+ return '-'
+
+ df['BENCHMARK VERSION'] = df['BENCHMARK VERSION'].apply(version_num_to_category)
+ df['timestamp'] = pd.to_datetime(df['timestamp'], format='%Y-%m-%d-%H-%M-%S')
+ df.sort_values(by='timestamp')
+
+
+ df_long = pd.melt(df, id_vars=['timestamp', 'BENCHMARK VERSION'],
+ value_vars=['mps', 'emps', 'fps', 'efps'],
+ var_name='perf_type', value_name='perf_value')
+ df_long.sort_values(by='timestamp')
+ df_long['timestamp'] = df_long['timestamp'].dt.strftime('%Y-%m-%d %H:%M')
+ print(df_long.tail())
+
+ # x_breaks = df_long['timestamp'][::2].to_list()
+ # print(f'{x_breaks=}')
+ plot = (p9.ggplot(df_long)
+ + p9.aes(x='timestamp', y='perf_value', color='BENCHMARK VERSION')
+ # + p9.geom_text(p9.aes(label='BENCHMARK VERSION'), size=9)
+ + p9.geom_point()
+ + p9.facet_wrap('~perf_type', scales='free_y')
+ + p9.theme_minimal()
+ + p9.labs(x='time', y='perf', color='Bench.\nVersion',
+ title='performance over time')
+ # + p9.scale_color_discrete(name='Bench.\nV')
+ # + p9.scale_x_discrete(breaks=x_breaks)
+ # + p9.scale_x_datetime(breaks=x_breaks)
+ + p9.theme(
+ figure_size=(15, 9),
+ axis_text_x=p9.element_text(angle=70, hjust=1),
+ )
+ )
+
+ print(plot)
+
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file
diff --git a/benchmarks/micro_benchmarks/bench_diff.py b/benchmarks/micro_benchmarks/bench_diff.py
new file mode 100755
index 0000000..846f06f
--- /dev/null
+++ b/benchmarks/micro_benchmarks/bench_diff.py
@@ -0,0 +1,87 @@
+#! python3
+import sys
+import json
+import re
+import numpy as np
+from termcolor import colored, cprint
+
+if __name__ == '__main__':
+ bench_types = dict()
+ path = sys.argv[1]
+ bench_file = sys.argv[2]
+ with open(f"{path}/{bench_file}", 'r') as f:
+ j = json.load(f)
+ reg = re.compile(r'(^[a-zA-Z][\w\s]*)(base|experimental)', re.IGNORECASE)
+
+ for benchmark in j['benchmarks']:
+ full_name = benchmark["name"]
+ name = reg.match(full_name).group(1)
+ suffix = reg.match(full_name).group(2)
+ if not bench_types.get(name):
+ bench_types[name] = dict()
+ if benchmark['aggregate_name'] == 'mean':
+ bench_types[name][suffix] = {"mean": benchmark['real_time']}
+ if benchmark['aggregate_name'] == 'stddev':
+ bench_types[name][suffix]["std_dev"] = benchmark['real_time']
+
+ diff_log_name = "bench_diff.log"
+ diff_log_file = open(f"{path}/bench_diff.log", "w")
+ diff_json_file = open(f"{path}/bench_diff.json", "w")
+
+ diff_dict = dict()
+
+ string_stencil = f"{'test name':>30}{' ':<3}{'base time':<9} | {'expr time':<9} | {'rel impr':<9}\n"
+ line_width=len(string_stencil)
+ string_stencil += f"{' ':>30}{' ':<3}{'base std':<9} | {'expr std':<9} | {'av_error/tot_diff':<9}\n"
+ string_stencil += line_width*'_'+'\n'
+ for benchmark_name, benchmark in bench_types.items():
+ base_time = benchmark['base']['mean']
+ base_std = benchmark['base']['std_dev']
+ experimental_time = benchmark['experimental']['mean']
+ experimental_std = benchmark['base']['std_dev']
+ diff = base_time - experimental_time
+ av_noise = (base_std + experimental_std)/2
+ rel_change = diff/base_time
+ rel_noise = av_noise / abs(diff)
+
+ if not diff_dict.get(benchmark_name):
+ diff_dict[benchmark_name] = dict()
+
+ diff_dict[benchmark_name]["base time"] = base_time
+ diff_dict[benchmark_name]["base std"] = base_std
+
+ diff_dict[benchmark_name]["experimental time"] = experimental_time
+ diff_dict[benchmark_name]["experimental std"] = experimental_std
+
+ diff_dict[benchmark_name]["diff"] = diff
+ diff_dict[benchmark_name]["av_noise"] = av_noise
+
+ diff_dict[benchmark_name]["rel change"] = rel_change
+ diff_dict[benchmark_name]["rel error"] = rel_noise
+
+ attr = ['bold', 'dark']
+ if rel_noise < 1:
+ if diff > 0:
+ ticker = '∆∆∆'
+ color = 'cyan'
+ if rel_noise<0.1:
+ color = 'green'
+ attr.append('blink')
+ else:
+ ticker = 'VVV'
+ color = 'red'
+ attr.append('blink')
+ else:
+ ticker = '---'
+ color = 'yellow'
+
+ res = f'{np.round(rel_change,3)} {colored(ticker, color, attrs=attr)}'
+ string_stencil += f"{benchmark_name:>30}{' ':<3}{np.round(base_time,3):<9} | {np.round(experimental_time,3):<9} | {res:<9}\n"
+ string_stencil += f"{' ':>30}{' ':<3}{np.round(base_std,3):<9} | {np.round(experimental_std,3):<9} | {np.round(rel_noise,3):<9}\n"
+ string_stencil += int(line_width/4)*' -- '+'\n'
+
+ print(string_stencil, file=diff_log_file)
+ print(string_stencil)
+
+ out_bench = {'context': j['context'], 'diffs': diff_dict }
+ json.dump(out_bench, diff_json_file, indent = 6)
diff --git a/benchmarks/micro_benchmarks/vec3_bench/CMakeLists.txt b/benchmarks/micro_benchmarks/vec3_bench/CMakeLists.txt
new file mode 100755
index 0000000..eb036d3
--- /dev/null
+++ b/benchmarks/micro_benchmarks/vec3_bench/CMakeLists.txt
@@ -0,0 +1,48 @@
+cmake_minimum_required(VERSION 3.25)
+project(vec3_microbenchmark)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CXX_STANDARD_REQUIRED ON)
+
+set(CMAKE_BUILD_TYPE Release)
+
+add_compile_definitions(TEST_ASSET_PATH="${CMAKE_CURRENT_SOURCE_DIR}/assets/")
+if(${CMAKE_CXX_COMPILER_ID} STREQUAL MSVC)
+ message("building for MSVC")
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O2")
+else()
+ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
+endif()
+message(${CMAKE_CXX_COMPILER_ID})
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -g")
+
+include_directories(.)
+include_directories(external)
+
+FetchContent_MakeAvailable(Catch2)
+
+list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
+include(CTest)
+include(Catch)
+
+# Create an interface library to suppress warnings for Catch2
+
+message("testing FLIPPY")
+include_directories(../flippy)
+include_directories(../flippy/external)
+
+message(${catch2_SOURCE_DIR})
+add_executable(
+ ${PROJECT_NAME} vec3_bench.cpp # node_neighbour_insertion_delition_bench.cpp
+)
+
+catch_discover_tests(${PROJECT_NAME})
+
+target_link_libraries(${PROJECT_NAME} PRIVATE Catch2::Catch2WithMain)
+target_include_directories(
+ ${PROJECT_NAME} SYSTEM PRIVATE ${catch2_SOURCE_DIR}/src
+ ${catch2_BINARY_DIR}/src)
+target_compile_options(${PROJECT_NAME} PRIVATE -Wall -Wextra -Wpedantic
+ -Wshadow -Wconversion -Werror)
+enable_testing()
+add_test(${PROJECT_NAME} ${PROJECT_NAME})
diff --git a/benchmarks/micro_benchmarks/vec3_bench/node_neighbour_insertion_delition_bench.cpp b/benchmarks/micro_benchmarks/vec3_bench/node_neighbour_insertion_delition_bench.cpp
new file mode 100755
index 0000000..c4dfe3a
--- /dev/null
+++ b/benchmarks/micro_benchmarks/vec3_bench/node_neighbour_insertion_delition_bench.cpp
@@ -0,0 +1,65 @@
+#include "flippy.hpp"
+#include
+#include
+#include
+
+using REAL = double;
+using INDEX = unsigned int;
+
+#define SETUP \
+ fp::vec3 v1{0., 1., 1.}; \
+ fp::vec3 v2{1., 0., 7.}; \
+ fp::vec3 v3{1., 0., 0.}; \
+ fp::vec3 v4{0., 2., 3.}; \
+ fp::vec3 v5{0., 0., 9.}; \
+ fp::vec3 v6{1., 1., 3.}; \
+ fp::vec3 v7{0., 3., 3.}; \
+ fp::vec3 v8{1., 0., 3.}; \
+ Node node0; \
+ node0.id = 0; \
+ node0.nn_ids = {1, 2, 3, 4, 5, 6, 7, 8}; \
+ node0.nn_distances = {v1, v2, v3, v4, v5, v6, v7, v8};
+
+#define BASE_EXPERIMENTAL_COMPARE(bench) \
+ BENCHMARK_ADVANCED("base")(Catch::Benchmark::Chronometer meter) { \
+ typedef fp::Node Node; \
+ SETUP \
+ bench \
+ }; \
+ BENCHMARK_ADVANCED("experimental")(Catch::Benchmark::Chronometer meter) { \
+ typedef fp::experimental::Node Node; \
+ SETUP \
+ bench \
+ };
+
+TEST_CASE("Benchmark node", "[!benchmark]") {
+
+ INDEX new_id{11};
+ fp::vec3 new_pos{0., 0., 0};
+
+ SECTION("emplace at 0") {
+#define emplace_at_0 \
+ meter.measure([&] { \
+ return node0.emplace_nn_id(new_id, new_pos, 0); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(emplace_at_0)
+ };
+
+ SECTION("emplace at the end") {
+#define emplace_at_the_end \
+ auto end = static_cast(node0.nn_ids.size() - 1); \
+ meter.measure([&] { \
+ return node0.emplace_nn_id(new_id, new_pos, end); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(emplace_at_the_end)
+ };
+
+ SECTION("emplace in the middle") {
+#define emplace_in_the_middle \
+ auto end = static_cast(node0.nn_ids.size() / 2); \
+ meter.measure([&] { \
+ return node0.emplace_nn_id(new_id, new_pos, end); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(emplace_in_the_middle)
+ };
+}
diff --git a/benchmarks/micro_benchmarks/vec3_bench/vec3_bench.cpp b/benchmarks/micro_benchmarks/vec3_bench/vec3_bench.cpp
new file mode 100755
index 0000000..ce116d3
--- /dev/null
+++ b/benchmarks/micro_benchmarks/vec3_bench/vec3_bench.cpp
@@ -0,0 +1,130 @@
+#include "../experimental/vec3_valarray.h"
+#include "experimental/vec3_valarray.h"
+#include "vec3.hpp"
+#include
+#include
+#include
+
+using REAL = double;
+constexpr int numTrials = 1;
+constexpr REAL min = -1e6, max = 1e6;
+
+#define TWO_VEC_SETUP \
+ const std::vector rand = GENERATE(take(numTrials, chunk(6, random(min, max)))); \
+ [[maybe_unused]] vec3 v0{rand[0], rand[1], rand[2]}; \
+ [[maybe_unused]] vec3 v1{rand[3], rand[4], rand[5]};
+
+#define BASE_EXPERIMENTAL_COMPARE(bench) \
+ BENCHMARK_ADVANCED("base")(Catch::Benchmark::Chronometer meter) { \
+ typedef fp::vec3 vec3; \
+ TWO_VEC_SETUP \
+ bench \
+ }; \
+ BENCHMARK_ADVANCED("experimental")(Catch::Benchmark::Chronometer meter) { \
+ typedef fp::experimental::vec3 vec3; \
+ TWO_VEC_SETUP \
+ bench \
+ };
+
+TEST_CASE("Benchmark vec3", "[!benchmark]") {
+
+ SECTION("named creation") {
+ const std::vector rand = GENERATE(take(numTrials, chunk(6, random(min, max))));
+ BENCHMARK_ADVANCED("base")(Catch::Benchmark::Chronometer meter) {
+ using fp::vec3;
+ meter.measure([&rand] {
+ vec3 v1{rand[3], rand[4], rand[5]};
+ return v1;
+ });
+ };
+
+ BENCHMARK_ADVANCED("experimental")(Catch::Benchmark::Chronometer meter) {
+ using fp::experimental::vec3;
+ meter.measure([&rand] {
+ vec3 v1{rand[3], rand[4], rand[5]};
+ return v1;
+ });
+ };
+ }
+
+ SECTION("unnamed creation") {
+ const std::vector rand = GENERATE(take(numTrials, chunk(6, random(min, max))));
+ BENCHMARK_ADVANCED("base") { return fp::vec3{rand[3], rand[4], rand[5]}; };
+ BENCHMARK_ADVANCED("experimental") { return fp::experimental::vec3{rand[3], rand[4], rand[5]}; };
+ }
+
+ SECTION("dot") {
+#define dot_bench \
+ meter.measure([&] { \
+ return v0.dot(v1); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(dot_bench)
+ }
+
+ SECTION("norm") {
+#define norm_bench \
+ meter.measure([&] { \
+ return v0.norm(); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(norm_bench)
+ }
+
+ SECTION("norm_square") {
+#define norm_square_bench \
+ meter.measure([&] { \
+ return v0.norm_square(); \
+ });
+ BASE_EXPERIMENTAL_COMPARE(norm_square_bench)
+ }
+
+ SECTION("two vector subtract") {
+#define subtraction_bench \
+ meter.measure([&] { \
+ return v0 - v1; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(subtraction_bench)
+ };
+
+ SECTION("vector minus") {
+#define vector_minus_bench \
+ meter.measure([&] { \
+ return -v1; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(vector_minus_bench)
+ };
+
+ SECTION("two vector add") {
+#define add_bench \
+ meter.measure([&] { \
+ return v0 + v1; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(add_bench)
+ };
+
+ SECTION("multiply left") {
+#define multiply_left_bench \
+ auto r = rand[4]; \
+ meter.measure([&] { \
+ return r * v1; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(multiply_left_bench)
+ };
+
+ SECTION("multiply right") {
+#define multiply_right_bench \
+ auto r = rand[4]; \
+ meter.measure([&] { \
+ return v1 * r; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(multiply_left_bench)
+ };
+
+ SECTION("complex expression") {
+#define complex_expression_bench \
+ auto r = rand[4]; \
+ meter.measure([&] { \
+ return (v1 * r + v0 * v1.dot(v0)).cross(v0) - v0; \
+ });
+ BASE_EXPERIMENTAL_COMPARE(complex_expression_bench)
+ };
+}
diff --git a/benchmarks/rbc_bench/CMakeLists.txt b/benchmarks/rbc_bench/CMakeLists.txt
new file mode 100755
index 0000000..e495b2b
--- /dev/null
+++ b/benchmarks/rbc_bench/CMakeLists.txt
@@ -0,0 +1,17 @@
+cmake_minimum_required(VERSION 3.16)
+project(flippy_rbc_benchmark)
+set(CMAKE_CXX_STANDARD 20)
+#set(CMAKE_BUILD_TYPE Debug)
+set(CMAKE_BUILD_TYPE Release)
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -fno-omit-frame-pointer")
+
+add_definitions(-DLOGGING_ON=1)
+add_definitions(-DPRINTING_ON=1)
+add_definitions(-DALLOCATIONTRACKING=1)
+
+include_directories(../benchmark_logger)
+include_directories(../external)
+include_directories(../../flippy)
+add_executable(${PROJECT_NAME} main.cpp)
+target_compile_options(${PROJECT_NAME} PRIVATE -Wall -Wextra -Wpedantic -Wshadow -Wconversion -Werror -ffast-math)
\ No newline at end of file
diff --git a/benchmarks/rbc_bench/main.cpp b/benchmarks/rbc_bench/main.cpp
new file mode 100755
index 0000000..f5aa5a2
--- /dev/null
+++ b/benchmarks/rbc_bench/main.cpp
@@ -0,0 +1,215 @@
+#if ALLOCATIONTRACKING
+
+#include
+#include
+#include
+
+namespace {
+struct AllocMetrics {
+ std::int64_t allocation = 0;
+ std::uint64_t allocated_size = 0;
+ std::int64_t de_allocation = 0;
+ std::int64_t leaks = 0;
+};
+} // namespace
+
+static AllocMetrics AM;
+
+void *operator new(size_t size) {
+ AM.allocation += 1;
+ AM.allocated_size += size;
+ // std::cout<<"allocating "<<(long double)size/1000.<<" kB\n";
+ return std::malloc(size);
+}
+void operator delete(void *memory) noexcept {
+ AM.de_allocation += 1;
+ std::free(memory);
+}
+
+void operator delete(void *memory, std::size_t) noexcept {
+ AM.de_allocation += 1;
+ std::free(memory);
+}
+
+#endif
+
+#include
+#include
+
+struct EnergyParameters {
+ fp::Real kappa, K_V, K_A, V_t, A_t;
+};
+using Triangulation = fp::Triangulation;
+using MCU = fp::MonteCarloUpdater;
+
+// This is the energy function that is used by flippy's built-in updater to decide if a move was energetically favorable
+// or not
+fp::Real surface_energy(const fp::Node &node,
+ const fp::Triangulation &trg,
+ const EnergyParameters &prms,
+ const std::vector &changed_neighborhood) {
+ fp::Real V = trg.global_geometry().volume;
+ fp::Real A = trg.global_geometry().area;
+ fp::Real dV = V - prms.V_t;
+ fp::Real dA = A - prms.A_t;
+
+ fp::Real E_bending = fp::changed_neighborhood_bending_energy(node, trg.nodes(), changed_neighborhood);
+
+ fp::Real energy = prms.kappa * E_bending + prms.K_V * dV * dV / prms.V_t + prms.K_A * dA * dA / prms.A_t;
+ return energy;
+}
+
+int main(int /*argc*/, char **argv) {
+ char LOG_FILE[300]{};
+ {
+ // triangulation iteration number of nodes
+ auto json_path = std::string(argv[1]);
+ fp::Json config = fp::json_read(json_path);
+ // N_node=12+30*n+20*n*(n-1)/2 where n is the same as n_trng
+ unsigned int n_triang = config["n_triang"].get();
+ fp::Real l_min = config["l_min"].get();
+ fp::Real kappa = config["kappa"].get(); /*kBT*/
+ fp::Real K_A = config["K_A"].get(); /*kBT/volume*/
+ fp::Real K_V = config["K_V"].get(); /*kBT/area*/
+ fp::Real red_vol = config["red_vol"].get();
+ int max_mc_steps = config["max_mc_steps"].get();
+ std::string save_dir = config["save_dir"];
+
+ // estimate of a typical bond length in the initial triangulation and then create a sphere such that the initial
+ // bond length is close to minimal. This formula is derived from the equidistant sub-triangulation of an
+ // icosahedron, where geodesic distances are used as a distance measure.
+ fp::Real R = fp::min_radius_with_non_overlapping_beads(l_min, n_triang);
+ fp::Real l_max = 1.8 * l_min; // if you make l_max closer to l_min
+ // bond_flip acceptance rate will go down
+ fp::Real r_Verlet = 2 * l_max;
+
+ fp::Real Vt = red_vol * fp::sphere_vol(R);
+ fp::Real At = fp::sphere_area(R);
+ EnergyParameters prms{.kappa = kappa, .K_V = K_V, .K_A = K_A, .V_t = Vt, .A_t = At};
+ // side length of a voxel from which the displacement of the node is drawn
+ fp::Real linear_displ = l_min / 8;
+ // max number of iteration steps (depending on the strength of your CPU, this should take anywhere from a couple
+ // of seconds to a couple of minutes
+
+ std::random_device random_number_generator_seed;
+ // create a random number generator and seed it with the current time
+ std::mt19937 rng(random_number_generator_seed());
+
+ auto guv = Triangulation(n_triang, R, r_Verlet);
+ auto mc_updater = MCU(guv, prms, surface_energy, rng, l_min, l_max);
+
+ fp::vec3 displ{};
+ std::uniform_real_distribution displ_distr(-linear_displ, linear_displ);
+
+ // squish the sphere in the z-direction to break the initial symmetry. This speeds up the convergence to a
+ // biconcave shape greatly.
+ guv.scale_node_coordinates(1, 1, 0.8);
+
+ fp::Json data_init = guv.make_egg_data();
+ fp::json_dump(save_dir + "test_run_init", data_init);
+
+ auto globals_saver = [&](const Triangulation &trg) {
+ std::unordered_map out;
+
+ out["volume"] = std::to_string(trg.global_geometry().volume);
+ out["area"] = std::to_string(trg.global_geometry().area);
+ out["volume_rel"] = std::to_string(trg.global_geometry().volume / Vt);
+ out["area_rel"] = std::to_string(trg.global_geometry().area / At);
+ return out;
+ };
+
+ auto stream_particle = [&](const fp::Node &node, const fp::Triangulation &) {
+ std::string s{"1"};
+ fp::vec3 pos = node.pos;
+ static const std::string empty{" "};
+ s += empty + std::to_string(pos.x) + empty + std::to_string(pos.y) + empty + std::to_string(pos.z) + empty +
+ std::to_string(node.curvature_vec.norm()) + "\n";
+ return s;
+ };
+ std::vector properties{
+ fp::experimental::xyzProperty{
+ .name = "species", .xyz_type = fp::experimental::XYZ_TYPE::S, .column_count = 1},
+ fp::experimental::xyzProperty{.name = "pos", .xyz_type = fp::experimental::XYZ_TYPE::R, .column_count = 3},
+ fp::experimental::xyzProperty{
+ .name = "curv", .xyz_type = fp::experimental::XYZ_TYPE::R, .column_count = 1}};
+ auto xyz_saver = fp::experimental::xyzDataSaver(save_dir + "data.xyz", properties, stream_particle);
+ std::vector shuffled_ids;
+ shuffled_ids.reserve(guv.size());
+
+ // create a vector that contains all node ids. We can shuffle this vector in each MC step to iterate randomly
+ // through the nodes
+ for (const auto &node : guv.nodes()) {
+ shuffled_ids.push_back(node.id);
+ }
+
+ fp::Real probability_target = 0.5f;
+ fp::Real adaptation_magnitude = 0.1f;
+ std::optional e_diff = 0;
+
+ fp::Real V0 = guv.global_geometry().volume;
+ fp::Real A0 = guv.global_geometry().area;
+ auto displ_updater = fp::DynamicDisplacementUpdater(linear_displ, adaptation_magnitude, probability_target);
+ auto logger = FlippyBenchmarkLogger