From 7a1fe25e37e71a20df20784e40a2e450456d0081 Mon Sep 17 00:00:00 2001 From: Ankith Date: Thu, 30 May 2024 14:51:53 +0530 Subject: [PATCH 01/62] Switch to scikit_build_core, remove setup.py - Moves the metadata from setup.py to pyproject.toml --- backend/pyproject.toml | 49 ++++++++++----- backend/setup.py | 139 ----------------------------------------- 2 files changed, 32 insertions(+), 156 deletions(-) delete mode 100644 backend/setup.py diff --git a/backend/pyproject.toml b/backend/pyproject.toml index fd92373c6..a7620e73c 100644 --- a/backend/pyproject.toml +++ b/backend/pyproject.toml @@ -1,24 +1,39 @@ -[build-system] -requires = [ - "setuptools>=42", - "wheel", - "ninja", - "cmake>=3.12", +[project] +name = "elasticapp" +version = "0.0.3" +description = "A CPP accelerated backend for PyElastica kernels" +# readme = "README.rst" # for long description +requires-python = ">=3.10" +license = {text = "MIT License"} + +# this is a list that can be expanded +authors = [{name = "Seung Hyun Kim", email = "skim0119@gmail.com"}] + +# classifiers can be added here, later +classifiers = [ ] -build-backend = "setuptools.build_meta" -[tool.mypy] -files = "setup.py" -python_version = "3.7" -strict = true -show_error_codes = true -enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"] -warn_unreachable = true +[project.urls] +"Homepage" = "https://github.com/GazzolaLab/PyElastica/" +"Bug Reports" = "https://github.com/GazzolaLab/PyElastica/issues" +"Source" = "https://github.com/GazzolaLab/PyElastica/" +"Documentation" = "https://docs.cosseratrods.org/" -[[tool.mypy.overrides]] -module = ["ninja"] -ignore_missing_imports = true +[build-system] +requires = ["scikit-build-core>=0.3.3", "pybind11"] +build-backend = "scikit_build_core.build" +# [tool.mypy] +# files = "setup.py" +# python_version = "3.7" +# strict = true +# show_error_codes = true +# enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"] +# warn_unreachable = true +# +# [[tool.mypy.overrides]] +# module = ["ninja"] +# ignore_missing_imports = true [tool.pytest.ini_options] minversion = "6.0" diff --git a/backend/setup.py b/backend/setup.py deleted file mode 100644 index 80240487e..000000000 --- a/backend/setup.py +++ /dev/null @@ -1,139 +0,0 @@ -import os -import re -import subprocess -import sys -from pathlib import Path - -from setuptools import Extension, setup -from setuptools.command.build_ext import build_ext - -# Convert distutils Windows platform specifiers to CMake -A arguments -PLAT_TO_CMAKE = { - "win32": "Win32", - "win-amd64": "x64", - "win-arm32": "ARM", - "win-arm64": "ARM64", -} - - -# A CMakeExtension needs a sourcedir instead of a file list. -# The name must be the _single_ output extension from the CMake build. -# If you need multiple extensions, see scikit-build. -class CMakeExtension(Extension): - def __init__(self, name: str, sourcedir: str = "") -> None: - super().__init__(name, sources=[]) - self.sourcedir = os.fspath(Path(sourcedir).resolve()) - - -class CMakeBuild(build_ext): - def build_extension(self, ext: CMakeExtension) -> None: - # Must be in this form due to bug in .resolve() only fixed in Python 3.10+ - ext_fullpath = Path.cwd() / self.get_ext_fullpath(ext.name) - extdir = ext_fullpath.parent.resolve() - - # Using this requires trailing slash for auto-detection & inclusion of - # auxiliary "native" libs - - debug = int(os.environ.get("DEBUG", 0)) if self.debug is None else self.debug - cfg = "Debug" if debug else "Release" - - # CMake lets you override the generator - we need to check this. - # Can be set with Conda-Build, for example. - cmake_generator = os.environ.get("CMAKE_GENERATOR", "") - - # Set Python_EXECUTABLE instead if you use PYBIND11_FINDPYTHON - # EXAMPLE_VERSION_INFO shows you how to pass a value into the C++ code - # from Python. - cmake_args = [ - f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={extdir}{os.sep}", - f"-DPYTHON_EXECUTABLE={sys.executable}", - f"-DCMAKE_BUILD_TYPE={cfg}", # not used on MSVC, but no harm - ] - build_args = [] - # Adding CMake arguments set as environment variable - # (needed e.g. to build for ARM OSx on conda-forge) - if "CMAKE_ARGS" in os.environ: - cmake_args += [item for item in os.environ["CMAKE_ARGS"].split(" ") if item] - - # In this example, we pass in the version to C++. You might not need to. - cmake_args += [f"-DEXAMPLE_VERSION_INFO={self.distribution.get_version()}"] - - if self.compiler.compiler_type != "msvc": - # Using Ninja-build since it a) is available as a wheel and b) - # multithreads automatically. MSVC would require all variables be - # exported for Ninja to pick it up, which is a little tricky to do. - # Users can override the generator with CMAKE_GENERATOR in CMake - # 3.15+. - if not cmake_generator or cmake_generator == "Ninja": - try: - import ninja - - ninja_executable_path = Path(ninja.BIN_DIR) / "ninja" - cmake_args += [ - "-GNinja", - f"-DCMAKE_MAKE_PROGRAM:FILEPATH={ninja_executable_path}", - ] - except ImportError: - pass - - else: - # Single config generators are handled "normally" - single_config = any(x in cmake_generator for x in {"NMake", "Ninja"}) - - # CMake allows an arch-in-generator style for backward compatibility - contains_arch = any(x in cmake_generator for x in {"ARM", "Win64"}) - - # Specify the arch if using MSVC generator, but only if it doesn't - # contain a backward-compatibility arch spec already in the - # generator name. - if not single_config and not contains_arch: - cmake_args += ["-A", PLAT_TO_CMAKE[self.plat_name]] - - # Multi-config generators have a different way to specify configs - if not single_config: - cmake_args += [ - f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}" - ] - build_args += ["--config", cfg] - - if sys.platform.startswith("darwin"): - # Cross-compile support for macOS - respect ARCHFLAGS if set - archs = re.findall(r"-arch (\S+)", os.environ.get("ARCHFLAGS", "")) - if archs: - cmake_args += ["-DCMAKE_OSX_ARCHITECTURES={}".format(";".join(archs))] - - # Set CMAKE_BUILD_PARALLEL_LEVEL to control the parallel build level - # across all generators. - if "CMAKE_BUILD_PARALLEL_LEVEL" not in os.environ: - # self.parallel is a Python 3 only way to set parallel jobs by hand - # using -j in the build_ext call, not supported by pip or PyPA-build. - if hasattr(self, "parallel") and self.parallel: - # CMake 3.12+ only. - build_args += [f"-j{self.parallel}"] - - build_temp = Path(self.build_temp) / ext.name - if not build_temp.exists(): - build_temp.mkdir(parents=True) - - subprocess.run( - ["cmake", ext.sourcedir, *cmake_args], cwd=build_temp, check=True - ) - subprocess.run( - ["cmake", "--build", ".", *build_args], cwd=build_temp, check=True - ) - - -# The information here can also be placed in setup.cfg - better separation of -# logic and declaration, and simpler if you include description/version in a file. -setup( - name="elasticapp", - version="0.0.3", - author="Seung Hyun Kim", - author_email="skim0119@gmail.com", - description="A CPP accelerated backend for PyElastica kernels", - long_description="", - ext_modules=[CMakeExtension("elasticapp")], - cmdclass={"build_ext": CMakeBuild}, - zip_safe=False, - python_requires=">=3.10", -) From 060f94f2b7bcb207787d28732a46da8bd93da0c6 Mon Sep 17 00:00:00 2001 From: Ankith Date: Thu, 30 May 2024 15:33:50 +0530 Subject: [PATCH 02/62] Add basic pybind11 example with meson buildsystem --- backend/meson.build | 25 +++++++++++++++++++++ backend/pyproject.toml | 4 ++-- backend/src/example_1.cpp | 46 +++++++++++++++++++++++++++++++++++++++ backend/src/meson.build | 8 +++++++ 4 files changed, 81 insertions(+), 2 deletions(-) create mode 100644 backend/meson.build create mode 100644 backend/src/example_1.cpp create mode 100644 backend/src/meson.build diff --git a/backend/meson.build b/backend/meson.build new file mode 100644 index 000000000..62c509157 --- /dev/null +++ b/backend/meson.build @@ -0,0 +1,25 @@ +project( + 'elasticapp', + ['cpp'], + version: '0.0.3', + default_options: ['cpp_std=c++17'], +) + +package = 'elasticapp' + +cc = meson.get_compiler('cpp') + +py = import('python').find_installation(pure: false) +py_dep = py.dependency() + +# find dependencies and create dep objects +pybind11_dep = declare_dependency( + include_directories: run_command( + py, + '-c', + 'import pybind11; print(pybind11.get_include());', + check: true, + ).stdout().strip(), +) + +subdir('src') diff --git a/backend/pyproject.toml b/backend/pyproject.toml index a7620e73c..f10127b1b 100644 --- a/backend/pyproject.toml +++ b/backend/pyproject.toml @@ -20,8 +20,8 @@ classifiers = [ "Documentation" = "https://docs.cosseratrods.org/" [build-system] -requires = ["scikit-build-core>=0.3.3", "pybind11"] -build-backend = "scikit_build_core.build" +requires = ["meson-python", "pybind11"] +build-backend = "mesonpy" # [tool.mypy] # files = "setup.py" diff --git a/backend/src/example_1.cpp b/backend/src/example_1.cpp new file mode 100644 index 000000000..088cc43c6 --- /dev/null +++ b/backend/src/example_1.cpp @@ -0,0 +1,46 @@ +// Example taken from: https://github.com/pybind/python_example/blob/master/src/main.cpp +// To demonstate the workability of the meson-based builder + +#include + +#define STRINGIFY(x) #x +#define MACRO_STRINGIFY(x) STRINGIFY(x) + +int add(int i, int j) { + return i + j; +} + +namespace py = pybind11; + +PYBIND11_MODULE(example_1, m) { + m.doc() = R"pbdoc( + Pybind11 example plugin + ----------------------- + + .. currentmodule:: python_example + + .. autosummary:: + :toctree: _generate + + add + subtract + )pbdoc"; + + m.def("add", &add, R"pbdoc( + Add two numbers + + Some other explanation about the add function. + )pbdoc"); + + m.def("subtract", [](int i, int j) { return i - j; }, R"pbdoc( + Subtract two numbers + + Some other explanation about the subtract function. + )pbdoc"); + +#ifdef VERSION_INFO + m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); +#else + m.attr("__version__") = "dev"; +#endif +} \ No newline at end of file diff --git a/backend/src/meson.build b/backend/src/meson.build new file mode 100644 index 000000000..1668a738e --- /dev/null +++ b/backend/src/meson.build @@ -0,0 +1,8 @@ +example_1 = py.extension_module( + 'example_1', + sources: ['example_1.cpp'], + dependencies: [py_dep, pybind11_dep], + link_language: 'cpp', + install: true, + subdir: package, +) From a587b6535cd18d66a0e6a7bba51d05bc1040a5c1 Mon Sep 17 00:00:00 2001 From: Ankith Date: Sun, 9 Jun 2024 23:44:04 +0530 Subject: [PATCH 03/62] Start integrating dependencies - brigand and cxxopts resolved via meson wraps - blaze, blaze_tensor and sleef are installed via a custom shell script that's called from the meson build --- .gitignore | 3 ++ backend/.gitignore | 10 ++++ backend/buildscripts/build-all.sh | 70 +++++++++++++++++++++++++++ backend/meson.build | 42 ++++++++++++++++ backend/subprojects/blaze.wrap | 8 +++ backend/subprojects/blaze_tensor.wrap | 7 +++ backend/subprojects/brigand.wrap | 8 +++ backend/subprojects/cxxopts.wrap | 14 ++++++ backend/subprojects/sleef.wrap | 8 +++ backend/subprojects/tbb.wrap | 8 +++ backend/subprojects/yaml-cpp.wrap | 8 +++ 11 files changed, 186 insertions(+) create mode 100644 backend/buildscripts/build-all.sh create mode 100644 backend/subprojects/blaze.wrap create mode 100644 backend/subprojects/blaze_tensor.wrap create mode 100644 backend/subprojects/brigand.wrap create mode 100644 backend/subprojects/cxxopts.wrap create mode 100644 backend/subprojects/sleef.wrap create mode 100644 backend/subprojects/tbb.wrap create mode 100644 backend/subprojects/yaml-cpp.wrap diff --git a/.gitignore b/.gitignore index 535b8659c..e799563d2 100644 --- a/.gitignore +++ b/.gitignore @@ -236,3 +236,6 @@ outcmaes/* # csv files *.csv + +# ./backend dependencies +deps diff --git a/backend/.gitignore b/backend/.gitignore index 46f42f8f3..9ac281d62 100644 --- a/backend/.gitignore +++ b/backend/.gitignore @@ -9,3 +9,13 @@ install_manifest.txt compile_commands.json CTestTestfile.cmake _deps + +subprojects/packagecache + +subprojects/blaze +subprojects/blaze_tensor +subprojects/brigand +subprojects/cxxopts-3.1.1 +subprojects/sleef +subprojects/tbb +subprojects/yaml-cpp diff --git a/backend/buildscripts/build-all.sh b/backend/buildscripts/build-all.sh new file mode 100644 index 000000000..7ae68cab9 --- /dev/null +++ b/backend/buildscripts/build-all.sh @@ -0,0 +1,70 @@ + +set -e -x + +OLD_CWD=$(pwd) + +export ELASTICA_DEP_PREFIX="$OLD_CWD/../deps" + +mkdir -p $ELASTICA_DEP_PREFIX && cd $ELASTICA_DEP_PREFIX + +export ELASTICA_INSTALL_PREFIX="$ELASTICA_DEP_PREFIX/installed" + +# With this we +# 1) Force install prefix to $ELASTICA_INSTALL_PREFIX +# 2) use lib directory within $ELASTICA_INSTALL_PREFIX (and not lib64) +# 3) make release binaries +# 4) build shared libraries +# 5) not have @rpath in the linked dylibs (needed on macs only) +# 6) tell cmake to search in $ELASTICA_INSTALL_PREFIX for sub dependencies +export ELASTICA_BASE_CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=$ELASTICA_INSTALL_PREFIX \ + -DCMAKE_INSTALL_LIBDIR:PATH=lib \ + -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_SHARED_LIBS=true \ + -DCMAKE_INSTALL_NAME_DIR=$ELASTICA_INSTALL_PREFIX/lib \ + -DCMAKE_PREFIX_PATH=$ELASTICA_INSTALL_PREFIX" + +mkdir -p $ELASTICA_INSTALL_PREFIX + +# set pkg_config_path so that system can find dependencies +export PKG_CONFIG_PATH="$ELASTICA_INSTALL_PREFIX/lib/pkgconfig:$PKG_CONFIG_PATH" + +# blaze +if [ ! -d blaze ]; then + git clone --depth 1 https://bitbucket.org/blaze-lib/blaze.git && cd blaze + + mkdir build + cmake -B build $ELASTICA_BASE_CMAKE_FLAGS \ + -DBLAZE_SHARED_MEMORY_PARALLELIZATION=OFF \ + -DUSE_LAPACK=OFF \ + -S . + cmake --build build --parallel $(nproc) + cmake --install build + cd .. +fi + +# blaze tensor +if [ ! -d blaze_tensor ]; then + git clone --depth 1 https://github.com/STEllAR-GROUP/blaze_tensor.git && cd blaze_tensor + + mkdir build + cmake -B build $ELASTICA_BASE_CMAKE_FLAGS \ + -Dblaze_DIR="$ELASTICA_INSTALL_PREFIX/share/blaze/cmake/" \ + -S . + cmake --build build --parallel $(nproc) + cmake --install build + cd .. +fi + +# sleef +if [ ! -d sleef ]; then + git clone --depth 1 https://github.com/shibatch/sleef && cd sleef + + mkdir build + cmake -B build $ELASTICA_BASE_CMAKE_FLAGS \ + -S . + cmake --build build --parallel $(nproc) + cmake --install build + cd .. +fi + +cd $OLD_CWD \ No newline at end of file diff --git a/backend/meson.build b/backend/meson.build index 62c509157..e9a1d4650 100644 --- a/backend/meson.build +++ b/backend/meson.build @@ -22,4 +22,46 @@ pybind11_dep = declare_dependency( ).stdout().strip(), ) +fs = import('fs') + +message('Running buildscripts/build-all.sh (might take a while)') +buildscripts = files('buildscripts/build-all.sh') +res = run_command('bash', buildscripts, check: false) +message(res.stdout().strip()) +if res.returncode() != 0 + error(res.stderr().strip()) +endif + +deps_installed = fs.parent(meson.current_source_dir()) / 'deps' / 'installed' + +deps_inc_dirs = [deps_installed / 'include'] +deps_lib_dirs = [deps_installed / 'lib'] + +# see https://github.com/mesonbuild/meson/issues/7943 +# The strategy for the below required dependencies are as follows +# - First, meson tries to resolve the dependency by searching on the system +# - if it doesn't find it, it fallbacks to the subproject wrap system or the +# ones installed in ./buildscripts + +blaze_dep = dependency('blaze') +blaze_tensor_dep = dependency('BlazeTensor') +sleef_dep = dependency('sleef', required: false) +if not sleef_dep.found() + sleef_dep = declare_dependency( + include_directories: deps_inc_dirs, + dependencies: cc.find_library('sleef', dirs: deps_lib_dirs), + ) +endif + +brigand_dep = dependency('brigand', fallback : 'brigand') +cxxopts_dep = dependency('cxxopts', fallback : 'cxxopts') + +# meson-python: error: Could not map installation path to an equivalent wheel directory: +# tbb_dep = dependency('tbb', fallback : 'tbb') +# yaml_cpp_dep = dependency('yaml-cpp', fallback : 'yaml-cpp') + +# TODO: are these required? +# setup_library "HighFive" "https://github.com/BlueBrain/HighFive" +# setup_library "spline" "https://github.com/tp5uiuc/spline.git" + subdir('src') diff --git a/backend/subprojects/blaze.wrap b/backend/subprojects/blaze.wrap new file mode 100644 index 000000000..f7d4cf05f --- /dev/null +++ b/backend/subprojects/blaze.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://bitbucket.org/blaze-lib/blaze.git +revision = head +depth = 1 +method = cmake + +[provide] +blaze = blaze_dep diff --git a/backend/subprojects/blaze_tensor.wrap b/backend/subprojects/blaze_tensor.wrap new file mode 100644 index 000000000..127baae0c --- /dev/null +++ b/backend/subprojects/blaze_tensor.wrap @@ -0,0 +1,7 @@ +[wrap-git] +url = https://github.com/STEllAR-GROUP/blaze_tensor.git +revision = head +method = cmake + +[provide] +BlazeTensor = BlazeTensor_dep diff --git a/backend/subprojects/brigand.wrap b/backend/subprojects/brigand.wrap new file mode 100644 index 000000000..d24c3f31a --- /dev/null +++ b/backend/subprojects/brigand.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/edouarda/brigand.git +revision = head +depth = 1 +method = cmake + +[provide] +brigand = brigand_dep diff --git a/backend/subprojects/cxxopts.wrap b/backend/subprojects/cxxopts.wrap new file mode 100644 index 000000000..d2c9d1b1c --- /dev/null +++ b/backend/subprojects/cxxopts.wrap @@ -0,0 +1,14 @@ +# Taken from meson wrap db +[wrap-file] +directory = cxxopts-3.1.1 +source_url = https://github.com/jarro2783/cxxopts/archive/v3.1.1.tar.gz +source_filename = cxxopts-3.1.1.tar.gz +source_hash = 523175f792eb0ff04f9e653c90746c12655f10cb70f1d5e6d6d9491420298a08 +patch_filename = cxxopts_3.1.1-1_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/cxxopts_3.1.1-1/get_patch +patch_hash = 3a38f72d4a9c394d32db47bbca6f00cb6ee330434b82b653f9d04bdea73170bb +source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/cxxopts_3.1.1-1/cxxopts-3.1.1.tar.gz +wrapdb_version = 3.1.1-1 + +[provide] +cxxopts = cxxopts_dep diff --git a/backend/subprojects/sleef.wrap b/backend/subprojects/sleef.wrap new file mode 100644 index 000000000..6281636c9 --- /dev/null +++ b/backend/subprojects/sleef.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/shibatch/sleef.git +revision = head +depth = 1 +method = cmake + +[provide] +sleef = sleef_dep diff --git a/backend/subprojects/tbb.wrap b/backend/subprojects/tbb.wrap new file mode 100644 index 000000000..f0ee601d7 --- /dev/null +++ b/backend/subprojects/tbb.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/oneapi-src/oneTBB.git +revision = head +depth = 1 +method = cmake + +[provide] +tbb = tbb_dep diff --git a/backend/subprojects/yaml-cpp.wrap b/backend/subprojects/yaml-cpp.wrap new file mode 100644 index 000000000..b1b886c6a --- /dev/null +++ b/backend/subprojects/yaml-cpp.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/jbeder/yaml-cpp.git +revision = head +depth = 1 +method = cmake + +[provide] +yaml-cpp = yaml_cpp_dep From 8c3658477f9df3c93d90403b553394d3e13d9863 Mon Sep 17 00:00:00 2001 From: Ankith Date: Thu, 13 Jun 2024 17:07:22 +0530 Subject: [PATCH 04/62] elasticapp _batch_matmul ports with benchmarking --- backend/benchmarking/matmul.py | 42 +++++++++++++ backend/src/_linalg.cpp | 107 +++++++++++++++++++++++++++++++++ backend/src/meson.build | 9 +++ 3 files changed, 158 insertions(+) create mode 100644 backend/benchmarking/matmul.py create mode 100644 backend/src/_linalg.cpp diff --git a/backend/benchmarking/matmul.py b/backend/benchmarking/matmul.py new file mode 100644 index 000000000..3235e0ef6 --- /dev/null +++ b/backend/benchmarking/matmul.py @@ -0,0 +1,42 @@ +from elasticapp._linalg import batch_matmul_naive, batch_matmul_blaze +from elastica._linalg import _batch_matmul +import numpy +import time + +# warm up jit for fair comparison +random_1 = numpy.random.random((3, 3, 1)) +random_2 = numpy.random.random((3, 3, 1)) +out1 = _batch_matmul(random_1, random_2) + + +def benchmark_batchsize(funcs: list, batches: list[int], num_iterations: int = 1000): + ret: dict = {} + for batch_size in batches: + random_a = numpy.random.random((3, 3, batch_size)) + random_b = numpy.random.random((3, 3, batch_size)) + + ret[batch_size] = {} + for func in funcs: + start = time.perf_counter() + for _ in range(num_iterations): + func(random_a, random_b) + + ret[batch_size][func.__name__] = ( + time.perf_counter() - start + ) / num_iterations + + return ret + + +results = benchmark_batchsize( + [batch_matmul_naive, batch_matmul_blaze, _batch_matmul], [2**i for i in range(14)] +) +for size, data in results.items(): + pyelastica = data["_batch_matmul"] + elasticapp = data["batch_matmul_naive"] + elasticapp_blaze = data["batch_matmul_blaze"] + print(f"{size = }") + print(f"{pyelastica = }") + print(f"{elasticapp = }, ratio: {elasticapp / pyelastica}") + print(f"{elasticapp_blaze = }, ratio: {elasticapp_blaze / pyelastica}") + print() diff --git a/backend/src/_linalg.cpp b/backend/src/_linalg.cpp new file mode 100644 index 000000000..b628ebd2a --- /dev/null +++ b/backend/src/_linalg.cpp @@ -0,0 +1,107 @@ +#include +#include +#include + +namespace py = pybind11; + +using blaze::StaticMatrix; + +/* Simple rewrite of elastica._linalg._batch_matmul */ +py::array_t batch_matmul_naive( + py::array_t first_matrix_collection, + py::array_t second_matrix_collection) +{ + auto s1 = first_matrix_collection.shape(0); + auto s2 = first_matrix_collection.shape(1); + auto max_k = first_matrix_collection.shape(2); + auto output_matrix = py::array_t{{s1, s2, max_k}}; + + auto a = first_matrix_collection.unchecked<3>(); + auto b = second_matrix_collection.unchecked<3>(); + auto c = output_matrix.mutable_unchecked<3>(); + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t m = 0; m < 3; m++) + { + for (py::ssize_t k = 0; k < max_k; k++) + { + c(i, m, k) += a(i, j, k) * b(j, m, k); + } + } + } + } + return output_matrix; +} + +/* blaze implementation of elastica._linalg._batch_matmul */ +py::array_t batch_matmul_blaze( + py::array_t first_matrix_collection, + py::array_t second_matrix_collection) +{ + auto s1 = first_matrix_collection.shape(0); + auto s2 = first_matrix_collection.shape(1); + auto max_k = first_matrix_collection.shape(2); + auto output_matrix = py::array_t{{s1, s2, max_k}}; + + auto a = first_matrix_collection.unchecked<3>(); + auto b = second_matrix_collection.unchecked<3>(); + auto c = output_matrix.mutable_unchecked<3>(); + + StaticMatrix blaze_a; + StaticMatrix blaze_b; + StaticMatrix blaze_c; + for (py::ssize_t k = 0; k < max_k; k++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + blaze_a(i, j) = a(i, j, k); + blaze_b(i, j) = b(i, j, k); + } + } + blaze_c = blaze_a * blaze_b; + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + c(i, j, k) = blaze_c(i, j); + } + } + } + + return output_matrix; +} + +PYBIND11_MODULE(_linalg, m) +{ + m.doc() = R"pbdoc( + elasticapp _linalg + --------------- + + .. currentmodule:: _linalg + + .. autosummary:: + :toctree: _generate + + batch_matmul_naive + )pbdoc"; + + m.def("batch_matmul_naive", &batch_matmul_naive, R"pbdoc( + This is batch matrix matrix multiplication function. Only batch + of 3x3 matrices can be multiplied. + )pbdoc"); + + m.def("batch_matmul_blaze", &batch_matmul_blaze, R"pbdoc( + This is batch matrix matrix multiplication function. Only batch + of 3x3 matrices can be multiplied. + )pbdoc"); + +#ifdef VERSION_INFO + m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); +#else + m.attr("__version__") = "dev"; +#endif +} \ No newline at end of file diff --git a/backend/src/meson.build b/backend/src/meson.build index 1668a738e..9fb245448 100644 --- a/backend/src/meson.build +++ b/backend/src/meson.build @@ -6,3 +6,12 @@ example_1 = py.extension_module( install: true, subdir: package, ) + +_linalg = py.extension_module( + '_linalg', + sources: ['_linalg.cpp'], + dependencies: [py_dep, pybind11_dep, blaze_dep], + link_language: 'cpp', + install: true, + subdir: package, +) From 3e05a9e4bfe4270d37abc409c92723b77b1af2f8 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Fri, 14 Jun 2024 18:14:45 +0900 Subject: [PATCH 05/62] add math utility libraries --- backend/src/CMakeLists.txt | 94 +++++++ backend/src/Utilities/CMakeLists.txt | 109 ++++++++ backend/src/Utilities/DefineTypes.h | 12 + backend/src/Utilities/Math.hpp | 15 ++ .../Math/BlazeDetail/BlazeCalculus.hpp | 91 +++++++ .../Math/BlazeDetail/BlazeLinearAlgebra.hpp | 245 ++++++++++++++++++ .../Math/BlazeDetail/BlazeRotation.hpp | 152 +++++++++++ backend/src/Utilities/Math/CMakeLists.txt | 23 ++ backend/src/Utilities/Math/Rot3.hpp | 21 ++ backend/src/Utilities/Math/Types.hpp | 50 ++++ backend/src/Utilities/Math/Vec3.hpp | 21 ++ 11 files changed, 833 insertions(+) create mode 100644 backend/src/CMakeLists.txt create mode 100644 backend/src/Utilities/CMakeLists.txt create mode 100644 backend/src/Utilities/DefineTypes.h create mode 100644 backend/src/Utilities/Math.hpp create mode 100644 backend/src/Utilities/Math/BlazeDetail/BlazeCalculus.hpp create mode 100644 backend/src/Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp create mode 100644 backend/src/Utilities/Math/BlazeDetail/BlazeRotation.hpp create mode 100644 backend/src/Utilities/Math/CMakeLists.txt create mode 100644 backend/src/Utilities/Math/Rot3.hpp create mode 100644 backend/src/Utilities/Math/Types.hpp create mode 100644 backend/src/Utilities/Math/Vec3.hpp diff --git a/backend/src/CMakeLists.txt b/backend/src/CMakeLists.txt new file mode 100644 index 000000000..03136235a --- /dev/null +++ b/backend/src/CMakeLists.txt @@ -0,0 +1,94 @@ +# Distributed under the MIT License. See LICENSE for details. + +#add_subdirectory(DataStructures) +#add_subdirectory(ErrorHandling) +#add_library(${PROJECT_NAME}::ErrorHandling ALIAS ErrorHandling) + +#add_subdirectory(ModuleSettings) +#add_library(${PROJECT_NAME}::ModuleSettings ALIAS ModuleSettings) + +#add_subdirectory(Options) +#add_library(${PROJECT_NAME}::Options ALIAS Options) + +#add_subdirectory(Parallel) +#add_library(${PROJECT_NAME}::Parallel ALIAS Parallel) + +#add_subdirectory(PythonBindings) +#if (${ELASTICA_BUILD_PYTHON_BINDINGS}) +# add_library(${PROJECT_NAME}::PythonBindings ALIAS PythonBindings) +#endif () + +#add_subdirectory(Simulator) +#add_library(${PROJECT_NAME}::Simulator ALIAS Simulator) + +add_subdirectory(Systems) +add_library(${PROJECT_NAME}::Systems ALIAS Systems) + +#add_subdirectory(Time) +#add_library(${PROJECT_NAME}::Time ALIAS Time) + +add_subdirectory(Utilities) +add_library(${PROJECT_NAME}::Utilities ALIAS Utilities) + +add_library(${PROJECT_NAME} INTERFACE) +target_include_directories(${PROJECT_NAME} INTERFACE + $ + $ + $) + +set(ELASTICA_TOPLEVEL_LIBS + "ErrorHandling;ModuleSettings;Options;Parallel;Simulator;Systems;Time;Utilities") + +foreach (ELASTICA_LIB ${ELASTICA_TOPLEVEL_LIBS}) + target_link_libraries(${PROJECT_NAME} + INTERFACE + ${ELASTICA_LIB} + ) +endforeach (ELASTICA_LIB ${ELASTICA_TOPLEVEL_LIBS}) + +# set to global +set_property(GLOBAL PROPERTY ELASTICA_TOPLEVEL_LIBS_PROPERTY ${ELASTICA_TOPLEVEL_LIBS}) +unset(ELASTICA_TOPLEVEL_LIBS) + +# target_link_libraries(${PROJECT_NAME} INTERFACE libs) + +# In case each library needs to be exported separately +# +# In case each library needs to be exported separately +# foreach (lib ${libs}) +# target_include_directories(${lib} +# $ +# $ +# $) +# endforeach (lib ${libs}) + +#set_property(TARGET ElasticaFlags +# APPEND PROPERTY +# INTERFACE_COMPILE_OPTIONS +# $<$:-stdlib=libc++>) +# +#check_and_add_cxx_link_flag(-stdlib=libc++) +#check_and_add_cxx_link_flag(-libc++experimental) +#check_and_add_cxx_link_flag(-stdlib=libstdc++) + +#target_link_libraries(elastica +# PUBLIC +# Simulator +# Environments +# Utilities +# ErrorHandling +# IO +# Systems +# ElasticaFlags +# ) + +#message(${ElasticaFlags}) +#if ("${CMAKE_CXX_COMPILER_ID}" MATCHES "GNU") +# message("Bitch") +# set(CXX_FILESYSTEM_LIBRARIES "stdc++fs") +# target_link_libraries(elastica PUBLIC ${CXX_FILESYSTEM_LIBRARIES}) +#elseif("${CMAKE_CXX_COMPILER_ID}" MATCHES "Clang") +# set(CXX_FILESYSTEM_LIBRARIES "") +## set(CXX_FILESYSTEM_LIBRARIES "c++experimental") +#endif () +#target_link_libraries(elastica PUBLIC ${CXX_FILESYSTEM_LIBRARIES}) diff --git a/backend/src/Utilities/CMakeLists.txt b/backend/src/Utilities/CMakeLists.txt new file mode 100644 index 000000000..b9bcaad64 --- /dev/null +++ b/backend/src/Utilities/CMakeLists.txt @@ -0,0 +1,109 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY Utilities) + +add_elastica_library(${LIBRARY}) + +#add_subdirectory(Geometry) +#add_subdirectory(Logging) +add_subdirectory(Math) +#add_subdirectory(MemoryMeter) +#add_subdirectory(NamedType) +#add_subdirectory(Singleton) +#add_subdirectory(Timing) +#add_subdirectory(TypeTraits) +#add_subdirectory(ValueTraits) + +# elastica_target_sources( +# ${LIBRARY} +# PRIVATE +# Demangle.cpp +# SystemClock.cpp +# PrettyType.cpp +# WrapText.cpp +# WidthStream.cpp +# ) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + HEADERS + # Approx.hpp + # AsConst.hpp + # Byte.hpp + # Clamp.hpp + # CountingIterator.hpp + # CRTP.hpp + DefineTypes.h + # Demangle.hpp + # Defer.hpp + # End.hpp + # FakeVirtual.hpp + # ForceInline.hpp + # Geometry.hpp + # Get.hpp + # Generators.hpp + # GSL.hpp + # IgnoreUnused.hpp + # Invoke.hpp + # KarmarkarKarpPartition.hpp + # ListRequires.hpp + # Logging.hpp + # MakeArray.hpp + # MakeCopyable.hpp + # MakeFromTuple.hpp + # MakeNamedFunctor.hpp + # MakeSignalingNan.hpp + # MakeString.hpp + Math.hpp + # MemoryMeter.hpp + # MemoryPool.hpp + # NamedTemplate.hpp + # NamedType.hpp + # NumericAt.hpp + # NumericRange.hpp + # NonCopyable.hpp + # NonCreatable.hpp + # NoSuchType.hpp + # Overloader.hpp + # PrettyType.hpp + # PrintHelpers.hpp + # ProtocolHelpers.hpp + # Rampup.hpp + # Registration.hpp + # Requires.hpp + # Singleton.hpp + # StaticConst.hpp + # StaticWarning.hpp + # StdHelpers.hpp + # StdVectorHelpers.hpp + # SystemClock.hpp + # TaggedTuple.hpp + # Thresholds.hpp + # TimeStamps.h + # Timing.h + # TMPL.hpp + # TMPLDebugging.hpp + # TypeTraits.hpp + # Unroll.hpp + # WidthStream.hpp + # WrapText.hpp +) + +target_link_libraries( + ${LIBRARY} + PUBLIC + Blaze + Brigand + ErrorHandling +) + +#set(LIBRARY_SOURCES Utilities.cpp ) +# +#add_elastica_library(${LIBRARY} ${LIBRARY_SOURCES}) +# +#target_link_libraries(${LIBRARY} INTERFACE ErrorHandling) +# + +# add_subdirectory(AppSupport) +# add_subdirectory(ConvertCase) diff --git a/backend/src/Utilities/DefineTypes.h b/backend/src/Utilities/DefineTypes.h new file mode 100644 index 000000000..bd9a7555e --- /dev/null +++ b/backend/src/Utilities/DefineTypes.h @@ -0,0 +1,12 @@ +#pragma once +#include + +namespace elastica { + + using real_t = double; +// typedef REAL real; +#define ALIGNMENT 32 + + using id_t = uint64_t; + +} // namespace elastica diff --git a/backend/src/Utilities/Math.hpp b/backend/src/Utilities/Math.hpp new file mode 100644 index 000000000..698ba17c5 --- /dev/null +++ b/backend/src/Utilities/Math.hpp @@ -0,0 +1,15 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +// Free operations +// #include "Math/Invert.hpp" +// #include "Math/Square.hpp" +// #include "Math/SqrLength.hpp" + +// Containers +#include "Math/Vec3.hpp" + +//#include diff --git a/backend/src/Utilities/Math/BlazeDetail/BlazeCalculus.hpp b/backend/src/Utilities/Math/BlazeDetail/BlazeCalculus.hpp new file mode 100644 index 000000000..07cc5404b --- /dev/null +++ b/backend/src/Utilities/Math/BlazeDetail/BlazeCalculus.hpp @@ -0,0 +1,91 @@ +#pragma once + +#include "blaze/Blaze.h" +#include "blaze_tensor/Blaze.h" + +namespace elastica{ + //************************************************************************** + /*!\brief Simple 2-point difference rule with zero at end points. + // + // \details + // Discrete 2-point difference in \elastica of a function f:[a,b]-> R, i.e + // D f[a,b] -> df[a,b] where f satisfies the conditions + // f(a) = f(b) = 0.0. Operates from rod's elemental space to nodal space. + // + // \example + // The following shows a typical use of the difference kernel + // with the expected (correct) result also shown. + // \snippet test_calculus.cpp difference_kernel_example + // + // \param[out] out_matrix(3, n_nodes) difference values + // \param[in] in_matrix(3, n_elems) vector batch \n + // where n_nodes = n_elems + 1 + // + // \return void/None + // + // \see fill later? + */ + template // blaze Matrix expression type 2 + void two_point_difference_kernel(MT1& out_matrix, const MT2& in_matrix) { + constexpr std::size_t dimension(3UL); + assert(in_matrix.rows() == dimension); + assert(out_matrix.rows() == dimension); + const std::size_t n_elems = in_matrix.columns(); + const std::size_t n_nodes = n_elems + 1UL; + //assert(out_matrix.columns() == n_nodes); + + blaze::column(out_matrix, 0UL) = blaze::column(in_matrix, 0UL); + blaze::column(out_matrix, n_nodes - 1UL) = + -blaze::column(in_matrix, n_elems - 1UL); + blaze::submatrix(out_matrix, 0UL, 1UL, dimension, n_elems - 1UL) = + blaze::submatrix(in_matrix, 0UL, 1UL, dimension, n_elems - 1UL) - + blaze::submatrix(in_matrix, 0UL, 0UL, dimension, n_elems - 1UL); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Simple trapezoidal quadrature rule with zero at end points. + // + // \details + // Discrete integral of a function in \elastica + // \f$ : [a,b] \rightarrow \mathbf{R}, \int_{a}^{b} f \rightarrow \mathbf{R} \f$ + // where f satisfies the conditions f(a) = f(b) = 0.0. + // Operates from rod's elemental space to nodal space. + // + // \example + // The following shows a typical use of the quadrature kernel + // with the expected (correct) result also shown. + // \snippet test_calculus.cpp quadrature_kernel_example + // + // \param[out] out_matrix(3, n_nodes) quadrature values + // \param[in] in_matrix(3, n_elems) vector batch \n + // where n_nodes = n_elems + 1 + // + // \return void/None + // + // \see fill later? + */ + template // blaze Matrix expression type 2 + void quadrature_kernel(MT1& out_matrix, const MT2& in_matrix) { + constexpr std::size_t dimension(3UL); + const std::size_t n_elems = in_matrix.columns(); + assert(in_matrix.rows() == dimension); + assert(out_matrix.rows() == dimension); + const std::size_t n_nodes = n_elems + 1UL; + assert(out_matrix.columns() == n_nodes); + using ValueType = typename MT1::ElementType; + + blaze::column(out_matrix, 0UL) = + ValueType(0.5) * blaze::column(in_matrix, 0UL); + blaze::column(out_matrix, n_nodes - 1UL) = + ValueType(0.5) * blaze::column(in_matrix, n_elems - 1UL); + blaze::submatrix(out_matrix, 0UL, 1UL, dimension, n_elems - 1UL) = + ValueType(0.5) * + (blaze::submatrix(in_matrix, 0UL, 1UL, dimension, n_elems - 1UL) + + blaze::submatrix(in_matrix, 0UL, 0UL, dimension, n_elems - 1UL)); + } + //************************************************************************** + +} // namespace elastica diff --git a/backend/src/Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp b/backend/src/Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp new file mode 100644 index 000000000..1c3060573 --- /dev/null +++ b/backend/src/Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp @@ -0,0 +1,245 @@ +#pragma once + +#include "blaze/Blaze.h" +#include "blaze_tensor/Blaze.h" + +/* some of them have for loops until blaze_tensor comes up with 3D tensor */ +/* products imp */ +/* things would have been much simpler if there was einsum() in blaze */ +/* NOTE: % operator corresponds to Schur product for matrices and tensors */ +namespace elastica { + //************************************************************************** + /*!\brief Vector Diference + // + // \param[out] output_vector(3, n_elems) + // \param[in] input_vector(3, n_elems) + // + // \return void/None + // + // \see TODO: fill later + */ + template + inline auto difference_kernel(const V& in_vector){ + constexpr std::size_t dimension(3UL); + assert(in_vector.rows() == dimension); + const std::size_t n_nodes = in_vector.columns(); + const std::size_t n_elems = n_nodes - 1UL; + + return blaze::submatrix(in_vector, 0UL, 1UL, dimension, n_elems) - + blaze::submatrix(in_vector, 0UL, 0UL, dimension, n_elems); + } + + //************************************************************************** + /*!\brief Batchwise matrix-vector product. + // + // \details + // Computes a batchwise matrix-vector product given in indical notation: \n + // matvec_batch{ik} = matrix_batch{ijk} * vector_batch{jk} + // + // \example + // The following shows a typical use of the batch_matvec function + // with the expected (correct) result also shown. + // \snippet test_linalg.cpp batch_matvec_example + // + // \param[out] matvec_batch(3, n_elems) matrix-vector product + // \param[in] matrix_batch(3, 3, n_elems) + // \param[in] vector_batch(3, n_elems) + // + // \return void/None + // + // \see fill later? + */ + template // blaze Matrix expression type 2 + void batch_matvec(MT1& matvec_batch, + const TT& matrix_batch, + const MT2& vector_batch) { + constexpr std::size_t dimension(3UL); + //const std::size_t n_elems = matrix_batch.columns(); + + // Written for debugging purpose +// const std::size_t vbatch_columns = vector_batch.columns(); +// const std::size_t vbatch_rows = vector_batch.rows(); +// const std::size_t mbatch_pages = matrix_batch.pages(); +// const std::size_t mbatch_columns = matrix_batch.columns(); +// const std::size_t mbatch_rows = matrix_batch.rows(); +// +// const std::size_t size1 = blaze::pageslice(matrix_batch, 0UL).columns(); +// const std::size_t size0 = blaze::pageslice(matrix_batch, 0UL).rows(); + // assert(matvec_batch.columns() == n_elems); + // assert(matvec_batch.rows() == dimension); + // assert(matrix_batch.pages() == dimension); + // assert(matrix_batch.rows() == dimension); + // assert(vector_batch.columns() == n_elems); + // assert(vector_batch.rows() == dimension); + + for (std::size_t i(0UL); i < dimension; ++i) { + auto val1 = blaze::pageslice(matrix_batch, i); + auto val2 = blaze::sum(val1 % vector_batch); + blaze::row(matvec_batch, i) = val2; + } + } + //************************************************************************** + + //************************************************************************** + /*!\brief Batchwise matrix-matrix product. + // + // \details + // Computes a batchwise matrix-matrix product given in + // indical notation: \n + // matmul_batch{ilk} = first_matrix_batch{ijk} * second_matrix_batch{jlk} + // + // \example + // The following shows a typical use of the batch_matmul function + // with the expected (correct) result also shown. + // \snippet test_linalg.cpp batch_matmul_example + // + // \param[out] matmul_batch(3, 3, n_elems) matrix-matrix product + // \param[in] first_matrix_batch(3, 3, n_elems) + // \param[in] second_matrix_batch(3, 3, n_elems) + // + // \return void/None + // + // \see fill later? + */ + template // blaze Tensor expression type 3 + void batch_matmul(TT1& matmul_batch, + const TT2& first_matrix_batch, + const TT3& second_matrix_batch) { + constexpr std::size_t dimension(3UL); + const std::size_t n_elems = first_matrix_batch.columns(); + assert(matmul_batch.pages() == dimension); + assert(matmul_batch.rows() == dimension); + assert(matmul_batch.columns() == n_elems); + assert(first_matrix_batch.pages() == dimension); + assert(first_matrix_batch.rows() == dimension); + assert(second_matrix_batch.pages() == dimension); + assert(second_matrix_batch.rows() == dimension); + assert(second_matrix_batch.columns() == n_elems); + + for (std::size_t i(0UL); i < dimension; ++i) + for (std::size_t j(0UL); j < dimension; ++j) + /* loop over dimensions, lesser iterations, bigger slices */ + blaze::row(blaze::pageslice(matmul_batch, i), j) = + blaze::sum( + blaze::pageslice(first_matrix_batch, i) % + blaze::trans(blaze::rowslice(second_matrix_batch, j))); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Batchwise vector-vector cross product. + // + // \details + // Computes a batchwise vector-vector cross product given in + // indical notation: \n + // cross_batch{il} = LCT{ijk} * first_vector_batch{jl} * + // second_vector_batch{kl} \n + // where LCT is the Levi-Civita Tensor + // + // \example + // The following shows a typical use of the batch_cross function + // with the expected (correct) result also shown. + // \snippet test_linalg.cpp batch_cross_example + // + // \param[out] cross_batch(3, n_elems) vector-vector cross product + // \param[in] first_vector_batch(3, n_elems) + // \param[in] second_vector_batch(3, n_elems) + // + // \return void/None + // + // \see fill later? + */ + template // blaze Matrix expression type 3 + void batch_cross(MT1& cross_batch, + const MT2& first_vector_batch, + const MT3& second_vector_batch) { + constexpr std::size_t dimension(3UL); + //const std::size_t n_elems = first_vector_batch.columns(); + assert(cross_batch.rows() == dimension); + //assert(cross_batch.columns() == n_elems); + assert(first_vector_batch.rows() == dimension); + assert(second_vector_batch.rows() == dimension); + //assert(second_vector_batch.columns() == n_elems); + + for (std::size_t i(0UL); i < dimension; ++i) + /* loop over dimensions, lesser iterations, bigger slices */ + /* remainder operator % cycles the indices across dimension*/ + blaze::row(cross_batch, i) = + blaze::row(first_vector_batch, (i + 1UL) % dimension) * + blaze::row(second_vector_batch, (i + 2UL) % dimension) - + blaze::row(first_vector_batch, (i + 2UL) % dimension) * + blaze::row(second_vector_batch, (i + 1UL) % dimension); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Batchwise vector-vector dot product. + // + // \details + // Computes a batchwise vector-vector dot product given in indical + // notation: \n + // dot_batch{j} = first_vector_batch{ij} * second_vector_batch{ij} + // + // \example + // The following shows a typical use of the batch_dot function + // with the expected (correct) result also shown. + // \snippet test_linalg.cpp batch_dot_example + // + // \param[in] first_vector_batch(3, n_elems) + // \param[in] second_vector_batch(3, n_elems) + // + // \return dot_batch(n_elems): ET object pointer to the dot + // product (for lazy evaluation). + // + // \see fill later? + */ + template // blaze Matrix expression type 2 + auto batch_dot(const MT1& first_vector_batch, + const MT2& second_vector_batch) { + constexpr std::size_t dimension(3UL); + const std::size_t n_elems = first_vector_batch.columns(); + assert(first_vector_batch.rows() == dimension); + assert(second_vector_batch.rows() == dimension); + assert(second_vector_batch.columns() == n_elems); + + return blaze::trans( + blaze::sum(first_vector_batch % second_vector_batch)); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Batchwise vector L2 norm. + // + // \details + // Computes a batchwise vector L2 norm given in indical notation: \n + // norm_batch{j} = (vector_batch{ij} * vector_batch{ij})^0.5 + // + // \example + // The following shows a typical use of the batch_norm function + // with the expected (correct) result also shown. + // \snippet test_linalg.cpp batch_norm_example + // + // \param[in] vector_batch(3, n_elems) + // + // \return norm_batch(n_elems): ET object pointer to the vector L2 norm + // (for lazy evaluation). + // + // \see fill later? + */ + template // blaze vector expression type + auto batch_norm(const MT& vector_batch) { + constexpr std::size_t dimension(3UL); + //const std::size_t n_elems = vector_batch.columns(); + assert(vector_batch.rows() == dimension); + + return blaze::sqrt(blaze::trans(blaze::sum(vector_batch % vector_batch))); + } + +} // namespace rod diff --git a/backend/src/Utilities/Math/BlazeDetail/BlazeRotation.hpp b/backend/src/Utilities/Math/BlazeDetail/BlazeRotation.hpp new file mode 100644 index 000000000..26eef583a --- /dev/null +++ b/backend/src/Utilities/Math/BlazeDetail/BlazeRotation.hpp @@ -0,0 +1,152 @@ +#pragma once + +#include "blaze/Blaze.h" +#include "blaze_tensor/Blaze.h" + +namespace elastica { + // LOG AND EXP DEFINED AS PER ELASTICA'S ROTATIONAL DIRECTIONS! + + //************************************************************************** + /*!\brief Batchwise matrix logarithmic operator. + // + // \details + // Batchwise for rotation matrix R computes the corresponding rotation + // axis vector {theta (u)} using the matrix log() operator: \n + // if theta == 0: + // theta (u) = 0 \n + // else: + // theta (u) = -theta * inv_skew[R - transpose(R)] / sin(theta) \n + // where theta = acos(0.5 * (trace(R) - 1)) and inv_skew[] corresponds + // to an inverse skew symmetric mapping from a skew symmetric matrix M + // to vector V as: + // 
+       |0 -z y|        |x|
+   M = |z 0 -x| to V = |y|
+       |-y x 0|        |z|
+
+ // + // \example + // The following shows a typical use of the log_batch function + // with the expected (correct) result also shown. + // \snippet test_rotations.cpp log_batch_example + // + // \param[out] void/None + // \param[in] rot_matrix_batch(3, 3, n_elems) rotation matrix batch + // + // \return rot_axis_vector_batch(3, n_elems) rotation axis vector batch + // + // \see fill later? + */ + template // blaze Tensor expression type + void batch_inv_rotate(MT& rot_axis_vector_batch, const TT& rot_matrix_batch) { + constexpr std::size_t dimension(3UL); + const std::size_t n_elems = rot_matrix_batch.columns(); + using ValueType = typename TT::ElementType; + assert(rot_matrix_batch.pages() == dimension); + assert(rot_matrix_batch.rows() == dimension); + assert(rot_axis_vector_batch.rows() == dimension); + assert(rot_axis_vector_batch.columns() == n_elems); + + /* hardcoded to avoid external memory passing */ + /* also now operates only on views */ + auto rot_matrix_trace_batch = + blaze::row(blaze::pageslice(rot_matrix_batch, 0UL), 0UL) + + blaze::row(blaze::pageslice(rot_matrix_batch, 1UL), 1UL) + + blaze::row(blaze::pageslice(rot_matrix_batch, 2UL), 2UL); + /* theta = acos((tr(R) - 1) / 2) */ + auto theta_batch = blaze::acos(ValueType(0.5) * + (rot_matrix_trace_batch - ValueType(1.0) - + ValueType(1e-12))); // TODO refactor 1e-12 + auto R_minus_RT = + rot_matrix_batch - blaze::trans(rot_matrix_batch, {1, 0, 2}); // TODO Something is wrong + /* theta (u) = -theta * inv_skew([R - RT]) / 2 sin(theta) */ + for (std::size_t i(0UL); i < dimension; ++i) { + auto inv_skew_symmetric_map = blaze::row( + blaze::pageslice(R_minus_RT, (i + 2UL) % dimension), + (i + 1UL) % dimension); // % for the cyclic rotation of indices + blaze::row(rot_axis_vector_batch, i) = + ValueType(-0.5) * theta_batch * inv_skew_symmetric_map / + (blaze::sin(theta_batch) + + ValueType(1e-14)); // TODO refactor 1e-14 + } + } + //************************************************************************** + + //************************************************************************** + /*!\brief Batchwise matrix exponential operator. + // + // \details + // Batchwise for rotation axis vector {theta u} computes the corresponding + // rotation matrix R using the matrix exp() operator (Rodrigues formula): \n + // R = I - sin(theta) * U + (1 - cos(theta)) * U * U \n + // Here a different tensorial form is implemented as follows: + // if i == j: + // R{ij} = cos(theta) + (1 - cos(theta)) * (u{i})^2 \n + // else: + // R{ij} = (1 - cos(theta)) * u{i} * u{j} - sin(theta) * + // skew_sym[u]{ij} \n + // where I is the identity matrix and skew_sym[] corresponds to a skew + // symmetric mapping from a vector V to a skew symmetric matrix M as: + // 
+       |x|        |0 -z y|
+   V = |y| to M = |z 0 -x|
+       |z|        |-y x 0|
+
+ // + // \example + // The following shows a typical use of the exp_batch function + // with the expected (correct) result also shown. + // \snippet test_rotations.cpp exp_batch_example + // + // \param[out] rot_matrix_batch(3, 3, n_elems) rotation matrix batch + // \param[in] rot_axis_vector_batch(3, n_elems) rotation axis vector batch + // + // \return void/None + // + // \see fill later? + */ + template // blaze Matrix expression type + void exp_batch(TT& rot_matrix_batch, const MT& rot_axis_vector_batch) { + constexpr std::size_t dimension(3UL); + //const std::size_t n_elems = rot_axis_vector_batch.columns(); + using ValueType = typename MT::ElementType; + assert(rot_axis_vector_batch.rows() == dimension); + assert(rot_matrix_batch.pages() == dimension); + assert(rot_matrix_batch.rows() == dimension); + //assert(rot_matrix_batch.columns() == n_elems); + + auto theta_batch = blaze::sqrt( + blaze::sum(blaze::pow(rot_axis_vector_batch, 2))); + /* TODO refactor 1e-14 */ + auto unit_rot_axis_vector_batch = + rot_axis_vector_batch % + blaze::expand(blaze::pow(theta_batch + ValueType(1e-14), -1), + dimension); + for (std::size_t i(0UL); i < dimension; ++i) { + // if i == j: + // R{ij} = cos(theta) + (1 - cos(theta)) * (u{i})^2 + blaze::row(blaze::pageslice(rot_matrix_batch, i), i) = + blaze::cos(theta_batch) + + (ValueType(1.0) - blaze::cos(theta_batch)) * + blaze::pow(blaze::row(unit_rot_axis_vector_batch, i), 2); + // else: + // R{ij} = (1 - cos(theta)) * u{i} * u{j} - sin(theta) * + // skew_sym[u]{ij} + for (std::size_t j : {i + 1UL, i + 2UL}) { + auto skew_symmetric_map = + std::pow(-1, + j - i) * // Sign-bit to check order of entries + blaze::row(unit_rot_axis_vector_batch, + dimension - i - (j % dimension)); // % source index; + blaze::row(blaze::pageslice(rot_matrix_batch, i), j % dimension) = + (ValueType(1.0) - blaze::cos(theta_batch)) * + blaze::row(unit_rot_axis_vector_batch, i) * + blaze::row(unit_rot_axis_vector_batch, j % dimension) - + blaze::sin(theta_batch) * skew_symmetric_map; + } + } + } + //************************************************************************** + +} // namespace elastica diff --git a/backend/src/Utilities/Math/CMakeLists.txt b/backend/src/Utilities/Math/CMakeLists.txt new file mode 100644 index 000000000..d3a595afa --- /dev/null +++ b/backend/src/Utilities/Math/CMakeLists.txt @@ -0,0 +1,23 @@ +set(LIBRARY Utilities) + +elastica_target_sources( + ${LIBRARY} + PRIVATE +) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + HEADERS + Rot3.hpp + Vec3.hpp + Types.hpp +) + +add_subdirectory(Python) + +target_link_libraries( + ${LIBRARY} + PRIVATE + Blaze +) diff --git a/backend/src/Utilities/Math/Rot3.hpp b/backend/src/Utilities/Math/Rot3.hpp new file mode 100644 index 000000000..c1e840b5d --- /dev/null +++ b/backend/src/Utilities/Math/Rot3.hpp @@ -0,0 +1,21 @@ +//============================================================================== +/*! +// \file +// \brief Header for a 3D matrix for use in \elastica interface +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +// +#include "Utilities/Math/Types.hpp" diff --git a/backend/src/Utilities/Math/Types.hpp b/backend/src/Utilities/Math/Types.hpp new file mode 100644 index 000000000..6c84669e7 --- /dev/null +++ b/backend/src/Utilities/Math/Types.hpp @@ -0,0 +1,50 @@ +//============================================================================== +/*! +// \file +// \brief Types belonging to the Math module +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +// +#include "Utilities/DefineTypes.h" + +namespace elastica { + + //**************************************************************************** + //! \brief 3D vector for use in \elastica interfaces + //! \ingroup math + using Vec3 = blaze::StaticVector; + //**************************************************************************** + + //**************************************************************************** + //! \brief 3D matrix for use in \elastica interfaces + //! \ingroup math + using Rot3 = blaze::StaticMatrix; + //**************************************************************************** + + ////////////////////////////////////////////////////////////////////////////// + // + // Forward declarations + // + ////////////////////////////////////////////////////////////////////////////// + + //**************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + struct RotationConvention; + /*! \endcond */ + //**************************************************************************** +} // namespace elastica diff --git a/backend/src/Utilities/Math/Vec3.hpp b/backend/src/Utilities/Math/Vec3.hpp new file mode 100644 index 000000000..f73a7d5be --- /dev/null +++ b/backend/src/Utilities/Math/Vec3.hpp @@ -0,0 +1,21 @@ +//============================================================================== +/*! +// \file +// \brief Header for a 3D vector for use in \elastica interface +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +// +#include "Utilities/Math/Types.hpp" From b3221cd2c49e24014db8ec59d051c0a63a2726c3 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Fri, 14 Jun 2024 18:46:04 +0900 Subject: [PATCH 06/62] add blaze traits --- backend/src/Systems/CMakeLists.txt | 62 ++++ .../src/Systems/CosseratRods/CMakeLists.txt | 65 ++++ .../CosseratRods/Traits/CMakeLists.txt | 41 +++ .../CosseratRods/Traits/DataOpsTraits.hpp | 87 +++++ .../DataType/BlazeBackend/DataTypeTraits.hpp | 21 ++ .../DataType/BlazeBackend/MatrixTag.hpp | 281 ++++++++++++++ .../DataType/BlazeBackend/TensorTag.hpp | 344 ++++++++++++++++++ .../DataType/BlazeBackend/VectorTag.hpp | 274 ++++++++++++++ .../Traits/DataType/CMakeLists.txt | 40 ++ .../Traits/DataType/Protocols.hpp | 197 ++++++++++ .../CosseratRods/Traits/DataType/Rank.hpp | 40 ++ .../Traits/DataType/ScalarTag.hpp | 235 ++++++++++++ .../src/Systems/CosseratRods/Traits/Types.hpp | 35 ++ 13 files changed, 1722 insertions(+) create mode 100644 backend/src/Systems/CMakeLists.txt create mode 100644 backend/src/Systems/CosseratRods/CMakeLists.txt create mode 100644 backend/src/Systems/CosseratRods/Traits/CMakeLists.txt create mode 100644 backend/src/Systems/CosseratRods/Traits/DataOpsTraits.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/DataTypeTraits.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/MatrixTag.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/TensorTag.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/VectorTag.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/CMakeLists.txt create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/Protocols.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/Rank.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/DataType/ScalarTag.hpp create mode 100644 backend/src/Systems/CosseratRods/Traits/Types.hpp diff --git a/backend/src/Systems/CMakeLists.txt b/backend/src/Systems/CMakeLists.txt new file mode 100644 index 000000000..59bcc17ed --- /dev/null +++ b/backend/src/Systems/CMakeLists.txt @@ -0,0 +1,62 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY Systems) + +add_elastica_library(${LIBRARY} INTERFACE) + +# add_subdirectory(Block) +# add_library(${LIBRARY}::Block ALIAS Block) + +# add_subdirectory(common) +# add_library(${LIBRARY}::Common ALIAS SystemsCommon) + +add_subdirectory(CosseratRods) +add_library(${LIBRARY}::CosseratRods ALIAS CosseratRods) + +# add_subdirectory(Customization) +# add_library(${LIBRARY}::Customization ALIAS SystemsCustomization) + +# add_subdirectory(Python) + +# add_subdirectory(RigidBody) +# add_library(${LIBRARY}::RigidBody ALIAS RigidBody) + +# add_subdirectory(States) +# add_library(${LIBRARY}::States ALIAS States) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + # Block.hpp + CosseratRods.hpp + # Protocols.hpp + # RigidBody.hpp + # Sphere.hpp + # Systems.hpp + # Tags.hpp + # Types.hpp +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + # Block.hpp + CosseratRods.hpp + # Protocols.hpp + # RigidBody.hpp + # Sphere.hpp + # Systems.hpp + # Tags.hpp + # Types.hpp +) + +target_link_libraries( + ${LIBRARY} + INTERFACE + # Block + # SystemsCommon + CosseratRods + # SystemsCustomization + # RigidBody + # States +) diff --git a/backend/src/Systems/CosseratRods/CMakeLists.txt b/backend/src/Systems/CosseratRods/CMakeLists.txt new file mode 100644 index 000000000..78cd4c5b1 --- /dev/null +++ b/backend/src/Systems/CosseratRods/CMakeLists.txt @@ -0,0 +1,65 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY CosseratRods) +add_elastica_library(${LIBRARY} INTERFACE) + +# set(LIBRARY_SOURCES +# Access.hpp +# Aliases.hpp +# Block.hpp +# BlockInitializer.hpp +# BlockSlice.hpp +# BlockView.hpp +# Components.hpp +# CosseratRodPlugin.hpp +# CosseratRods.hpp +# CosseratRodTraits.hpp +# Initializers.hpp +# Protocols.hpp +# Serialize.hpp # TODO : refactor +# Specialize.hpp +# Tags.hpp +# Types.hpp +# TypeTraits.hpp +# Utility.hpp +# ) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + ${LIBRARY_SOURCES} +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + ${LIBRARY_SOURCES} +) + +#add_subdirectory(Aliases) +#add_library(${LIBRARY}::Aliases ALIAS CosseratRodAliases) + +#add_subdirectory(Components) +#add_library(${LIBRARY}::Components ALIAS CosseratRodComponents) + +#add_subdirectory(Initializers) +#add_library(${LIBRARY}::Initializers ALIAS CosseratRodInitializers) + +add_subdirectory(Traits) +add_library(${LIBRARY}::Traits ALIAS CosseratRodTraits) + +#add_subdirectory(Utility) +#add_library(${LIBRARY}::Utility ALIAS CosseratRodUtility) + +target_link_libraries( + ${LIBRARY} + INTERFACE + CosseratRodAliases + CosseratRodComponents + CosseratRodInitializers + CosseratRodTraits + CosseratRodUtility + Utilities +) + +add_subdirectory(Python) diff --git a/backend/src/Systems/CosseratRods/Traits/CMakeLists.txt b/backend/src/Systems/CosseratRods/Traits/CMakeLists.txt new file mode 100644 index 000000000..8aad7e44c --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/CMakeLists.txt @@ -0,0 +1,41 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY CosseratRodTraits) + +add_elastica_library(${LIBRARY} INTERFACE) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + DataOpsTraits.hpp + # PlacementTrait.hpp + Protocols.hpp + Types.hpp +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + DataOpsTraits.hpp + # PlacementTrait.hpp + Protocols.hpp + Types.hpp +) + +add_subdirectory(DataType) +add_library(${LIBRARY}::DataTraits ALIAS CosseratRodDataTraits) + +# add_subdirectory(PlacementTraits) +# add_library(${LIBRARY}::PlacementTraits ALIAS CosseratRodPlacementTraits) + +add_subdirectory(Operations) +add_library(${LIBRARY}::OperationTraits ALIAS CosseratRodOperationTraits) + +target_link_libraries( + ${LIBRARY} + INTERFACE + CosseratRodDataTraits + # CosseratRodPlacementTraits + CosseratRodOperationTraits +) + diff --git a/backend/src/Systems/CosseratRods/Traits/DataOpsTraits.hpp b/backend/src/Systems/CosseratRods/Traits/DataOpsTraits.hpp new file mode 100644 index 000000000..cb0d0e816 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataOpsTraits.hpp @@ -0,0 +1,87 @@ +//============================================================================== +/*! +// \file +// \brief Datatraits for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +/// +#include "Systems/CosseratRods/Traits/Types.hpp" +#include "Utilities/DefineTypes.h" +/// +//#include "Systems/CosseratRods/Traits/DataTraits/Rank.hpp" +// Not used, since we directly use blaze's implementation of a vector +//#include "Systems/CosseratRods/Traits/DataType/ScalarTag.hpp" +/// + +#define ELASTICA_USE_BLAZE 1 // todo : from CMAKE + +#if defined(ELASTICA_USE_BLAZE) +#include "DataType/BlazeBackend/DataTypeTraits.hpp" +#include "Operations/BlazeBackend/OpsTraits.hpp" +#include "Operations/BlazeBackend/OptimizationLevel.hpp" +#endif // ELASTICA_USE_BLAZE + +/// +#include "Utilities/NonCreatable.hpp" + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Traits implementing data-structures for Cosserat rods in \elastica + * \ingroup cosserat_rod_traits + * + * \details + * DataOpsTraits is the customization point for altering the data-structures + * to be used within the Cosserat rod hierarchy implemented using @ref + * blocks in \elastica. It defines (domain-specific) types corresponding to + * a Cosserat rod (such as a matrix, vector etc.), and is intended for use + * as a template parameter in CosseratRodTraits. + * + * \see elastica::cosserat_rod::CosseratRodTraits + */ + struct DataOpsTraits : private NonCreatable { + //**Type definitions****************************************************** + //! Real type + using real_type = real_t; + //! Type of index + using index_type = std::size_t; + // using OnRod = ScalarTag; + //! Type to track indices on a Cosserat rod, with shape (1, ) + using Index = VectorTag; + //! Type to track scalars on a Cosserat rod, with shape (1, ) + using Scalar = VectorTag; + //! Type to track vectors on a Cosserat rod, with shape (d) + using Vector = MatrixTag; + //! Type to track matrices on a Cosserat rod, with shape (d, d) + using Matrix = TensorTag; + //************************************************************************ + + //**Operation definitions************************************************* + using Operations = OpsTraits; + //************************************************************************ + }; + //************************************************************************** + + } // namespace cosserat_rod +} // namespace elastica \ No newline at end of file diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/DataTypeTraits.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/DataTypeTraits.hpp new file mode 100644 index 000000000..2681724a6 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/DataTypeTraits.hpp @@ -0,0 +1,21 @@ +//============================================================================== +/*! +// \file +// \brief Blaze data-traits for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include "Systems/CosseratRods/Traits/DataType/BlazeBackend/MatrixTag.hpp" +#include "Systems/CosseratRods/Traits/DataType/BlazeBackend/TensorTag.hpp" +#include "Systems/CosseratRods/Traits/DataType/BlazeBackend/VectorTag.hpp" diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/MatrixTag.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/MatrixTag.hpp new file mode 100644 index 000000000..f91f2a146 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/MatrixTag.hpp @@ -0,0 +1,281 @@ +//============================================================================== +/*! +// \file +// \brief Tag indicating a matrix (Rank 1) for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include "Configuration/Systems.hpp" // For index checking +// assert +#include "ErrorHandling/Assert.hpp" +#include "Systems/CosseratRods/Traits/DataType/Protocols.hpp" +#include "Systems/CosseratRods/Traits/DataType/Rank.hpp" + +// matrix types +#include +// #include +#include +#include +#include + +// this comes from the simulator module, which seems like an anti-pattern +// should be in Systems instead +#include "Simulator/Frames.hpp" +#include "Utilities/ProtocolHelpers.hpp" + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Tag for a Rank 2 tensor within Cosserat rods in \elastica + * \ingroup cosserat_rod_traits + * + * \details + * MatrixTag is used within the Cosserat rod components of \elastica to + * indicate a vector. + * + * \tparam Real a floating point type + */ + template + struct MatrixTag : public ::tt::ConformsTo { + //**Type definitions****************************************************** + //! The main type for a TaggedTuple + using type = + blaze::DynamicMatrix /*, TagType*/>; + //! Data type for better readability + using data_type = type; + //! The type of slice for the data type + using slice_type = blaze::Submatrix_; + //! The type of const slice for the data type + using const_slice_type = + const blaze::Submatrix_; + //! The type of a ghost element for the data type + using ghost_type = blaze::StaticVector; + //! The type of a reference to the underlying data type + using reference_type = blaze::Column_; + //! The type of const reference to the underlying data type + using const_reference_type = const blaze::Column_; + //! The type of reference to the slice type + using reference_slice_type = + blaze::Subvector_, blaze::unaligned>; + //! The type of const reference to the slice type + using const_reference_slice_type = + const blaze::Subvector_, + blaze::unaligned>; + //! The rank of the data type + using rank = Rank<2U>; + //************************************************************************ + + //**Utility functions***************************************************** + /*!\name Utility functions */ + //@{ + + //************************************************************************ + /*!\brief Obtain the ghost value + * + * \details + * Obtains a default value for putting in ghost elements. This function + * can be overriden in case the user wants to. This should not be a + * Real (because implicit conversion fills entire data-structure with + * that particular value) but rather a ghost_type object (like a matrix + * for example) + */ + static inline constexpr auto ghost_value() noexcept -> ghost_type { + constexpr Real zero(0.0); + return ghost_type{zero, zero, zero}; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is typically used in the block to generate BlockSlices or access + * indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(data_type& data, std::size_t index) { + return blaze::column(data, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type at a given index + * + * \details + * Overload to obtain a slice of const data at `index`. + * It is typically used in the block to generate ConstBlockSlices or + * access indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(data_type const& data, + std::size_t index) { + return blaze::column(data, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain subslice of the slice type at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * It is typically used in the block to generate slices that are then + * filled in by ghost values. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type& slice, std::size_t index) { + return blaze::column(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain sub-slice of the const slice type at a given index + * + * \details + * Overload to obtain a subslice of const slice at `index`. + * + * It is sometimes needed in a Plugin where the context is const, but the + * data of the underlying structure may well be modified. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type const& slice, + std::size_t index) { + return blaze::column(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain subslice of a temporary slice type at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * + * It is typically used when lazily generating individual slices from a + * lazily generated slice : for example blockslice.get_position(i), where + * get_position() lazily generates a slice, which then gets further + * sliced. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type&& slice, + std::size_t index) { + return blaze::column(static_cast(slice), index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain sub-slice of the const_slice type at a given index + * + * \details + * Overload to obtain a subslice of const_slice at `index`. + * It is needed in a Plugin where the slice itself is const, such as a + * ConstBlockSlice + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(const_slice_type const& slice, + std::size_t index) { + return blaze::column(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Resize the data type to contain at least `size` entries + * + * \param data The data to be resized + * \param new_size New size of data + */ + static inline void resize(data_type& data, std::size_t new_size) { + ELASTICA_ASSERT(data.columns() <= new_size, + "Contract violation, block shrinks"); + return data.resize(Frames::Dimension, new_size, true); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline decltype(auto) slice(data_type& data, + std::size_t start_index, + std::size_t size) { + return blaze::submatrix(data, 0UL, start_index, Frames::Dimension, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from const data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline decltype(auto) slice(data_type const& data, + std::size_t start_index, + std::size_t size) { + return blaze::submatrix(data, 0UL, start_index, Frames::Dimension, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //@} + //************************************************************************ + }; + //************************************************************************** + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/TensorTag.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/TensorTag.hpp new file mode 100644 index 000000000..0dd00491d --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/TensorTag.hpp @@ -0,0 +1,344 @@ +//============================================================================== +/*! +// \file +// \brief Tag indicating a Rank 3 tensor for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include "Configuration/Systems.hpp" // For index checking + +// assert +#include "ErrorHandling/Assert.hpp" + +// Rank +#include "Systems/CosseratRods/Traits/DataType/Protocols.hpp" +#include "Systems/CosseratRods/Traits/DataType/Rank.hpp" + +// Tensor types +#include +#include +#include +#include +#include +#include + +// this comes from the simulator module, which seems like an anti-pattern +// should be in Systems instead +#include "Simulator/Frames.hpp" +#include "Utilities/Math/Vec3.hpp" +#include "Utilities/ProtocolHelpers.hpp" + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Tag for a Rank 3 tensor within Cosserat rods in \elastica + * \ingroup cosserat_rod_traits + * + * \details + * TensorTag is used within the Cosserat rod components of \elastica to + * indicate a tensor. + * + * \tparam Real a floating point type + */ + template + struct TensorTag : public ::tt::ConformsTo { + //**Type definitions****************************************************** + //! The main type for a TaggedTuple + using type = blaze::DynamicTensor; + //! Data type for better readability + using data_type = type; + //! The type of slice for the data type + using slice_type = blaze::Subtensor_; + //! The type of const slice for the data type + using const_slice_type = + const blaze::Subtensor_; + /* + * Developer note: + * Blaze tensor cant initialize with a matrix type + * so the ghosts for now are tensor types which involve + * more work than necessary + * + * using ghost_type = blaze::DynamicMatrix; + * blaze::StaticMatrix; + */ + //! The type of a ghost element for the data type + using ghost_type = + blaze::StaticMatrix; + //! The type of a reference to the underlying data type + using reference_type = blaze::ColumnSlice_; + //! The type of const reference to the underlying data type + using const_reference_type = const blaze::ColumnSlice_; + //! The type of reference to the slice type + using reference_slice_type = + blaze::Submatrix_, blaze::unaligned>; + //! The type of const reference to the slice type + using const_reference_slice_type = + const blaze::Submatrix_, + blaze::unaligned>; + //! The rank of the data type + using rank = Rank<3U>; + //************************************************************************ + + //**Utility functions***************************************************** + /*!\name Utility functions */ + //@{ + + //************************************************************************ + /*!\brief Obtain the ghost value + * + * \details + * Obtains a default value for putting in ghost elements. This function + * can be overriden in case the user wants to. This should not be a + * Real (because implicit conversion fills entire data-structure with + * that particular value) but rather a ghost_type object (like a matrix + * for example) + */ + static inline auto ghost_value() noexcept -> ghost_type { + // cannot be constexpr as it involves allocation :/ + constexpr Real zero(0.0); + constexpr Real identity(1.0); + return ghost_type{{{identity}, {zero}, {zero}}, + {{zero}, {identity}, {zero}}, + {{zero}, {zero}, {identity}}}; + } + //************************************************************************ + + // should really be in the Slice Tag class, but it ends up being + // confusing as what we are really doing is taking a slice of a + // slice but its abstracted away + // template + // static inline constexpr decltype(auto) slice(MT& matrix, + // std::size_t i) + + // TODO decltype expression is same as SliceType + // TODO This is internal function, need not be exposed + // + // + + //************************************************************************ + /*!\brief Obtain slice of the data type at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is typically used in the block to generate BlockSlices or access + * indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(data_type& data, std::size_t index) { + return blaze::columnslice(data, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type at a given index + * + * \details + * Overload to obtain a slice of const data at `index`. + * It is typically used in the block to generate ConstBlockSlices or + * access indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(data_type const& data, + std::size_t index) { + return blaze::columnslice(data, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the slice at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * It is typically used in the block to generate slices that are then + * filled in by ghost values. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type& slice, std::size_t index) { + return blaze::columnslice(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + // return blaze::subtensor(slice, 0UL, 0UL, index, Frames::Dimension, + // Frames::Dimension, 1UL); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-const slice of the const slice at a given index + * + * \details + * Overload to obtain a subslice of const slice at `index`. + * + * It is sometimes needed in a Plugin where the context is const, but the + * data of the underlying structure may well be modified. + * + * \param slice Slice to be further const sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type const& slice, + std::size_t index) { + return blaze::columnslice(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the slice at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * + * It is typically used when lazily generating individual slices from a + * lazily generated slice : for example blockslice.get_position(i), where + * get_position() lazily generates a slice, which then gets further + * sliced. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type&& slice, + std::size_t index) { + return blaze::columnslice(static_cast(slice), index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain sub-slice of the const_slice type at a given index + * + * \details + * Overload to obtain a subslice of const_slice at `index`. + * It is needed in a Plugin where the slice itself is const, such as a + * ConstBlockSlice + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(const_slice_type const& slice, + std::size_t index) { + return blaze::columnslice(slice, index, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Resize the data type to contain at least `size` entries + * + * \param data The data to be resized + * \param new_size New size of data + */ + static inline void resize(data_type& data, std::size_t new_size) { + ELASTICA_ASSERT(data.columns() <= new_size, + "Contract violation, block shrinks"); + // (o, m, n) + return data.resize(Frames::Dimension, Frames::Dimension, new_size, + true); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. With the default size, it is typically used in the + * block to generate slices that are to be filled in by ghost values. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline decltype(auto) slice(data_type& data, + std::size_t start_index, + std::size_t size) { + // (page, row, column), (o, m, n) + return blaze::subtensor(data, 0UL, 0UL, start_index, Frames::Dimension, + Frames::Dimension, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from const data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. With the default size, it is typically used in the + * block to generate slices that are to be filled in by ghost values. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline decltype(auto) slice(data_type const& data, + std::size_t start_index, + std::size_t size) { + // (page, row, column), (o, m, n) + return blaze::subtensor(data, 0UL, 0UL, start_index, Frames::Dimension, + Frames::Dimension, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Assign diagonal elements + * + * @return + */ + //************************************************************************ + static inline void diagonal_assign(slice_type& data, std::size_t index, + Vec3 const& diag) { + auto submatrix = blaze::columnslice( + data, index, blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + submatrix(0UL, 0UL) = static_cast(diag[0UL]); + submatrix(1UL, 1UL) = static_cast(diag[1UL]); + submatrix(2UL, 2UL) = static_cast(diag[2UL]); + } + + static inline void diagonal_assign(data_type& data, std::size_t index, + Vec3 const& diag) { + auto submatrix = blaze::columnslice( + data, index, blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + submatrix(0UL, 0UL) = static_cast(diag[0UL]); + submatrix(1UL, 1UL) = static_cast(diag[1UL]); + submatrix(2UL, 2UL) = static_cast(diag[2UL]); + } + + //@} + //************************************************************************ + }; + //************************************************************************** + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/VectorTag.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/VectorTag.hpp new file mode 100644 index 000000000..38e77d333 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/BlazeBackend/VectorTag.hpp @@ -0,0 +1,274 @@ +//============================================================================== +/*! +// \file +// \brief Tag indicating a vector (Rank 1) for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include "Configuration/Systems.hpp" // For index checking +// assert +#include "ErrorHandling/Assert.hpp" +#include "Systems/CosseratRods/Traits/DataType/Protocols.hpp" +#include "Systems/CosseratRods/Traits/DataType/Rank.hpp" + +// Vector types +#include +#include +#include +#include + +#include "Utilities/ProtocolHelpers.hpp" + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Tag for a Rank 1 tensor within Cosserat rods in \elastica + * \ingroup cosserat_rod_traits + * + * \details + * VectorTag is used within the Cosserat rod components of \elastica to + * indicate a vector. + * + * \tparam Real a floating point type + */ + template + struct VectorTag : public ::tt::ConformsTo { + //**Type definitions****************************************************** + //! The main type for a TaggedTuple + /* + */ + using type = + blaze::DynamicVector /*, TagType*/>; + //! Data type for better readability + using data_type = type; + //! The type of slice for the data type + using slice_type = blaze::Subvector_; + //! The type of a const slice for the data type + using const_slice_type = + const blaze::Subvector_; + //! The type of a ghost element for the data type + using ghost_type = Real; + //! The type of a reference to the underlying data type + using reference_type = Real&; + //! The type of const reference to the underlying data type + using const_reference_type = Real const&; + //! The type of reference to the slice type + using reference_slice_type = Real&; + //! The type of const reference to the slice type + using const_reference_slice_type = Real const&; + //! The rank of the data type + using rank = Rank<1U>; + //************************************************************************ + + //**Utility functions***************************************************** + /*!\name Utility functions */ + //@{ + + //************************************************************************ + /*!\brief Obtain the ghost value + * + * \details + * Obtains a default value for putting in ghost elements. This function + * can be overriden in case the user wants to. This should not be a + * Real (because implicit conversion fills entire data-structure with + * that particular value) but rather a ghost_type object (like a matrix + * for example) + */ + static inline constexpr auto ghost_value() noexcept -> ghost_type { + return ghost_type(0.0); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is typically used in the block to generate BlockSlices or access + * indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(data_type& data, std::size_t index) + -> reference_type { + return data[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type at a given index + * + * \details + * Overload to obtain a slice of const data at `index`. + * It is typically used in the block to generate ConstBlockSlices or + * access indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(data_type const& data, std::size_t index) + -> const_reference_type { + return data[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the slice at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * It is typically used in the block to generate slices that are then + * filled in by ghost values. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(slice_type& slice, std::size_t index) + -> reference_type { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the const slice at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is sometimes needed in a Plugin where the context is const, but the + * data of the underlying structure may well be modified. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(slice_type const& slice, std::size_t index) + -> const_reference_type { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the slice at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * + * It is typically used when lazily generating individual slices from a + * lazily generated slice : for example blockslice.get_position(i), where + * get_position() lazily generates a slice, which then gets further + * sliced. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(slice_type&& slice, + std::size_t index) { + return static_cast(slice)[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain sub-slice of the const_slice type at a given index + * + * \details + * Overload to obtain a subslice of const_slice at `index`. + * It is needed in a Plugin where the slice itself is const, such as a + * ConstBlockSlice + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(const_slice_type const& slice, + std::size_t index) { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Resize the data type to contain at least `size` entries + * + * \param data The data to be resized + * \param new_size New size of data + */ + static inline void resize(data_type& data, std::size_t new_size) { + ELASTICA_ASSERT(data.size() <= new_size, + "Contract violation, block shrinks"); + return data.resize(new_size, true); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline auto slice(data_type& data, std::size_t start_index, + std::size_t size) -> slice_type { + return blaze::subvector(data, start_index, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + static inline auto slice(slice_type& data, std::size_t start_index, + std::size_t size) -> slice_type { + return blaze::subvector(data, start_index, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from const data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline auto slice(data_type const& data, std::size_t start_index, + std::size_t size) -> const_slice_type { + // unaligned always + return blaze::subvector(data, start_index, size, + blaze::ELASTICA_DEFAULT_SYSTEM_INDEX_CHECK); + } + //************************************************************************ + + //@} + //************************************************************************ + }; + //************************************************************************** + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/CMakeLists.txt b/backend/src/Systems/CosseratRods/Traits/DataType/CMakeLists.txt new file mode 100644 index 000000000..837a312a2 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/CMakeLists.txt @@ -0,0 +1,40 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY CosseratRodDataTraits) + +add_elastica_library(${LIBRARY} INTERFACE) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + BlazeBackend/DataTypeTraits.hpp + BlazeBackend/MatrixTag.hpp + BlazeBackend/TensorTag.hpp + BlazeBackend/VectorTag.hpp + Protocols.hpp + Rank.hpp + ScalarTag.hpp +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + BlazeBackend/DataTypeTraits.hpp + BlazeBackend/MatrixTag.hpp + BlazeBackend/TensorTag.hpp + BlazeBackend/VectorTag.hpp + Protocols.hpp + Rank.hpp + ScalarTag.hpp +) + + +target_link_libraries( + ${LIBRARY} + INTERFACE + ModuleSettings + Utilities + Blaze + BlazeTensor +) + diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/Protocols.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/Protocols.hpp new file mode 100644 index 000000000..500c89297 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/Protocols.hpp @@ -0,0 +1,197 @@ +//============================================================================== +/*! +// \file +// \brief Data-trait protocols +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include // for size_t +#include // for declval + +#include "Utilities/TypeTraits.hpp" + +namespace elastica { + + namespace cosserat_rod { + + namespace protocols { + + //======================================================================== + // + // CLASS DEFINITION + // + //======================================================================== + + //************************************************************************ + /*!\brief Data-trait protocol. + * \ingroup cosserat_rod_traits + * + * Class to enforce adherence to a Cosserat Rod Data-trait protocol. Any + * valid data trait class within the \elastica library should (publicly) + * inherit from this class to indicate it qualifies as a data-trait. Only + * in case a class is derived publicly from this base class, the + * tt::conforms_to and tt::assert_conforms_to type traits recognizes the + * class as valid data traits. + * + * Requires that a conforming type `ConformingType` has these nested types + * \snippet this expected_types + * and these static member functions, + * \snippet this expected_static_functions + * + * \example + * The following shows an example of minimal conformance to this protocol. + * With the setup shown here + * \snippet Mocks/CosseratRodTraits.hpp vector_datatrait + * we ensure conformance as + * \snippet DataType/Test_Protocols.cpp datatrait_protocol_eg + */ + struct DataTrait { + //********************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Auxiliary helper struct for enforcing protocols. + // \ingroup protocols + */ + template + struct test { + public: + /// [expected_types] + //**Type definitions************************************************** + //! The main type for a TaggedTuple + using data_type = typename ConformingType::type; + //! The type of slice for the data type + using slice_type = typename ConformingType::slice_type; + //! The type of const slice for the data type + using const_slice_type = typename ConformingType::const_slice_type; + //! The type of a reference to the underlying data type + using reference_type = typename ConformingType::reference_type; + //! The type of const reference to the underlying data type + using const_reference_type = + typename ConformingType::const_reference_type; + //! The type of reference to the slice type + using reference_slice_type = + typename ConformingType::reference_slice_type; + //! The type of const reference to the slice type + using const_reference_slice_type = + typename ConformingType::const_reference_slice_type; + //! The type of a ghost element for the data type + using ghost_type = typename ConformingType::ghost_type; + // //! The rank of the the data type + // using rank = typename ConformingType::rank; + //******************************************************************** + /// [expected_types] + + /// [expected_static_functions] + // Function prototypes needing no arguments + using ghost_value_return_type = + decltype(ConformingType::ghost_value()); + static_assert(cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`ghost_type ghost_value()`)error"); + + // Function prototypes needing two arguments + using index = std::size_t; + using slice_at_index_return_type = decltype(ConformingType::slice( + std::declval(), std::declval())); + static_assert( + cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`auto slice(data_type& data, std::size_t index)` +)error"); + + using const_slice_at_index_return_type = + decltype(ConformingType::slice(std::declval(), + std::declval())); + static_assert(cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`auto slice(data_type const& data, std::size_t index)` +)error"); + + using slice_of_slice_at_index_return_type = + decltype(ConformingType::slice(std::declval(), + std::declval())); + static_assert(cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`auto slice(slice_type& slice, std::size_t index)` +)error"); + + using slice_of_const_slice_at_index_return_type = + decltype(ConformingType::slice(std::declval(), + std::declval())); + static_assert( + cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`auto slice(slice_type const& slice, std::size_t index)` +)error"); + + using const_slice_of_const_slice_at_index_return_type = + decltype(ConformingType::slice( + std::declval(), + std::declval())); + static_assert( + cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`auto slice(const_slice_type const& slice, std::size_t index)` +)error"); + + using new_dofs = std::size_t; + using resize_return_type = decltype(ConformingType::resize( + std::declval(), std::declval())); + static_assert(cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`void resize(data_type& data, std::size_t new_dofs)`)error"); + + // Function prototypes needing three arguments + using start = std::size_t; + using size = std::size_t; + using slice_return_type = decltype(ConformingType::slice( + std::declval(), std::declval(), + std::declval())); + static_assert(cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`slice_type slice(data_type& data, std::size_t start, std::size_t size)`)error"); + + using const_slice_return_type = decltype(ConformingType::slice( + std::declval(), std::declval(), + std::declval())); + static_assert( + cpp17::is_same_v, + R"error( +Not a conforming data trait, doesn't properly implement static function +`const_slice_type slice(data_type const& data, std::size_t start, std::size_t size)`)error"); + + /// [expected_static_functions] + }; + /*! \endcond */ + //********************************************************************** + }; + //************************************************************************ + + } // namespace protocols + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/Rank.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/Rank.hpp new file mode 100644 index 000000000..85ecc7eda --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/Rank.hpp @@ -0,0 +1,40 @@ +//============================================================================== +/*! +// \file +// \brief Marks a rank of a tensor +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +#include + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Marks a rank of a tensor + // \ingroup cosserat_rod + */ + template + struct Rank : std::integral_constant {}; + /*! \endcond */ + //************************************************************************** + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/DataType/ScalarTag.hpp b/backend/src/Systems/CosseratRods/Traits/DataType/ScalarTag.hpp new file mode 100644 index 000000000..6e22fa8ff --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/DataType/ScalarTag.hpp @@ -0,0 +1,235 @@ +//============================================================================== +/*! +// \file +// \brief Tag indicating a vector (Rank 0) for Cosserat Rods within \elastica +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +// Vector types +#include // ref wrapper +#include + +// assert +#include "ErrorHandling/Assert.hpp" +#include "Systems/CosseratRods/Traits/DataTraits/Protocols.hpp" +#include "Systems/CosseratRods/Traits/DataTraits/Rank.hpp" +#include "Utilities/NoSuchType.hpp" + +namespace elastica { + + namespace cosserat_rod { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Tag for a Rank 1 tensor within Cosserat rods in \elastica + * \ingroup cosserat_rod_traits + * + * \details + * ScalarTag is used within the Cosserat rod components of \elastica to + * indicate a vector. + * + * \tparam RealOrIndex a floating point (or) index type + */ + template + struct ScalarTag : public ::tt::ConformsTo { + //**Type definitions****************************************************** + //! The main type for a TaggedTuple + using type = std::vector; + //! Data type for better readability + using data_type = type; + //! The type of slice for the data type + using slice_type = std::vector>; + //! The type of a const slice for the data type + using const_slice_type = + const std::vector>; + //! The type of a ghost element for the data type + using ghost_type = NoSuchType; + //! The type of a reference to the underlying data type + using reference_type = RealOrIndex&; + //! The type of const reference to the underlying data type + using const_reference_type = RealOrIndex const&; + //! The type of reference to the slice type + using reference_slice_type = RealOrIndex&; + //! The type of const reference to the slice type + using const_reference_slice_type = RealOrIndex const&; + //! The rank of the data type + using rank = Rank<0U>; + //************************************************************************ + + //**Utility functions***************************************************** + /*!\name Utility functions */ + //@{ + + //************************************************************************ + /*!\brief Obtain the ghost value + * + * \details + * Obtains a default value for putting in ghost elements. This function + * can be overriden in case the user wants to. This should not be a + * RealOrIndex (because implicit conversion fills entire data-structure + * with that particular value) but rather a ghost_type object (like a + * matrix for example) + */ + static inline constexpr auto ghost_value() noexcept -> ghost_type { + return {}; + } + //********************************************************************** + + //************************************************************************ + /*!\brief Obtain slice of the data type at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is typically used in the block to generate BlockSlices or access + * indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(data_type& data, std::size_t index) + -> reference_type { + return data[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type at a given index + * + * \details + * Overload to obtain a slice of const data at `index`. + * It is typically used in the block to generate ConstBlockSlices or + * access indexed data + * + * \param data Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(data_type const& data, std::size_t index) + -> const_reference_type { + return data[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the slice at a given index + * + * \details + * Overload to obtain a subslice of slice at `index`. + * It is typically used in the block to generate slices that are then + * filled in by ghost values. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(slice_type& slice, std::size_t index) + -> reference_type { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain a sub-slice of the const slice at a given index + * + * \details + * Overload to obtain a slice of data at `index`. + * It is sometimes needed in a Plugin where the context is const, but the + * data of the underlying structure may well be modified. + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline auto slice(slice_type const& slice, std::size_t index) + -> const_reference_type { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain sub-slice of the const_slice type at a given index + * + * \details + * Overload to obtain a subslice of const_slice at `index`. + * It is needed in a Plugin where the slice itself is const, such as a + * ConstBlockSlice + * + * \param slice Slice to be further sliced + * \param index Index of slicing + */ + static inline decltype(auto) slice(const_slice_type const& slice, + std::size_t index) { + return slice[index]; + } + //************************************************************************ + + //************************************************************************ + /*!\brief Resize the data type to contain at least `size` entries + * + * \param data The data to be resized + * \param new_size New size of data + */ + static inline void resize(data_type& data, std::size_t new_size) { + ELASTICA_ASSERT(data.size() <= new_size, + "Contract violation, block shrinks"); + return data.resize(new_size); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline auto slice(data_type& data, std::size_t start_index, + std::size_t size) -> slice_type { + return slice_type(&data[start_index], &data[start_index + size]); + } + //************************************************************************ + + //************************************************************************ + /*!\brief Obtain slice of the const data type through a set of indices + * + * \details + * Overload to obtain a slice of `size` amount from const data at `start`. + * This is typically used in the block to generate slices that are to be + * filled in by the class hierarchies, when filling in a block + * incrementally. + * + * \param data Data to be sliced + * \param start_index Start index + * \param size Size of the data + */ + static inline auto slice(data_type const& data, std::size_t start_index, + std::size_t size) -> const_slice_type { + // unaligned always + return const_slice_type(&data[start_index], &data[start_index + size]); + } + //************************************************************************ + + //@} + //************************************************************************ + }; + //************************************************************************** + + } // namespace cosserat_rod + +} // namespace elastica diff --git a/backend/src/Systems/CosseratRods/Traits/Types.hpp b/backend/src/Systems/CosseratRods/Traits/Types.hpp new file mode 100644 index 000000000..f34ca9c75 --- /dev/null +++ b/backend/src/Systems/CosseratRods/Traits/Types.hpp @@ -0,0 +1,35 @@ +//============================================================================== +/*! +// \file +// \brief All trait types belonging to a Cosserat rod +// +// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved +// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved +// +// Distributed under the MIT License. +// See LICENSE.txt for details. +*/ +//============================================================================== + +#pragma once + +namespace elastica { + + namespace cosserat_rod { + + //////////////////////////////////////////////////////////////////////////// + // + // Forward declarations of cosserat rod trait types + // + //////////////////////////////////////////////////////////////////////////// + + //************************************************************************ + /*! \cond ELASTICA_INTERNAL */ + struct DataOpsTraits; + // struct PlacementTrait; + /*! \endcond */ + //************************************************************************ + + } // namespace cosserat_rod + +} // namespace elastica From ecedea500f06e32de6dc93aa7f11435a41a39d96 Mon Sep 17 00:00:00 2001 From: Ankith Date: Sun, 23 Jun 2024 17:37:13 +0530 Subject: [PATCH 07/62] Wrap and test everything in BlazeLinearAlgebra.hpp --- backend/benchmarking/matmul.py | 6 +- backend/src/_linalg.cpp | 271 ++++++++++++++++++++++++++++++++- backend/src/meson.build | 2 +- backend/tests/test_linalg.py | 113 ++++++++++++++ 4 files changed, 388 insertions(+), 4 deletions(-) create mode 100644 backend/tests/test_linalg.py diff --git a/backend/benchmarking/matmul.py b/backend/benchmarking/matmul.py index 3235e0ef6..04e4a2b18 100644 --- a/backend/benchmarking/matmul.py +++ b/backend/benchmarking/matmul.py @@ -1,4 +1,4 @@ -from elasticapp._linalg import batch_matmul_naive, batch_matmul_blaze +from elasticapp._linalg import batch_matmul_naive, batch_matmul_blaze, batch_matmul from elastica._linalg import _batch_matmul import numpy import time @@ -29,14 +29,16 @@ def benchmark_batchsize(funcs: list, batches: list[int], num_iterations: int = 1 results = benchmark_batchsize( - [batch_matmul_naive, batch_matmul_blaze, _batch_matmul], [2**i for i in range(14)] + [batch_matmul_naive, batch_matmul_blaze, batch_matmul, _batch_matmul], [2**i for i in range(14)] ) for size, data in results.items(): pyelastica = data["_batch_matmul"] elasticapp = data["batch_matmul_naive"] elasticapp_blaze = data["batch_matmul_blaze"] + elasticapp_blaze_v2 = data["batch_matmul"] print(f"{size = }") print(f"{pyelastica = }") print(f"{elasticapp = }, ratio: {elasticapp / pyelastica}") print(f"{elasticapp_blaze = }, ratio: {elasticapp_blaze / pyelastica}") + print(f"{elasticapp_blaze_v2 = }, ratio: {elasticapp_blaze_v2 / pyelastica}") print() diff --git a/backend/src/_linalg.cpp b/backend/src/_linalg.cpp index b628ebd2a..651e29177 100644 --- a/backend/src/_linalg.cpp +++ b/backend/src/_linalg.cpp @@ -2,8 +2,12 @@ #include #include +#include "Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp" + namespace py = pybind11; +using blaze::DynamicMatrix; +using blaze::DynamicTensor; using blaze::StaticMatrix; /* Simple rewrite of elastica._linalg._batch_matmul */ @@ -75,6 +79,234 @@ py::array_t batch_matmul_blaze( return output_matrix; } +py::array_t difference_kernel(py::array_t vector_batch) +{ + const auto v0 = vector_batch.shape(0); + const auto num_elems = vector_batch.shape(1); + + assert(v0 == 3UL); + auto output_arr = py::array_t{{v0, num_elems - 1}}; + + auto vector_batch_unchecked = vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicMatrix blaze_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); + } + } + auto blaze_output = elastica::difference_kernel(blaze_vector_batch); + for (py::ssize_t j = 0; j < num_elems - 1; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_matvec( + py::array_t matrix_collection, py::array_t vector_collection) +{ + const auto v0 = vector_collection.shape(0); + const auto num_elems = vector_collection.shape(1); + assert(v0 == 3UL); + + assert(matrix_collection.shape(0) == 3UL); + assert(matrix_collection.shape(1) == 3UL); + assert(matrix_collection.shape(2) == num_elems); + + auto output_arr = py::array_t{{v0, num_elems}}; + + auto matrix_collection_unchecked = matrix_collection.unchecked<3>(); + auto vector_collection_unchecked = vector_collection.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicTensor blaze_matrix_collection(v0, v0, num_elems); + DynamicMatrix blaze_vector_collection(v0, num_elems); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_matrix_collection(i, j, k) = matrix_collection_unchecked(i, j, k); + } + blaze_vector_collection(j, k) = vector_collection_unchecked(j, k); + } + } + + DynamicMatrix blaze_output(v0, num_elems); + elastica::batch_matvec(blaze_output, blaze_matrix_collection, blaze_vector_collection); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_matmul( + py::array_t first_matrix_batch, py::array_t second_matrix_batch) +{ + const auto m0 = first_matrix_batch.shape(0); + const auto m1 = first_matrix_batch.shape(1); + const auto num_elems = first_matrix_batch.shape(2); + assert(m0 == 3UL); + assert(m1 == 3UL); + + assert(second_matrix_batch.shape(0) == 3UL); + assert(second_matrix_batch.shape(1) == 3UL); + assert(second_matrix_batch.shape(2) == num_elems); + + auto output_arr = py::array_t{{m0, m1, num_elems}}; + + auto first_matrix_batch_unchecked = first_matrix_batch.unchecked<3>(); + auto second_matrix_batch_unchecked = second_matrix_batch.unchecked<3>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<3>(); + + DynamicTensor blaze_first_matrix_batch(m0, m1, num_elems); + DynamicTensor blaze_second_matrix_batch(m0, m1, num_elems); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_matrix_batch(i, j, k) = first_matrix_batch_unchecked(i, j, k); + blaze_second_matrix_batch(i, j, k) = second_matrix_batch_unchecked(i, j, k); + } + } + } + DynamicTensor blaze_output(m0, m1, num_elems); + elastica::batch_matmul(blaze_output, blaze_first_matrix_batch, blaze_second_matrix_batch); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j, k) = blaze_output(i, j, k); + } + } + } + return output_arr; +} + +py::array_t batch_cross( + py::array_t first_vector_batch, py::array_t second_vector_batch) +{ + const auto v0 = first_vector_batch.shape(0); + const auto num_elems = first_vector_batch.shape(1); + assert(v0 == 3UL); + + assert(second_vector_batch.shape(0) == 3UL); + assert(second_vector_batch.shape(1) == num_elems); + + auto output_arr = py::array_t{{v0, num_elems}}; + + auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); + auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicMatrix blaze_first_vector_batch(v0, num_elems); + DynamicMatrix blaze_second_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); + blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); + } + } + + DynamicMatrix blaze_output(v0, num_elems); + elastica::batch_cross(blaze_output, blaze_first_vector_batch, blaze_second_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_dot( + py::array_t first_vector_batch, py::array_t second_vector_batch) +{ + const auto v0 = first_vector_batch.shape(0); + const auto num_elems = first_vector_batch.shape(1); + assert(v0 == 3UL); + + assert(second_vector_batch.shape(0) == 3UL); + assert(second_vector_batch.shape(1) == num_elems); + + auto output_arr = py::array_t{num_elems}; + + auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); + auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); + + DynamicMatrix blaze_first_vector_batch(v0, num_elems); + DynamicMatrix blaze_second_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); + blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); + } + } + + auto blaze_output = elastica::batch_dot(blaze_first_vector_batch, blaze_second_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + output_arr_unchecked(j) = blaze_output[j]; + } + + return output_arr; +} + +py::array_t batch_norm(py::array_t vector_batch) +{ + const auto v0 = vector_batch.shape(0); + const auto num_elems = vector_batch.shape(1); + + assert(v0 == 3UL); + + auto output_arr = py::array_t{num_elems}; + + auto vector_batch_unchecked = vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); + + DynamicMatrix blaze_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); + } + } + + auto blaze_output = elastica::batch_norm(blaze_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + output_arr_unchecked(j) = blaze_output[j]; + } + + return output_arr; +} + PYBIND11_MODULE(_linalg, m) { m.doc() = R"pbdoc( @@ -87,6 +319,13 @@ PYBIND11_MODULE(_linalg, m) :toctree: _generate batch_matmul_naive + batch_matmul_blaze + difference_kernel + batch_matvec + batch_matmul + batch_cross + batch_dot + batch_norm )pbdoc"; m.def("batch_matmul_naive", &batch_matmul_naive, R"pbdoc( @@ -99,9 +338,39 @@ PYBIND11_MODULE(_linalg, m) of 3x3 matrices can be multiplied. )pbdoc"); + m.def("difference_kernel", &difference_kernel, R"pbdoc( + Vector Difference + )pbdoc"); + + m.def("batch_matvec", &batch_matvec, R"pbdoc( + Computes a batchwise matrix-vector product given in indical notation: + matvec_batch{ik} = matrix_batch{ijk} * vector_batch{jk} + )pbdoc"); + + m.def("batch_matmul", &batch_matmul, R"pbdoc( + Computes a batchwise matrix-matrix product given in indical notation: + matmul_batch{ilk} = first_matrix_batch{ijk} * second_matrix_batch{jlk} + )pbdoc"); + + m.def("batch_cross", &batch_cross, R"pbdoc( + Computes a batchwise vector-vector cross product given in indical notation: + cross_batch{il} = LCT{ijk} * first_vector_batch{jl} * second_vector_batch{kl} + where LCT is the Levi-Civita Tensor + )pbdoc"); + + m.def("batch_dot", &batch_dot, R"pbdoc( + Computes a batchwise vector-vector dot product given in indical notation: + dot_batch{j} = first_vector_batch{ij} * second_vector_batch{ij} + )pbdoc"); + + m.def("batch_norm", &batch_norm, R"pbdoc( + Computes a batchwise vector L2 norm given in indical notation: + norm_batch{j} = (vector_batch{ij} * vector_batch{ij})^0.5 + )pbdoc"); + #ifdef VERSION_INFO m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); #else m.attr("__version__") = "dev"; #endif -} \ No newline at end of file +} diff --git a/backend/src/meson.build b/backend/src/meson.build index 9fb245448..690dd2722 100644 --- a/backend/src/meson.build +++ b/backend/src/meson.build @@ -10,7 +10,7 @@ example_1 = py.extension_module( _linalg = py.extension_module( '_linalg', sources: ['_linalg.cpp'], - dependencies: [py_dep, pybind11_dep, blaze_dep], + dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], link_language: 'cpp', install: true, subdir: package, diff --git a/backend/tests/test_linalg.py b/backend/tests/test_linalg.py new file mode 100644 index 000000000..1fba0332b --- /dev/null +++ b/backend/tests/test_linalg.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 + +# This file is based on pyelastica tests/test_math/test_linalg.py + +# System imports +import numpy as np +import pytest +from numpy.testing import assert_allclose +from elasticapp._linalg import ( + difference_kernel, + batch_matvec, + batch_matmul, + batch_cross, + batch_dot, + batch_norm, +) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_difference_kernel(blocksize: int): + input_vector_collection = np.random.randn(3, blocksize) + output_vector_collection = difference_kernel(input_vector_collection) + + correct_vector_collection = ( + input_vector_collection[:, 1:] - input_vector_collection[:, :-1] + ) + + assert_allclose(output_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_matvec(blocksize: int): + input_matrix_collection = np.random.randn(3, 3, blocksize) + input_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_matvec( + input_matrix_collection, input_vector_collection + ) + + correct_vector_collection = [ + np.dot(input_matrix_collection[..., i], input_vector_collection[..., i]) + for i in range(blocksize) + ] + correct_vector_collection = np.array(correct_vector_collection).T + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_matmul(blocksize: int): + input_first_matrix_collection = np.random.randn(3, 3, blocksize) + input_second_matrix_collection = np.random.randn(3, 3, blocksize) + + test_matrix_collection = batch_matmul( + input_first_matrix_collection, input_second_matrix_collection + ) + + correct_matrix_collection = np.empty((3, 3, blocksize)) + for i in range(blocksize): + correct_matrix_collection[..., i] = np.dot( + input_first_matrix_collection[..., i], + input_second_matrix_collection[..., i], + ) + + assert_allclose(test_matrix_collection, correct_matrix_collection) + + +# TODO : Generalize to two dimensions +@pytest.mark.parametrize("dim", [3]) +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_cross(dim, blocksize: int): + input_first_vector_collection = np.random.randn(dim, blocksize) + input_second_vector_collection = np.random.randn(dim, blocksize) + + test_vector_collection = batch_cross( + input_first_vector_collection, input_second_vector_collection + ) + correct_vector_collection = np.cross( + input_first_vector_collection, input_second_vector_collection, axisa=0, axisb=0 + ).T + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_dot(blocksize: int): + input_first_vector_collection = np.random.randn(3, blocksize) + input_second_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_dot( + input_first_vector_collection, input_second_vector_collection + ) + + correct_vector_collection = np.einsum( + "ij,ij->j", input_first_vector_collection, input_second_vector_collection + ) + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_norm(blocksize: int): + input_first_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_norm(input_first_vector_collection) + + correct_vector_collection = np.sqrt( + np.einsum( + "ij,ij->j", input_first_vector_collection, input_first_vector_collection + ) + ) + + assert_allclose(test_vector_collection, correct_vector_collection) From d6128171538356335f15e6ec269ac4f24437deff Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 24 Jun 2024 14:03:41 -0500 Subject: [PATCH 08/62] python-binding template for vector/matrix/tensor --- .../src/Utilities/Math/Python/Bindings.cpp | 31 +++ .../src/Utilities/Math/Python/BlazeMatrix.cpp | 106 ++++++++ .../src/Utilities/Math/Python/BlazeTensor.cpp | 111 +++++++++ .../src/Utilities/Math/Python/BlazeVector.cpp | 230 ++++++++++++++++++ .../src/Utilities/Math/Python/CMakeLists.txt | 28 +++ .../Utilities/Math/Python/SliceHelpers.hpp | 108 ++++++++ 6 files changed, 614 insertions(+) create mode 100644 backend/src/Utilities/Math/Python/Bindings.cpp create mode 100644 backend/src/Utilities/Math/Python/BlazeMatrix.cpp create mode 100644 backend/src/Utilities/Math/Python/BlazeTensor.cpp create mode 100644 backend/src/Utilities/Math/Python/BlazeVector.cpp create mode 100644 backend/src/Utilities/Math/Python/CMakeLists.txt create mode 100644 backend/src/Utilities/Math/Python/SliceHelpers.hpp diff --git a/backend/src/Utilities/Math/Python/Bindings.cpp b/backend/src/Utilities/Math/Python/Bindings.cpp new file mode 100644 index 000000000..c2ff0a134 --- /dev/null +++ b/backend/src/Utilities/Math/Python/Bindings.cpp @@ -0,0 +1,31 @@ +//****************************************************************************** +// Includes +//****************************************************************************** +#include + +namespace py = pybind11; + +namespace py_bindings { + void bind_blaze_vector(py::module& m); // NOLINT + // void bind_blaze_index_vector(py::module& m); // NOLINT + // void bind_blaze_subvector(py::module& m); // NOLINT + // void bind_blaze_index_subvector(py::module& m); // NOLINT + void bind_blaze_matrix(py::module& m); // NOLINT + // void bind_blaze_submatrix(py::module& m); // NOLINT + void bind_blaze_tensor(py::module& m); // NOLINT + // void bind_blaze_subtensor(py::module& m); // NOLINT +} // namespace py_bindings + +PYBIND11_MODULE(_PyArrays, m) { // NOLINT + m.doc() = R"pbdoc( + Bindings for Elastica++ array types + )pbdoc"; + py_bindings::bind_blaze_vector(m); + // py_bindings::bind_blaze_index_vector(m); + // py_bindings::bind_blaze_subvector(m); + // py_bindings::bind_blaze_index_subvector(m); + py_bindings::bind_blaze_matrix(m); + // py_bindings::bind_blaze_submatrix(m); + py_bindings::bind_blaze_tensor(m); + // py_bindings::bind_blaze_subtensor(m); +} diff --git a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp new file mode 100644 index 000000000..9c3fb654e --- /dev/null +++ b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp @@ -0,0 +1,106 @@ +//****************************************************************************** +// Includes +//****************************************************************************** + +// +// #include "PythonBindings/BoundChecks.hpp" +// +#include "Utilities/DefineTypes.h" +// #include "Utilities/MakeString.hpp" +// +#include "Utilities/Math/Python/SliceHelpers.hpp" +// +#include +#include +#include +#include +#include +#include +// +#include +#include +// +#include + +namespace py = pybind11; + +namespace py_bindings { + + //**************************************************************************** + /*!\brief Helps bind a matrix type in \elastica + * \ingroup python_bindings + */ + void bind_blaze_matrix(py::module& m) { // NOLINT + using Real = ::elastica::real_t; + using type = ::blaze::DynamicMatrix>; + + // Wrapper for basic type operations + py::class_(m, "Matrix", py::buffer_protocol()) + .def(py::init(), py::arg("rows"), + py::arg("columns")) + .def(py::init([](py::buffer buffer) { + py::buffer_info info = buffer.request(); + // Sanity-check the buffer + if (info.format != py::format_descriptor::format()) { + throw std::runtime_error( + "Incompatible format: expected a Real array."); + } + if (info.ndim != 2) { + throw std::runtime_error("Incompatible dimension."); + } + const auto rows = static_cast(info.shape[0]); + const auto columns = static_cast(info.shape[1]); + auto data = static_cast(info.ptr); + return type(rows, columns, data); + }), + py::arg("buffer")) + // Expose the data as a Python buffer so it can be cast into Numpy + // arrays + .def_buffer([](type& matrix) { + return py::buffer_info( + matrix.data(), + // Size of one scalar + sizeof(Real), py::format_descriptor::format(), + // Number of dimensions + 2, + // Size of the buffer + {matrix.rows(), matrix.columns()}, + // Stride for each index (in bytes). Data is stored + // in row-major layout (see `type.hpp`). + {sizeof(Real) * matrix.spacing(), sizeof(Real)}); + }) + .def_property_readonly( + "shape", + +[](const type& self) { + return std::tuple(self.rows(), + self.columns()); + }) + // __getitem__ and __setitem__ are the subscript operators (M[*,*]). + .def( + "__getitem__", + +[](const type& self, + const std::tuple& x) { + // matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); + return self(std::get<0>(x), std::get<1>(x)); + }) + .def( + "__getitem__", + +[](type& t, std::tuple const slice) { + return array_slice(t, std::move(slice)); + }) + // Need __str__ for converting to string/printing + .def( + "__str__", + +[](const type& self) { return std::string(MakeString{} << self); }) + .def( + "__setitem__", + +[](type& self, const std::tuple& x, + const Real val) { + // matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); + self(std::get<0>(x), std::get<1>(x)) = val; + }); + } + //**************************************************************************** + +} // namespace py_bindings diff --git a/backend/src/Utilities/Math/Python/BlazeTensor.cpp b/backend/src/Utilities/Math/Python/BlazeTensor.cpp new file mode 100644 index 000000000..57dc73a4d --- /dev/null +++ b/backend/src/Utilities/Math/Python/BlazeTensor.cpp @@ -0,0 +1,111 @@ +//****************************************************************************** +// Includes +//****************************************************************************** + + +// #include "PythonBindings/BoundChecks.hpp" +// +#include "Utilities/DefineTypes.h" +// #include "Utilities/MakeString.hpp" +// +#include "Utilities/Math/Python/SliceHelpers.hpp" +// +#include +#include +#include +#include +#include +#include +// +#include +#include +// +#include + +namespace py = pybind11; + +namespace py_bindings { + + //**************************************************************************** + /*!\brief Helps bind a tensor type in \elastica + * \ingroup python_bindings + */ + void bind_blaze_tensor(py::module& m) { // NOLINT + using Real = ::elastica::real_t; + using type = ::blaze::DynamicTensor; + + // Wrapper for basic type operations + py::class_(m, "Tensor", py::buffer_protocol()) + .def(py::init(), + py::arg("pages"), py::arg("rows"), py::arg("columns")) + .def(py::init([](py::buffer buffer) { + py::buffer_info info = buffer.request(); + // Sanity-check the buffer + if (info.format != py::format_descriptor::format()) { + throw std::runtime_error( + "Incompatible format: expected a Real array."); + } + if (info.ndim != 3) { + throw std::runtime_error("Incompatible dimension."); + } + const auto pages = static_cast(info.shape[0]); + const auto rows = static_cast(info.shape[1]); + const auto columns = static_cast(info.shape[2]); + auto data = static_cast(info.ptr); + return type(pages, rows, columns, data); + }), + py::arg("buffer")) + // Expose the data as a Python buffer so it can be cast into Numpy + // arrays + .def_buffer([](type& tensor) { + return py::buffer_info( + tensor.data(), + // Size of one scalar + sizeof(Real), py::format_descriptor::format(), + // Number of dimensions + 3, + // Size of the buffer + {tensor.pages(), tensor.rows(), tensor.columns()}, + // Stride for each index (in bytes). Data is stored + // in column-major layout (see `type.hpp`). + {sizeof(Real) * tensor.rows() * tensor.spacing(), + sizeof(Real) * tensor.spacing(), sizeof(Real)}); + }) + .def_property_readonly( + "shape", + +[](const type& self) { + return std::tuple( + self.pages(), self.rows(), self.columns()); + }) + // __getitem__ and __setitem__ are the subscript operators (M[*,*]). + .def( + "__getitem__", + +[](const type& self, + const std::tuple& x) { + // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), + std::get<2>(x)); + return self(std::get<0>(x), std::get<1>(x), std::get<2>(x)); + }) + .def( + "__getitem__", + +[](type& self, std::tuple slice) { + return array_slice(self, std::move(slice)); + }) + .def( + "__setitem__", + +[](type& self, + const std::tuple& x, + const Real val) { + // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), + std::get<2>(x)); + self(std::get<0>(x), std::get<1>(x), std::get<2>(x)) = val; + }) + // Need __str__ for converting to string/printing + .def( + "__str__", +[](const type& self) { + return std::string(MakeString{} << self); + }); + } + //**************************************************************************** + +} // namespace py_bindings diff --git a/backend/src/Utilities/Math/Python/BlazeVector.cpp b/backend/src/Utilities/Math/Python/BlazeVector.cpp new file mode 100644 index 000000000..adedad39f --- /dev/null +++ b/backend/src/Utilities/Math/Python/BlazeVector.cpp @@ -0,0 +1,230 @@ +//****************************************************************************** +// Includes +//****************************************************************************** + +// #include "PythonBindings/BoundChecks.hpp" +// +#include "Utilities/DefineTypes.h" +// +#include "Utilities/Math/Python/SliceHelpers.hpp" +// +#include +#include +#include +#include +#include +// +#include +#include +#include +// +#include +// + +namespace py = pybind11; + +namespace py_bindings { + + //**************************************************************************** + /*!\brief Helps bind a vector type in \elastica + * \ingroup python_bindings + */ + void bind_blaze_vector(py::module& m) { // NOLINT + + using Real = ::elastica::real_t; + using type = ::blaze::DynamicVector>; + // Wrapper for basic DataVector operations + auto py_vector = + py::class_(m, "Vector", py::buffer_protocol()) + .def(py::init(), py::arg("size")) + .def(py::init(), py::arg("size"), + py::arg("fill")) + .def(py::init([](std::vector const& values) { + type result(values.size()); + std::copy(values.begin(), values.end(), result.begin()); + return result; + }), + py::arg("values")) + .def(py::init([](py::buffer buffer) { + py::buffer_info info = buffer.request(); + // Sanity-check the buffer + if (info.format != py::format_descriptor::format()) { + throw std::runtime_error( + "Incompatible format: expected a Real array"); + } + if (info.ndim != 1) { + throw std::runtime_error("Incompatible dimension."); + } + const auto size = static_cast(info.shape[0]); + auto data = static_cast(info.ptr); + type result(size); + std::copy_n(data, result.size(), result.begin()); + return result; + }), + py::arg("buffer")) + // Expose the data as a Python buffer so it can be cast into Numpy + // arrays + .def_buffer([](type& t) { + return py::buffer_info(t.data(), + // Size of one scalar + sizeof(Real), + py::format_descriptor::format(), + // Number of dimensions + 1, + // Size of the buffer + {t.size()}, + // Stride for each index (in bytes) + {sizeof(Real)}); + }) + .def( + "__iter__", + [](const type& t) { + return py::make_iterator(t.begin(), t.end()); + }, + // Keep object alive while iterator exists + py::keep_alive<0, 1>()) + // __len__ is for being able to write len(my_data_vector) in python + .def("__len__", [](const type& t) { return t.size(); }) + .def_property_readonly( + "shape", + +[](const type& t) { + return std::tuple(t.size()); + }) + // __getitem__ and __setitem__ are the subscript operators + // (operator[]). To define (and overload) operator() use __call__ + .def( + "__getitem__", + +[](const type& t, const size_t i) { + // bounds_check(t, i); + return t[i]; + }) + .def( + "__getitem__", + +[](type& t, const py::slice slice) { + return array_slice(t, std::move(slice)); + }) + .def( + "__setitem__", +[](type& t, const size_t i, const Real v) { + // bounds_check(t, i); + t[i] = v; + }); + + // Need __str__ for converting to string + py_vector + .def("__str__") + // repr allows you to output the object in an interactive python + // terminal using obj to get the "string REPResenting the object". + .def("__repr__") + .def(py::self += py::self) + // Need to do math explicitly converting to DataVector because we don't + // want to represent all the possible expression template types + .def( + "abs", +[](const type& t) { return type{abs(t)}; }) + .def( + "acos", +[](const type& t) { return type{acos(t)}; }) + .def( + "acosh", +[](const type& t) { return type{acosh(t)}; }) + .def( + "asin", +[](const type& t) { return type{asin(t)}; }) + .def( + "asinh", +[](const type& t) { return type{asinh(t)}; }) + .def( + "atan", +[](const type& t) { return type{atan(t)}; }) + .def( + "atan2", + +[](const type& y, const type& x) { return type{atan2(y, x)}; }) + .def( + "atanh", +[](const type& t) { return type{atanh(t)}; }) + .def( + "cbrt", +[](const type& t) { return type{cbrt(t)}; }) + .def( + "cos", +[](const type& t) { return type{cos(t)}; }) + .def( + "cosh", +[](const type& t) { return type{cosh(t)}; }) + .def( + "erf", +[](const type& t) { return type{erf(t)}; }) + .def( + "erfc", +[](const type& t) { return type{erfc(t)}; }) + .def( + "exp", +[](const type& t) { return type{exp(t)}; }) + .def( + "exp2", +[](const type& t) { return type{exp2(t)}; }) + .def( + "exp10", +[](const type& t) { return type{exp10(t)}; }) + // .def( + // "fabs", +[](const type& t) { return type{fabs(t)}; }) + .def( + "hypot", + +[](const type& x, const type& y) { return type{hypot(x, y)}; }) + .def( + "invcbrt", +[](const type& t) { return type{invcbrt(t)}; }) + .def( + "invsqrt", +[](const type& t) { return type{invsqrt(t)}; }) + .def( + "log", +[](const type& t) { return type{log(t)}; }) + .def( + "log2", +[](const type& t) { return type{log2(t)}; }) + .def( + "log10", +[](const type& t) { return type{log10(t)}; }) + .def( + "max", +[](const type& t) { return Real{max(t)}; }) + .def( + "min", +[](const type& t) { return Real{min(t)}; }) + .def( + "pow", + +[](const type& base, double exp) { return type{pow(base, exp)}; }) + .def( + "sin", +[](const type& t) { return type{sin(t)}; }) + .def( + "sinh", +[](const type& t) { return type{sinh(t)}; }) + .def( + "sqrt", +[](const type& t) { return type{sqrt(t)}; }) + .def( + "tan", +[](const type& t) { return type{tan(t)}; }) + .def( + "tanh", +[](const type& t) { return type{tanh(t)}; }) + .def( + "__pow__", +[](const type& base, + const double exp) { return type{pow(base, exp)}; }) + .def( + "__add__", +[](const type& self, + const Real other) { return type{self + other}; }) + .def( + "__radd__", +[](const type& self, + const Real other) { return type{other + self}; }) + .def( + "__sub__", +[](const type& self, + const Real other) { return type{self - other}; }) + .def( + "__rsub__", +[](const type& self, + const Real other) { return type{other - self}; }) + .def( + "__mul__", +[](const type& self, + const Real other) { return type{self * other}; }) + .def( + "__rmul__", +[](const type& self, + const Real other) { return type{other * self}; }) + // Need __div__ for python 2 and __truediv__ for python 3. + .def( + "__div__", +[](const type& self, + const Real other) { return type{self / other}; }) + .def( + "__truediv__", +[](const type& self, + const Real other) { return type{self / other}; }) + .def( + "__rdiv__", +[](const type& self, + const Real other) { return type{other / self}; }) + .def( + "__rtruediv__", + +[](const type& self, const Real other) { + return type{other / self}; + }) + .def(py::self == py::self) + .def(py::self != py::self) + .def( + "__neg__", +[](const type& t) { return type{-t}; }); + } + //**************************************************************************** + +} // namespace py_bindings diff --git a/backend/src/Utilities/Math/Python/CMakeLists.txt b/backend/src/Utilities/Math/Python/CMakeLists.txt new file mode 100644 index 000000000..970d41e9c --- /dev/null +++ b/backend/src/Utilities/Math/Python/CMakeLists.txt @@ -0,0 +1,28 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY "PyArrays") + +elastica_python_add_module( + Arrays + LIBRARY_NAME ${LIBRARY} + SOURCES + Bindings.cpp + BlazeVector.cpp + # BlazeSubvector.cpp + # BlazeIndexVector.cpp + # BlazeIndexSubVector.cpp + BlazeMatrix.cpp + # BlazeSubmatrix.cpp + BlazeTensor.cpp + # BlazeSubtensor.cpp +) + +elastica_python_link_libraries( + ${LIBRARY} + PRIVATE + Blaze + BlazeTensor + pybind11::module + Utilities +) diff --git a/backend/src/Utilities/Math/Python/SliceHelpers.hpp b/backend/src/Utilities/Math/Python/SliceHelpers.hpp new file mode 100644 index 000000000..a810a9938 --- /dev/null +++ b/backend/src/Utilities/Math/Python/SliceHelpers.hpp @@ -0,0 +1,108 @@ +//****************************************************************************** +// Includes +//****************************************************************************** + +#pragma once +// +// #include "Utilities/MakeString.hpp" +// +#include +#include +#include +#include +#include +#include +// +#include +// +#include +#include +#include +#include + +namespace py_bindings { + + struct SliceInfo { + //! Start of slice + std::size_t start; + //! Length of slice + std::size_t slicelength; + }; + + //**************************************************************************** + /*!\brief Checks validitity of slice along Axis, raises error if invalid + * \ingroup python_bindings + */ + template + auto check_slice(const std::size_t limit, pybind11::slice const slice) { + std::size_t start, stop, step, slicelength; + if (!slice.compute(limit, &start, &stop, &step, &slicelength)) + throw pybind11::error_already_set(); + if (step != 1) + throw std::runtime_error(std::string( + "step !=1 unsupported along axis " << Axis)); + + return SliceInfo{start, slicelength}; + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Takes 1D slices of an array + * \ingroup python_bindings + */ + template + auto array_slice(T& t, pybind11::slice slice) { + constexpr std::size_t axis = 0UL; + auto slice_info = check_slice(t.size(), std::move(slice)); + return ::blaze::subvector(t, slice_info.start, slice_info.slicelength); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Takes 2D slices of an array + * \ingroup python_bindings + */ + template + auto array_slice(T& t, std::tuple slices) { + constexpr std::size_t row_slice = 0UL; + constexpr std::size_t col_slice = 1UL; + auto slice_info = std::make_tuple( + check_slice(t.rows(), + std::get(std::move(slices))), + check_slice(t.columns(), + std::get(std::move(slices)))); + return ::blaze::submatrix(t, std::get(slice_info).start, + std::get(slice_info).start, + std::get(slice_info).slicelength, + std::get(slice_info).slicelength); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Takes 3D slices of an array + * \ingroup python_bindings + */ + template + auto array_slice( + T& t, + std::tuple slices) { + constexpr std::size_t page_slice = 0UL; + constexpr std::size_t row_slice = 1UL; + constexpr std::size_t col_slice = 2UL; + auto slice_info = std::make_tuple( + check_slice(t.pages(), + std::get(std::move(slices))), + check_slice(t.rows(), + std::get(std::move(slices))), + check_slice(t.columns(), + std::get(std::move(slices)))); + return ::blaze::subtensor(t, std::get(slice_info).start, + std::get(slice_info).start, + std::get(slice_info).start, + std::get(slice_info).slicelength, + std::get(slice_info).slicelength, + std::get(slice_info).slicelength); + } + //**************************************************************************** + +} // namespace py_bindings From fb7f32d55abb1df98980e7988d87811913ba6726 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 26 Jun 2024 09:51:58 -0500 Subject: [PATCH 09/62] include SO3 backend implementation --- .../Systems/States/Expressions/backends.hpp | 14 ++ .../Expressions/backends/CMakeLists.txt | 27 +++ .../Expressions/backends/Declarations.hpp | 27 +++ .../States/Expressions/backends/README.md | 5 + .../States/Expressions/backends/blaze.hpp | 11 ++ .../Expressions/backends/blaze/CMakeLists.txt | 39 ++++ .../Expressions/backends/blaze/Resize.hpp | 46 +++++ .../backends/blaze/SO3PrimitiveAddAssign.hpp | 38 ++++ .../SO3PrimitiveAddAssign/BaseTemplate.hpp | 32 ++++ .../blaze/SO3PrimitiveAddAssign/SIMD.hpp | 181 ++++++++++++++++++ .../blaze/SO3PrimitiveAddAssign/Scalar.hpp | 108 +++++++++++ .../backends/blaze/SO3PrimitiveAssign.hpp | 58 ++++++ .../Expressions/backends/blaze/Size.hpp | 38 ++++ 13 files changed, 624 insertions(+) create mode 100644 backend/src/Systems/States/Expressions/backends.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/CMakeLists.txt create mode 100644 backend/src/Systems/States/Expressions/backends/Declarations.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/README.md create mode 100644 backend/src/Systems/States/Expressions/backends/blaze.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/CMakeLists.txt create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/Resize.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/SIMD.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp create mode 100644 backend/src/Systems/States/Expressions/backends/blaze/Size.hpp diff --git a/backend/src/Systems/States/Expressions/backends.hpp b/backend/src/Systems/States/Expressions/backends.hpp new file mode 100644 index 000000000..97f39d78e --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends.hpp @@ -0,0 +1,14 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +/// +#include "Systems/States/Expressions/backends/Declarations.hpp" +/// + +#define ELASTICA_USE_BLAZE 1 // todo : from CMAKE + +#if defined(ELASTICA_USE_BLAZE) +#include "Systems/States/Expressions/backends/blaze.hpp" +#endif // ELASTICA_USE_BLAZE diff --git a/backend/src/Systems/States/Expressions/backends/CMakeLists.txt b/backend/src/Systems/States/Expressions/backends/CMakeLists.txt new file mode 100644 index 000000000..c50e31671 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/CMakeLists.txt @@ -0,0 +1,27 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY StateExpressionBackends) + +add_elastica_library(${LIBRARY} INTERFACE) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + blaze.hpp + Declarations.hpp +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + blaze.hpp + Declarations.hpp +) + +add_subdirectory(blaze) + +target_link_libraries( + ${LIBRARY} + INTERFACE + SO3PrimitiveAddAssignKernel +) diff --git a/backend/src/Systems/States/Expressions/backends/Declarations.hpp b/backend/src/Systems/States/Expressions/backends/Declarations.hpp new file mode 100644 index 000000000..1a545acd0 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/Declarations.hpp @@ -0,0 +1,27 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include // std::size_t + +namespace elastica { + + namespace states { + + // assign version doesnt care about vectorization + template + auto SO3_primitive_assign(LHSMatrix& lhs, RHSMatrix const& rhs_matrix, + RHSVector const& rhs_vector) noexcept -> void; + // add assign version doesnt care about vectorization + template + auto SO3_primitive_add_assign(LHSMatrix& lhs, + RHSVector const& rhs_vector) noexcept -> void; + template + inline auto size_backend(Type const&) noexcept -> std::size_t; + template + inline auto resize_backend(Type&, std::size_t) -> void; + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/README.md b/backend/src/Systems/States/Expressions/backends/README.md new file mode 100644 index 000000000..6d2d95092 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/README.md @@ -0,0 +1,5 @@ +## backends + +To interface the requirements of SO3 and SE3 templates with any custom type that the user provides, we need a couple of +interface functions to connect the frontend (state templates) to the backend data type. This folder contains all such +interface functions for commonly used elastica types. diff --git a/backend/src/Systems/States/Expressions/backends/blaze.hpp b/backend/src/Systems/States/Expressions/backends/blaze.hpp new file mode 100644 index 000000000..221763f9a --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze.hpp @@ -0,0 +1,11 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include "Systems/States/Expressions/backends/Declarations.hpp" +/// +#include "Systems/States/Expressions/backends/blaze/Resize.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp" +#include "Systems/States/Expressions/backends/blaze/Size.hpp" diff --git a/backend/src/Systems/States/Expressions/backends/blaze/CMakeLists.txt b/backend/src/Systems/States/Expressions/backends/blaze/CMakeLists.txt new file mode 100644 index 000000000..b0abc964c --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/CMakeLists.txt @@ -0,0 +1,39 @@ +# Distributed under the MIT License. See LICENSE.txt for details. + +set(LIBRARY SO3PrimitiveAddAssignKernel) + +add_elastica_library(${LIBRARY} INTERFACE) + +elastica_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/elastica + SO3PrimitiveAddAssign/BaseTemplate.hpp + SO3PrimitiveAddAssign/Scalar.hpp + SO3PrimitiveAddAssign/SIMD.hpp + Resize.hpp + Size.hpp + SO3PrimitiveAddAssign.hpp + SO3PrimitiveAssign.hpp +) + +elastica_target_sources( + ${LIBRARY} + INTERFACE + SO3PrimitiveAddAssign/BaseTemplate.hpp + SO3PrimitiveAddAssign/Scalar.hpp + SO3PrimitiveAddAssign/SIMD.hpp + Resize.hpp + Size.hpp + SO3PrimitiveAddAssign.hpp + SO3PrimitiveAssign.hpp +) + + +target_link_libraries( + ${LIBRARY} + INTERFACE + Utilities + Blaze + Frames + ModuleSettings +) diff --git a/backend/src/Systems/States/Expressions/backends/blaze/Resize.hpp b/backend/src/Systems/States/Expressions/backends/blaze/Resize.hpp new file mode 100644 index 000000000..da2a7bc50 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/Resize.hpp @@ -0,0 +1,46 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include // std::size_t + +#include +#include + +// +// this comes from the simulator module, which seems like an anti-pattern +// should be in Systems instead +#include "Simulator/Frames.hpp" +// +#include "Systems/States/Expressions/backends/Declarations.hpp" + +namespace elastica { + + namespace states { + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + // the parts below are repeated from the Block module, but they are almost + // always never used in the states module since the size is always expected + // to be + template + inline auto resize_backend(blaze::DynamicTensor& data, + std::size_t new_size) -> void { + return data.resize(Frames::Dimension, Frames::Dimension, new_size, true); + } + + template // Type tag + inline auto resize_backend(blaze::DynamicMatrix& data, + std::size_t new_size) -> void { + return data.resize(Frames::Dimension, new_size, true); + } + /*! \endcond */ + //************************************************************************** + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign.hpp b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign.hpp new file mode 100644 index 000000000..35d6d3fa7 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign.hpp @@ -0,0 +1,38 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +// +#include "Systems/States/Expressions/backends/Declarations.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/SIMD.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp" +#include "Systems/States/Expressions/backends/blaze/Size.hpp" + +namespace elastica { + + namespace states { + + template + auto SO3_primitive_add_assign(blaze::DynamicTensor& lhs_matrix_batch, + RHSVectorBatch const& rhs_vector) noexcept + -> void { + /* + * This version creates memory every time + auto temp_to_synch(lhs_matrix_batch); // TODO: avoid making new memory + // (but not sure how) + // The rotation computation cannot be computed asynchronously during + // add_assign. + SO3_primitive_assign(lhs_matrix_batch, temp_to_synch, rhs_vector); + */ + + detail::SO3AddAssign()>:: + apply(lhs_matrix_batch, rhs_vector); + } + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp new file mode 100644 index 000000000..91df12873 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "Configuration/Kernels.hpp" +#include "ModuleSettings/Vectorization.hpp" + +namespace elastica { + + namespace states { + + namespace detail { + + enum class SO3AddAssignKind { scalar, simd }; + + template + struct SO3AddAssign; + + template + constexpr auto backend_choice() -> BackendEnumKind; + + template <> + constexpr auto backend_choice() -> SO3AddAssignKind { + return ELASTICA_USE_VECTORIZATION + ? SO3AddAssignKind:: + ELASTICA_COSSERATROD_LIB_SO3_ADDITION_BACKEND + : SO3AddAssignKind::scalar; + } + + } // namespace detail + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/SIMD.hpp b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/SIMD.hpp new file mode 100644 index 000000000..bb8c04ed5 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/SIMD.hpp @@ -0,0 +1,181 @@ +#pragma once + +#include +#include +// +#include +#include // for padding +// +#include +#include +// +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp" +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp" +#include "Systems/States/Expressions/backends/blaze/Size.hpp" +// +#include "ErrorHandling/Assert.hpp" +// +#include "Utilities/Unroll.hpp" + +namespace elastica { + + namespace states { + + namespace detail { + + template <> + struct SO3AddAssign { + template + static auto apply(blaze::DynamicTensor& lhs_matrix_batch, + RHSVectorBatch const& rhs_vector_batch) noexcept + -> void { + using ::elastica::states::size_backend; + const std::size_t n_dofs = size_backend(lhs_matrix_batch); + constexpr bool remainder(!::blaze::usePadding || + !blaze::IsPadded_v); + constexpr std::size_t SIMDSIZE(blaze::SIMDTrait::size); + + const std::size_t i_dof_pos( + remainder ? blaze::prevMultiple(n_dofs, SIMDSIZE) : n_dofs); + ELASTICA_ASSERT(i_dof_pos <= n_dofs, "Invalid end calculation"); + + // first do loops of size simd width + + // begin has (row, page) as the syntax + typename blaze::DynamicTensor::Iterator q_it[9UL] = { + lhs_matrix_batch.begin(0UL, 0UL), + lhs_matrix_batch.begin(1UL, 0UL), + lhs_matrix_batch.begin(2UL, 0UL), + lhs_matrix_batch.begin(0UL, 1UL), + lhs_matrix_batch.begin(1UL, 1UL), + lhs_matrix_batch.begin(2UL, 1UL), + lhs_matrix_batch.begin(0UL, 2UL), + lhs_matrix_batch.begin(1UL, 2UL), + lhs_matrix_batch.begin(2UL, 2UL)}; + + typename RHSVectorBatch::ConstIterator u_it[3UL] = { + rhs_vector_batch.begin(0UL), rhs_vector_batch.begin(1UL), + rhs_vector_batch.begin(2UL)}; + // + using SIMDType = typename blaze::DynamicTensor::SIMDType; + + // decide how much to unroll here + // unrolling more probably wont help because of register pressure + SIMDType q[9UL]; + SIMDType u[3UL]; + SIMDType c_alpha, alpha, beta; + +#if !(BLAZE_SVML_MODE || BLAZE_SLEEF_MODE) + // a temporary aligned register to pipe contents into + alignas(::blaze::AlignmentOf_v) T temp[SIMDSIZE]; +#endif + + std::size_t i_dof = 0UL; + + for (; i_dof < i_dof_pos; i_dof += SIMDSIZE) { + { + // 1. load data into registers + + // 1.1 load u + UNROLL_LOOP(3UL) + for (std::size_t dim = 0UL; dim < 3UL; ++dim) { + u[dim] = u_it[dim].load(); + } + + // 1.2 load q + UNROLL_LOOP(9UL) + for (std::size_t dim = 0UL; dim < 9UL; ++dim) { + q[dim] = q_it[dim].load(); + } + } + + { + // 2. compute the angle of rotation + // blaze dispatches to std + beta = ::blaze::sqrt(u[0UL] * u[0UL] + u[1UL] * u[1UL] + + u[2UL] * u[2UL]); + + // 3. normalize the axis of rotation + /*Clang-Tidy: Use range-based for loop instead*/ + UNROLL_LOOP(3UL) + for (std::size_t dim = 0UL; dim < 3UL; ++dim) { /* NOLINT */ + // bad access pattern + // TODO : refactor magic number + auto arg = (beta + ::blaze::set(T(1e-14))); + u[dim] /= arg; + } + + // compute trigonometric entities +#if BLAZE_SVML_MODE || BLAZE_SLEEF_MODE + c_alpha = ::blaze::cos(beta); + beta = ::blaze::sin(beta); +#else + // no intrinsics, so we unpack the SIMD vector, compute sin and + // cos on each via regular math functions and then repack it +#pragma unroll SIMDSIZE + for (std::size_t simd_idx = 0UL; simd_idx < SIMDSIZE; ++simd_idx) + temp[simd_idx] = ::blaze::cos(beta[simd_idx]); + + c_alpha = ::blaze::loada(temp); + +#pragma unroll SIMDSIZE + for (std::size_t simd_idx = 0UL; simd_idx < SIMDSIZE; ++simd_idx) + temp[simd_idx] = ::blaze::sin(beta[simd_idx]); + + beta = ::blaze::loada(temp); +#endif + alpha = ::blaze::set(T(1.0)) - c_alpha; + } + + { + // 4. compute the rotated batch + /*Clang-Tidy: Use range-based for loop instead*/ + UNROLL_LOOP(3UL) + for (std::size_t j_dim = 0UL; j_dim < 3UL; ++j_dim) { /* NOLINT */ + // Force evaluate it + const auto com = (q[j_dim] * u[0UL] + q[3UL + j_dim] * u[1UL] + + q[6UL + j_dim] * u[2UL]) + .eval(); + + q_it[0UL + j_dim].stream( + c_alpha * q[j_dim] + alpha * u[0UL] * com + + beta * (q[3UL + j_dim] * u[2UL] - q[6UL + j_dim] * u[1UL])); + + q_it[3UL + j_dim].stream( + c_alpha * q[3UL + j_dim] + alpha * u[1UL] * com + + beta * (q[6UL + j_dim] * u[0UL] - q[j_dim] * u[2UL])); + + q_it[6UL + j_dim].stream( + c_alpha * q[6UL + j_dim] + alpha * u[2UL] * com + + beta * (q[j_dim] * u[1UL] - q[3UL + j_dim] * u[0UL])); + + } // jdim + } + + { + // 5. advance all the iterators to the next SIMD lane + /*Clang-Tidy: Use range-based for loop instead*/ + UNROLL_LOOP(3UL) + for (std::size_t dim = 0UL; dim < 3UL; ++dim) /* NOLINT */ + u_it[dim] += SIMDSIZE; + + UNROLL_LOOP(9UL) + for (std::size_t dim = 0UL; dim < 9UL; ++dim) /* NOLINT */ + q_it[dim] += SIMDSIZE; + } + + } // i_dof + + // then do the last loops, peeling them off to serial scalar + // implementation UNTESTED + if (remainder) + SO3AddAssign::internal_apply( + lhs_matrix_batch, rhs_vector_batch, i_dof, n_dofs); + } + }; + + } // namespace detail + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp new file mode 100644 index 000000000..6d9e76223 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/Scalar.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include +#include +#include + +#include "Systems/States/Expressions/backends/blaze/SO3PrimitiveAddAssign/BaseTemplate.hpp" +#include "Systems/States/Expressions/backends/blaze/Size.hpp" +// +#include "Utilities/Unroll.hpp" + +namespace elastica { + + namespace states { + + namespace detail { + + template <> + struct SO3AddAssign { + template + static auto internal_apply(blaze::DynamicTensor& lhs_matrix_batch, + RHSVectorBatch const& rhs_vector_batch, + const std::size_t index_start, + const std::size_t index_fin) noexcept + -> void { + // 0. allocate memory for rolling window + T q[9UL]; + /* + * structure is + * q0 q3 q6 + * q1 q4 q7 + * q2 q5 q8 + */ + T u[3UL]; + // T theta, alpha, beta; // doesnt matter the extra memory + + for (std::size_t i_dof = index_start; i_dof < index_fin; ++i_dof) { + // 1. load data into registers + UNROLL_LOOP(3UL) + for (std::size_t dim = 0UL; dim < 3UL; ++dim) { + // bad access pattern + u[dim] = rhs_vector_batch(dim, i_dof); + } + + UNROLL_LOOP(3UL) // Needs to be timed + for (std::size_t k_dim = 0UL; k_dim < 3UL; ++k_dim) { + UNROLL_LOOP(3UL) + for (std::size_t j_dim = 0UL; j_dim < 3UL; ++j_dim) { + // bad access pattern + // convention according to SO3PrimitiveAssign + q[3UL * k_dim + j_dim] = lhs_matrix_batch(k_dim, j_dim, i_dof); + } + } + + // Now work only with registers + + // 2. compute the angle of rotation + // blaze dispatches to std + const T theta = ::std::sqrt(u[0UL] * u[0UL] + u[1UL] * u[1UL] + + u[2UL] * u[2UL]); + const T c_alpha = ::std::cos(theta); + const T alpha = T(1.0) - c_alpha; + const T beta = ::std::sin(theta); + + // 3. normalize the axis of rotation + /*Clang-Tidy: Use range-based for loop instead*/ + UNROLL_LOOP(3UL) + for (std::size_t dim = 0UL; dim < 3UL; ++dim) { /* NOLINT */ + // bad access pattern + u[dim] /= (theta + 1e-14); // TODO : refactor magic number + } + + /*Clang-Tidy: Use range-based for loop instead*/ + UNROLL_LOOP(3UL) + for (std::size_t j_dim = 0UL; j_dim < 3UL; ++j_dim) { /* NOLINT */ + const T com = q[j_dim] * u[0UL] + q[3UL + j_dim] * u[1UL] + + q[6UL + j_dim] * u[2UL]; + + lhs_matrix_batch(0UL, j_dim, i_dof) = + c_alpha * q[j_dim] + alpha * u[0UL] * com + + beta * (q[3UL + j_dim] * u[2UL] - q[6UL + j_dim] * u[1UL]); + + lhs_matrix_batch(1UL, j_dim, i_dof) = + c_alpha * q[3UL + j_dim] + alpha * u[1UL] * com + + beta * (q[6UL + j_dim] * u[0UL] - q[j_dim] * u[2UL]); + + lhs_matrix_batch(2UL, j_dim, i_dof) = + c_alpha * q[6UL + j_dim] + alpha * u[2UL] * com + + beta * (q[j_dim] * u[1UL] - q[3UL + j_dim] * u[0UL]); + } // j_dim + } // i_dof + } + + template + static inline auto apply( + blaze::DynamicTensor& lhs_matrix_batch, + RHSVectorBatch const& rhs_vector_batch) noexcept -> void { + using ::elastica::states::size_backend; + internal_apply(lhs_matrix_batch, rhs_vector_batch, 0UL, + size_backend(lhs_matrix_batch)); + } + }; + + } // namespace detail + + } // namespace states + +} // namespace elastica diff --git a/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp new file mode 100644 index 000000000..a103f1168 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/SO3PrimitiveAssign.hpp @@ -0,0 +1,58 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +#include "Systems/States/Expressions/backends/Declarations.hpp" +#include "Utilities/Math/BlazeDetail/BlazeRotation.hpp" + +namespace elastica { + + namespace states { + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + // clang-format off + template + auto SO3_primitive_assign(blaze::DynamicTensor& lhs_matrix_batch, + blaze::DynamicTensor const& rhs_matrix_batch, + RHSVectorBatch const& vector_batch) noexcept // vector bacth includes (omega*dt). Must include minus sign here + -> void { +// using ::elastica::exp_batch; +// exp_batch(lhs_matrix_batch, -vector_batch); +// lhs_matrix_batch = lhs_matrix_batch * rhs_matrix_batch; + // TODO: Double check. Place the function externally. Double check the sign + auto theta = blaze::sqrt(blaze::sum(vector_batch % vector_batch)); + auto alpha = 1.0 - blaze::cos(theta); + auto beta = blaze::sin(theta); + auto u0 = blaze::row(vector_batch, 0UL) / (theta+1e-14); + auto u1 = blaze::row(vector_batch, 1UL) / (theta+1e-14); + auto u2 = blaze::row(vector_batch, 2UL) / (theta+1e-14); + auto q0 = blaze::row(blaze::pageslice(rhs_matrix_batch, 0UL), 0UL); + auto q1 = blaze::row(blaze::pageslice(rhs_matrix_batch, 0UL), 1UL); + auto q2 = blaze::row(blaze::pageslice(rhs_matrix_batch, 0UL), 2UL); + auto q3 = blaze::row(blaze::pageslice(rhs_matrix_batch, 1UL), 0UL); + auto q4 = blaze::row(blaze::pageslice(rhs_matrix_batch, 1UL), 1UL); + auto q5 = blaze::row(blaze::pageslice(rhs_matrix_batch, 1UL), 2UL); + auto q6 = blaze::row(blaze::pageslice(rhs_matrix_batch, 2UL), 0UL); + auto q7 = blaze::row(blaze::pageslice(rhs_matrix_batch, 2UL), 1UL); + auto q8 = blaze::row(blaze::pageslice(rhs_matrix_batch, 2UL), 2UL); + // TODO: maybe replace blaze::pow(x,2) to x % x + blaze::row(blaze::pageslice(lhs_matrix_batch, 0UL), 0UL) = alpha*((-blaze::pow(u1, 2) - blaze::pow(u2, 2))*q0 + q3*u0*u1 + q6*u0*u2) + beta*( q3*u2 - q6*u1) + q0; + blaze::row(blaze::pageslice(lhs_matrix_batch, 0UL), 1UL) = alpha*((-blaze::pow(u1, 2) - blaze::pow(u2, 2))*q1 + q4*u0*u1 + q7*u0*u2) + beta*( q4*u2 - q7*u1) + q1; + blaze::row(blaze::pageslice(lhs_matrix_batch, 0UL), 2UL) = alpha*((-blaze::pow(u1, 2) - blaze::pow(u2, 2))*q2 + q5*u0*u1 + q8*u0*u2) + beta*( q5*u2 - q8*u1) + q2; + blaze::row(blaze::pageslice(lhs_matrix_batch, 1UL), 0UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u2, 2))*q3 + q0*u0*u1 + q6*u1*u2) + beta*(-q0*u2 + q6*u0) + q3; + blaze::row(blaze::pageslice(lhs_matrix_batch, 1UL), 1UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u2, 2))*q4 + q1*u0*u1 + q7*u1*u2) + beta*(-q1*u2 + q7*u0) + q4; + blaze::row(blaze::pageslice(lhs_matrix_batch, 1UL), 2UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u2, 2))*q5 + q2*u0*u1 + q8*u1*u2) + beta*(-q2*u2 + q8*u0) + q5; + blaze::row(blaze::pageslice(lhs_matrix_batch, 2UL), 0UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u1, 2))*q6 + q0*u0*u2 + q3*u1*u2) + beta*( q0*u1 - q3*u0) + q6; + blaze::row(blaze::pageslice(lhs_matrix_batch, 2UL), 1UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u1, 2))*q7 + q1*u0*u2 + q4*u1*u2) + beta*( q1*u1 - q4*u0) + q7; + blaze::row(blaze::pageslice(lhs_matrix_batch, 2UL), 2UL) = alpha*((-blaze::pow(u0, 2) - blaze::pow(u1, 2))*q8 + q2*u0*u2 + q5*u1*u2) + beta*( q2*u1 - q5*u0) + q8; + } + // clang-format on + + /*! \endcond */ + //************************************************************************** + } + +} diff --git a/backend/src/Systems/States/Expressions/backends/blaze/Size.hpp b/backend/src/Systems/States/Expressions/backends/blaze/Size.hpp new file mode 100644 index 000000000..a9b930039 --- /dev/null +++ b/backend/src/Systems/States/Expressions/backends/blaze/Size.hpp @@ -0,0 +1,38 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +#include +#include // std::size_t + +#include "Systems/States/Expressions/backends/Declarations.hpp" + +namespace elastica { + + namespace states { + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + template + inline auto size_backend( + blaze::DynamicTensor const& tensor_batch) noexcept -> std::size_t { + return tensor_batch.columns(); + } + + template // Type tag + inline auto size_backend( + blaze::DynamicMatrix const& matrix_batch) noexcept + -> std::size_t { + return matrix_batch.columns(); + } + /*! \endcond */ + //************************************************************************** + + } // namespace states + +} // namespace elastica From f66d0502bda668c3ed70e58bb6d1b3759ea722fe Mon Sep 17 00:00:00 2001 From: Ankith Date: Fri, 5 Jul 2024 14:46:17 +0530 Subject: [PATCH 10/62] Fix Utilities/Math/Python issues and make it build --- backend/meson.build | 2 +- backend/src/Utilities/Math/Python/BlazeMatrix.cpp | 6 +++--- backend/src/Utilities/Math/Python/BlazeTensor.cpp | 14 +++++++------- backend/src/Utilities/Math/Python/BlazeVector.cpp | 6 +++--- backend/src/Utilities/Math/Python/SliceHelpers.hpp | 2 +- backend/src/meson.build | 14 ++++++++++++++ 6 files changed, 29 insertions(+), 15 deletions(-) diff --git a/backend/meson.build b/backend/meson.build index e9a1d4650..4c5cce8f2 100644 --- a/backend/meson.build +++ b/backend/meson.build @@ -2,7 +2,7 @@ project( 'elasticapp', ['cpp'], version: '0.0.3', - default_options: ['cpp_std=c++17'], + default_options: ['cpp_std=c++14'], ) package = 'elasticapp' diff --git a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp index 9c3fb654e..56f55a110 100644 --- a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp +++ b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp @@ -90,9 +90,9 @@ namespace py_bindings { return array_slice(t, std::move(slice)); }) // Need __str__ for converting to string/printing - .def( - "__str__", - +[](const type& self) { return std::string(MakeString{} << self); }) + // .def( + // "__str__", + // +[](const type& self) { return std::string(MakeString{} << self); }) .def( "__setitem__", +[](type& self, const std::tuple& x, diff --git a/backend/src/Utilities/Math/Python/BlazeTensor.cpp b/backend/src/Utilities/Math/Python/BlazeTensor.cpp index 57dc73a4d..156f7f891 100644 --- a/backend/src/Utilities/Math/Python/BlazeTensor.cpp +++ b/backend/src/Utilities/Math/Python/BlazeTensor.cpp @@ -83,7 +83,7 @@ namespace py_bindings { +[](const type& self, const std::tuple& x) { // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), - std::get<2>(x)); + // std::get<2>(x)); return self(std::get<0>(x), std::get<1>(x), std::get<2>(x)); }) .def( @@ -97,14 +97,14 @@ namespace py_bindings { const std::tuple& x, const Real val) { // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), - std::get<2>(x)); + // std::get<2>(x)); self(std::get<0>(x), std::get<1>(x), std::get<2>(x)) = val; - }) - // Need __str__ for converting to string/printing - .def( - "__str__", +[](const type& self) { - return std::string(MakeString{} << self); }); + // Need __str__ for converting to string/printing + // .def( + // "__str__", +[](const type& self) { + // return std::string(MakeString{} << self); + // }); } //**************************************************************************** diff --git a/backend/src/Utilities/Math/Python/BlazeVector.cpp b/backend/src/Utilities/Math/Python/BlazeVector.cpp index adedad39f..86239102f 100644 --- a/backend/src/Utilities/Math/Python/BlazeVector.cpp +++ b/backend/src/Utilities/Math/Python/BlazeVector.cpp @@ -112,11 +112,11 @@ namespace py_bindings { // Need __str__ for converting to string py_vector - .def("__str__") + // .def("__str__") // repr allows you to output the object in an interactive python // terminal using obj to get the "string REPResenting the object". - .def("__repr__") - .def(py::self += py::self) + // .def("__repr__") + // .def(py::self += py::self) // Need to do math explicitly converting to DataVector because we don't // want to represent all the possible expression template types .def( diff --git a/backend/src/Utilities/Math/Python/SliceHelpers.hpp b/backend/src/Utilities/Math/Python/SliceHelpers.hpp index a810a9938..cc6a5b5fa 100644 --- a/backend/src/Utilities/Math/Python/SliceHelpers.hpp +++ b/backend/src/Utilities/Math/Python/SliceHelpers.hpp @@ -40,7 +40,7 @@ namespace py_bindings { throw pybind11::error_already_set(); if (step != 1) throw std::runtime_error(std::string( - "step !=1 unsupported along axis " << Axis)); + "step !=1 unsupported along axis ") + std::to_string(Axis)); return SliceInfo{start, slicelength}; } diff --git a/backend/src/meson.build b/backend/src/meson.build index 690dd2722..623d13960 100644 --- a/backend/src/meson.build +++ b/backend/src/meson.build @@ -15,3 +15,17 @@ _linalg = py.extension_module( install: true, subdir: package, ) + +_PyArrays = py.extension_module( + '_PyArrays', + sources: [ + 'Utilities/Math/Python/Bindings.cpp', + 'Utilities/Math/Python/BlazeVector.cpp', + 'Utilities/Math/Python/BlazeMatrix.cpp', + 'Utilities/Math/Python/BlazeTensor.cpp', + ], + dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], + link_language: 'cpp', + install: true, + subdir: package, +) From 3737d5c969e4ee5b82f362dbac9bf9f036d5a188 Mon Sep 17 00:00:00 2001 From: Ankith Date: Sun, 7 Jul 2024 18:03:50 +0530 Subject: [PATCH 11/62] Add elasticapp documentation --- backend/README.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 backend/README.md diff --git a/backend/README.md b/backend/README.md new file mode 100644 index 000000000..60b6c39ec --- /dev/null +++ b/backend/README.md @@ -0,0 +1,33 @@ +# Elasticapp backend for PyElastica + +This file serves as the documentation for the `elasticapp` backend. + +## Installation + +In the root of the PyElastica source tree, run the following command + +``` +pip install ./backend +``` + +This command will take care of installation of all build-time dependencies, compilation of C++ source files and install the a python package called `elasticapp`. + +## Testing + +Make sure you have `pytest` installed. In the root of the PyElastica source tree, run the following command + +``` +pytest backend/tests +``` + +## Benchmarking + +Standalone scripts for benchmarking purposes are available in `backend/benchmarking` folder. + +### Benchmarking `matmul` + +For benchmarking various `matmul` implementations, run + +``` +python3 backend/benchmarking/matmul.py +``` From 4e7a04b85c0270a8219af946fde70a71d1f39258 Mon Sep 17 00:00:00 2001 From: Ankith Date: Sun, 7 Jul 2024 22:24:25 +0530 Subject: [PATCH 12/62] port linalg to use Utils/Math/Python objects - Updated benchmarking and test scripts accordingly --- backend/benchmarking/matmul.py | 40 ++- backend/pyproject.toml | 3 + backend/src/_linalg.cpp | 315 +++--------------------- backend/src/_linalg_numpy.cpp | 376 +++++++++++++++++++++++++++++ backend/src/meson.build | 9 + backend/tests/test_linalg.py | 34 ++- backend/tests/test_linalg_numpy.py | 113 +++++++++ 7 files changed, 581 insertions(+), 309 deletions(-) create mode 100644 backend/src/_linalg_numpy.cpp create mode 100644 backend/tests/test_linalg_numpy.py diff --git a/backend/benchmarking/matmul.py b/backend/benchmarking/matmul.py index 04e4a2b18..87f439767 100644 --- a/backend/benchmarking/matmul.py +++ b/backend/benchmarking/matmul.py @@ -1,4 +1,10 @@ -from elasticapp._linalg import batch_matmul_naive, batch_matmul_blaze, batch_matmul +from elasticapp._linalg_numpy import ( + batch_matmul_naive, + batch_matmul_blaze, + batch_matmul, +) +from elasticapp._linalg import batch_matmul as batch_matmul_final +from elasticapp._PyArrays import Tensor from elastica._linalg import _batch_matmul import numpy import time @@ -16,12 +22,15 @@ def benchmark_batchsize(funcs: list, batches: list[int], num_iterations: int = 1 random_b = numpy.random.random((3, 3, batch_size)) ret[batch_size] = {} - for func in funcs: + for func_name, func, func_wrap in funcs: + random_a_w = func_wrap(random_a) if func_wrap else random_a + random_b_w = func_wrap(random_b) if func_wrap else random_b + start = time.perf_counter() for _ in range(num_iterations): - func(random_a, random_b) + func(random_a_w, random_b_w) - ret[batch_size][func.__name__] = ( + ret[batch_size][func_name] = ( time.perf_counter() - start ) / num_iterations @@ -29,16 +38,25 @@ def benchmark_batchsize(funcs: list, batches: list[int], num_iterations: int = 1 results = benchmark_batchsize( - [batch_matmul_naive, batch_matmul_blaze, batch_matmul, _batch_matmul], [2**i for i in range(14)] + [ + ("pyelastica", _batch_matmul, None), + ("elasticapp_naive", batch_matmul_naive, None), + ("elasticapp_blaze", batch_matmul_blaze, None), + ("elasticapp_blaze_copy", batch_matmul, None), + ("elasticapp_blaze_final", batch_matmul_final, Tensor), + ], + [2**i for i in range(14)], ) for size, data in results.items(): - pyelastica = data["_batch_matmul"] - elasticapp = data["batch_matmul_naive"] - elasticapp_blaze = data["batch_matmul_blaze"] - elasticapp_blaze_v2 = data["batch_matmul"] + pyelastica = data["pyelastica"] + elasticapp_naive = data["elasticapp_naive"] + elasticapp_blaze = data["elasticapp_blaze"] + elasticapp_blaze_copy = data["elasticapp_blaze_copy"] + elasticapp_blaze_final = data["elasticapp_blaze_final"] print(f"{size = }") print(f"{pyelastica = }") - print(f"{elasticapp = }, ratio: {elasticapp / pyelastica}") + print(f"{elasticapp_naive = }, ratio: {elasticapp_naive / pyelastica}") print(f"{elasticapp_blaze = }, ratio: {elasticapp_blaze / pyelastica}") - print(f"{elasticapp_blaze_v2 = }, ratio: {elasticapp_blaze_v2 / pyelastica}") + print(f"{elasticapp_blaze_copy = }, ratio: {elasticapp_blaze_copy / pyelastica}") + print(f"{elasticapp_blaze_final = }, ratio: {elasticapp_blaze_final / pyelastica}") print() diff --git a/backend/pyproject.toml b/backend/pyproject.toml index f10127b1b..35108c608 100644 --- a/backend/pyproject.toml +++ b/backend/pyproject.toml @@ -23,6 +23,9 @@ classifiers = [ requires = ["meson-python", "pybind11"] build-backend = "mesonpy" +[tool.meson-python.args] +setup = ['-Doptimization=3'] + # [tool.mypy] # files = "setup.py" # python_version = "3.7" diff --git a/backend/src/_linalg.cpp b/backend/src/_linalg.cpp index 651e29177..f726959d1 100644 --- a/backend/src/_linalg.cpp +++ b/backend/src/_linalg.cpp @@ -10,301 +10,50 @@ using blaze::DynamicMatrix; using blaze::DynamicTensor; using blaze::StaticMatrix; -/* Simple rewrite of elastica._linalg._batch_matmul */ -py::array_t batch_matmul_naive( - py::array_t first_matrix_collection, - py::array_t second_matrix_collection) -{ - auto s1 = first_matrix_collection.shape(0); - auto s2 = first_matrix_collection.shape(1); - auto max_k = first_matrix_collection.shape(2); - auto output_matrix = py::array_t{{s1, s2, max_k}}; - - auto a = first_matrix_collection.unchecked<3>(); - auto b = second_matrix_collection.unchecked<3>(); - auto c = output_matrix.mutable_unchecked<3>(); - for (py::ssize_t i = 0; i < 3; i++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - for (py::ssize_t m = 0; m < 3; m++) - { - for (py::ssize_t k = 0; k < max_k; k++) - { - c(i, m, k) += a(i, j, k) * b(j, m, k); - } - } - } - } - return output_matrix; -} +using ElasticaVector = ::blaze::DynamicVector>; +using ElasticaMatrix = ::blaze::DynamicMatrix>; +using ElasticaTensor = ::blaze::DynamicTensor; -/* blaze implementation of elastica._linalg._batch_matmul */ -py::array_t batch_matmul_blaze( - py::array_t first_matrix_collection, - py::array_t second_matrix_collection) +ElasticaMatrix difference_kernel(ElasticaMatrix &vector_batch) { - auto s1 = first_matrix_collection.shape(0); - auto s2 = first_matrix_collection.shape(1); - auto max_k = first_matrix_collection.shape(2); - auto output_matrix = py::array_t{{s1, s2, max_k}}; - - auto a = first_matrix_collection.unchecked<3>(); - auto b = second_matrix_collection.unchecked<3>(); - auto c = output_matrix.mutable_unchecked<3>(); - - StaticMatrix blaze_a; - StaticMatrix blaze_b; - StaticMatrix blaze_c; - for (py::ssize_t k = 0; k < max_k; k++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - blaze_a(i, j) = a(i, j, k); - blaze_b(i, j) = b(i, j, k); - } - } - blaze_c = blaze_a * blaze_b; - for (py::ssize_t i = 0; i < 3; i++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - c(i, j, k) = blaze_c(i, j); - } - } - } - - return output_matrix; -} - -py::array_t difference_kernel(py::array_t vector_batch) -{ - const auto v0 = vector_batch.shape(0); - const auto num_elems = vector_batch.shape(1); - - assert(v0 == 3UL); - auto output_arr = py::array_t{{v0, num_elems - 1}}; - - auto vector_batch_unchecked = vector_batch.unchecked<2>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); - - DynamicMatrix blaze_vector_batch(v0, num_elems); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); - } - } - auto blaze_output = elastica::difference_kernel(blaze_vector_batch); - for (py::ssize_t j = 0; j < num_elems - 1; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - output_arr_unchecked(i, j) = blaze_output(i, j); - } - } - - return output_arr; + return elastica::difference_kernel(vector_batch); } -py::array_t batch_matvec( - py::array_t matrix_collection, py::array_t vector_collection) +ElasticaMatrix batch_matvec( + ElasticaTensor &matrix_collection, ElasticaMatrix &vector_collection) { - const auto v0 = vector_collection.shape(0); - const auto num_elems = vector_collection.shape(1); - assert(v0 == 3UL); - - assert(matrix_collection.shape(0) == 3UL); - assert(matrix_collection.shape(1) == 3UL); - assert(matrix_collection.shape(2) == num_elems); - - auto output_arr = py::array_t{{v0, num_elems}}; - - auto matrix_collection_unchecked = matrix_collection.unchecked<3>(); - auto vector_collection_unchecked = vector_collection.unchecked<2>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); - - DynamicTensor blaze_matrix_collection(v0, v0, num_elems); - DynamicMatrix blaze_vector_collection(v0, num_elems); - for (py::ssize_t k = 0; k < num_elems; k++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_matrix_collection(i, j, k) = matrix_collection_unchecked(i, j, k); - } - blaze_vector_collection(j, k) = vector_collection_unchecked(j, k); - } - } - - DynamicMatrix blaze_output(v0, num_elems); - elastica::batch_matvec(blaze_output, blaze_matrix_collection, blaze_vector_collection); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - output_arr_unchecked(i, j) = blaze_output(i, j); - } - } - - return output_arr; + ElasticaMatrix blaze_output(matrix_collection.rows(), matrix_collection.columns()); + elastica::batch_matvec(blaze_output, matrix_collection, vector_collection); + return blaze_output; } -py::array_t batch_matmul( - py::array_t first_matrix_batch, py::array_t second_matrix_batch) +ElasticaTensor batch_matmul( + ElasticaTensor &first_matrix_batch, ElasticaTensor &second_matrix_batch) { - const auto m0 = first_matrix_batch.shape(0); - const auto m1 = first_matrix_batch.shape(1); - const auto num_elems = first_matrix_batch.shape(2); - assert(m0 == 3UL); - assert(m1 == 3UL); - - assert(second_matrix_batch.shape(0) == 3UL); - assert(second_matrix_batch.shape(1) == 3UL); - assert(second_matrix_batch.shape(2) == num_elems); - - auto output_arr = py::array_t{{m0, m1, num_elems}}; - - auto first_matrix_batch_unchecked = first_matrix_batch.unchecked<3>(); - auto second_matrix_batch_unchecked = second_matrix_batch.unchecked<3>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<3>(); - - DynamicTensor blaze_first_matrix_batch(m0, m1, num_elems); - DynamicTensor blaze_second_matrix_batch(m0, m1, num_elems); - for (py::ssize_t k = 0; k < num_elems; k++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_first_matrix_batch(i, j, k) = first_matrix_batch_unchecked(i, j, k); - blaze_second_matrix_batch(i, j, k) = second_matrix_batch_unchecked(i, j, k); - } - } - } - DynamicTensor blaze_output(m0, m1, num_elems); - elastica::batch_matmul(blaze_output, blaze_first_matrix_batch, blaze_second_matrix_batch); - for (py::ssize_t k = 0; k < num_elems; k++) - { - for (py::ssize_t j = 0; j < 3; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - output_arr_unchecked(i, j, k) = blaze_output(i, j, k); - } - } - } - return output_arr; + ElasticaTensor blaze_output(first_matrix_batch.pages(), first_matrix_batch.rows(), first_matrix_batch.columns()); + elastica::batch_matmul(blaze_output, first_matrix_batch, second_matrix_batch); + return blaze_output; } -py::array_t batch_cross( - py::array_t first_vector_batch, py::array_t second_vector_batch) +ElasticaMatrix batch_cross( + ElasticaMatrix &first_vector_batch, ElasticaMatrix &second_vector_batch) { - const auto v0 = first_vector_batch.shape(0); - const auto num_elems = first_vector_batch.shape(1); - assert(v0 == 3UL); - - assert(second_vector_batch.shape(0) == 3UL); - assert(second_vector_batch.shape(1) == num_elems); - - auto output_arr = py::array_t{{v0, num_elems}}; - - auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); - auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); - - DynamicMatrix blaze_first_vector_batch(v0, num_elems); - DynamicMatrix blaze_second_vector_batch(v0, num_elems); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); - blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); - } - } - - DynamicMatrix blaze_output(v0, num_elems); - elastica::batch_cross(blaze_output, blaze_first_vector_batch, blaze_second_vector_batch); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - output_arr_unchecked(i, j) = blaze_output(i, j); - } - } - - return output_arr; + ElasticaMatrix blaze_output(first_vector_batch.rows(), first_vector_batch.columns()); + elastica::batch_cross(blaze_output, first_vector_batch, second_vector_batch); + return blaze_output; } -py::array_t batch_dot( - py::array_t first_vector_batch, py::array_t second_vector_batch) +ElasticaVector batch_dot( + ElasticaMatrix &first_vector_batch, ElasticaMatrix &second_vector_batch) { - const auto v0 = first_vector_batch.shape(0); - const auto num_elems = first_vector_batch.shape(1); - assert(v0 == 3UL); - - assert(second_vector_batch.shape(0) == 3UL); - assert(second_vector_batch.shape(1) == num_elems); - - auto output_arr = py::array_t{num_elems}; - - auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); - auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); - - DynamicMatrix blaze_first_vector_batch(v0, num_elems); - DynamicMatrix blaze_second_vector_batch(v0, num_elems); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); - blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); - } - } - - auto blaze_output = elastica::batch_dot(blaze_first_vector_batch, blaze_second_vector_batch); - for (py::ssize_t j = 0; j < num_elems; j++) - { - output_arr_unchecked(j) = blaze_output[j]; - } - - return output_arr; + return elastica::batch_dot(first_vector_batch, second_vector_batch); } -py::array_t batch_norm(py::array_t vector_batch) +ElasticaVector batch_norm(ElasticaMatrix &vector_batch) { - const auto v0 = vector_batch.shape(0); - const auto num_elems = vector_batch.shape(1); - - assert(v0 == 3UL); - - auto output_arr = py::array_t{num_elems}; - - auto vector_batch_unchecked = vector_batch.unchecked<2>(); - auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); - - DynamicMatrix blaze_vector_batch(v0, num_elems); - for (py::ssize_t j = 0; j < num_elems; j++) - { - for (py::ssize_t i = 0; i < 3; i++) - { - blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); - } - } - - auto blaze_output = elastica::batch_norm(blaze_vector_batch); - for (py::ssize_t j = 0; j < num_elems; j++) - { - output_arr_unchecked(j) = blaze_output[j]; - } - - return output_arr; + return elastica::batch_norm(vector_batch); } PYBIND11_MODULE(_linalg, m) @@ -328,16 +77,6 @@ PYBIND11_MODULE(_linalg, m) batch_norm )pbdoc"; - m.def("batch_matmul_naive", &batch_matmul_naive, R"pbdoc( - This is batch matrix matrix multiplication function. Only batch - of 3x3 matrices can be multiplied. - )pbdoc"); - - m.def("batch_matmul_blaze", &batch_matmul_blaze, R"pbdoc( - This is batch matrix matrix multiplication function. Only batch - of 3x3 matrices can be multiplied. - )pbdoc"); - m.def("difference_kernel", &difference_kernel, R"pbdoc( Vector Difference )pbdoc"); diff --git a/backend/src/_linalg_numpy.cpp b/backend/src/_linalg_numpy.cpp new file mode 100644 index 000000000..061d3df07 --- /dev/null +++ b/backend/src/_linalg_numpy.cpp @@ -0,0 +1,376 @@ +#include +#include +#include + +#include "Utilities/Math/BlazeDetail/BlazeLinearAlgebra.hpp" + +namespace py = pybind11; + +using blaze::DynamicMatrix; +using blaze::DynamicTensor; +using blaze::StaticMatrix; + +/* Simple rewrite of elastica._linalg._batch_matmul */ +py::array_t batch_matmul_naive( + py::array_t first_matrix_collection, + py::array_t second_matrix_collection) +{ + auto s1 = first_matrix_collection.shape(0); + auto s2 = first_matrix_collection.shape(1); + auto max_k = first_matrix_collection.shape(2); + auto output_matrix = py::array_t{{s1, s2, max_k}}; + + auto a = first_matrix_collection.unchecked<3>(); + auto b = second_matrix_collection.unchecked<3>(); + auto c = output_matrix.mutable_unchecked<3>(); + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t m = 0; m < 3; m++) + { + for (py::ssize_t k = 0; k < max_k; k++) + { + c(i, m, k) += a(i, j, k) * b(j, m, k); + } + } + } + } + return output_matrix; +} + +/* blaze implementation of elastica._linalg._batch_matmul */ +py::array_t batch_matmul_blaze( + py::array_t first_matrix_collection, + py::array_t second_matrix_collection) +{ + auto s1 = first_matrix_collection.shape(0); + auto s2 = first_matrix_collection.shape(1); + auto max_k = first_matrix_collection.shape(2); + auto output_matrix = py::array_t{{s1, s2, max_k}}; + + auto a = first_matrix_collection.unchecked<3>(); + auto b = second_matrix_collection.unchecked<3>(); + auto c = output_matrix.mutable_unchecked<3>(); + + StaticMatrix blaze_a; + StaticMatrix blaze_b; + StaticMatrix blaze_c; + for (py::ssize_t k = 0; k < max_k; k++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + blaze_a(i, j) = a(i, j, k); + blaze_b(i, j) = b(i, j, k); + } + } + blaze_c = blaze_a * blaze_b; + for (py::ssize_t i = 0; i < 3; i++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + c(i, j, k) = blaze_c(i, j); + } + } + } + + return output_matrix; +} + +py::array_t difference_kernel(py::array_t vector_batch) +{ + const auto v0 = vector_batch.shape(0); + const auto num_elems = vector_batch.shape(1); + + assert(v0 == 3UL); + auto output_arr = py::array_t{{v0, num_elems - 1}}; + + auto vector_batch_unchecked = vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicMatrix blaze_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); + } + } + auto blaze_output = elastica::difference_kernel(blaze_vector_batch); + for (py::ssize_t j = 0; j < num_elems - 1; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_matvec( + py::array_t matrix_collection, py::array_t vector_collection) +{ + const auto v0 = vector_collection.shape(0); + const auto num_elems = vector_collection.shape(1); + assert(v0 == 3UL); + + assert(matrix_collection.shape(0) == 3UL); + assert(matrix_collection.shape(1) == 3UL); + assert(matrix_collection.shape(2) == num_elems); + + auto output_arr = py::array_t{{v0, num_elems}}; + + auto matrix_collection_unchecked = matrix_collection.unchecked<3>(); + auto vector_collection_unchecked = vector_collection.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicTensor blaze_matrix_collection(v0, v0, num_elems); + DynamicMatrix blaze_vector_collection(v0, num_elems); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_matrix_collection(i, j, k) = matrix_collection_unchecked(i, j, k); + } + blaze_vector_collection(j, k) = vector_collection_unchecked(j, k); + } + } + + DynamicMatrix blaze_output(v0, num_elems); + elastica::batch_matvec(blaze_output, blaze_matrix_collection, blaze_vector_collection); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_matmul( + py::array_t first_matrix_batch, py::array_t second_matrix_batch) +{ + const auto m0 = first_matrix_batch.shape(0); + const auto m1 = first_matrix_batch.shape(1); + const auto num_elems = first_matrix_batch.shape(2); + assert(m0 == 3UL); + assert(m1 == 3UL); + + assert(second_matrix_batch.shape(0) == 3UL); + assert(second_matrix_batch.shape(1) == 3UL); + assert(second_matrix_batch.shape(2) == num_elems); + + auto output_arr = py::array_t{{m0, m1, num_elems}}; + + auto first_matrix_batch_unchecked = first_matrix_batch.unchecked<3>(); + auto second_matrix_batch_unchecked = second_matrix_batch.unchecked<3>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<3>(); + + DynamicTensor blaze_first_matrix_batch(m0, m1, num_elems); + DynamicTensor blaze_second_matrix_batch(m0, m1, num_elems); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_matrix_batch(i, j, k) = first_matrix_batch_unchecked(i, j, k); + blaze_second_matrix_batch(i, j, k) = second_matrix_batch_unchecked(i, j, k); + } + } + } + DynamicTensor blaze_output(m0, m1, num_elems); + elastica::batch_matmul(blaze_output, blaze_first_matrix_batch, blaze_second_matrix_batch); + for (py::ssize_t k = 0; k < num_elems; k++) + { + for (py::ssize_t j = 0; j < 3; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j, k) = blaze_output(i, j, k); + } + } + } + return output_arr; +} + +py::array_t batch_cross( + py::array_t first_vector_batch, py::array_t second_vector_batch) +{ + const auto v0 = first_vector_batch.shape(0); + const auto num_elems = first_vector_batch.shape(1); + assert(v0 == 3UL); + + assert(second_vector_batch.shape(0) == 3UL); + assert(second_vector_batch.shape(1) == num_elems); + + auto output_arr = py::array_t{{v0, num_elems}}; + + auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); + auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<2>(); + + DynamicMatrix blaze_first_vector_batch(v0, num_elems); + DynamicMatrix blaze_second_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); + blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); + } + } + + DynamicMatrix blaze_output(v0, num_elems); + elastica::batch_cross(blaze_output, blaze_first_vector_batch, blaze_second_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + output_arr_unchecked(i, j) = blaze_output(i, j); + } + } + + return output_arr; +} + +py::array_t batch_dot( + py::array_t first_vector_batch, py::array_t second_vector_batch) +{ + const auto v0 = first_vector_batch.shape(0); + const auto num_elems = first_vector_batch.shape(1); + assert(v0 == 3UL); + + assert(second_vector_batch.shape(0) == 3UL); + assert(second_vector_batch.shape(1) == num_elems); + + auto output_arr = py::array_t{num_elems}; + + auto first_vector_batch_unchecked = first_vector_batch.unchecked<2>(); + auto second_vector_batch_unchecked = second_vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); + + DynamicMatrix blaze_first_vector_batch(v0, num_elems); + DynamicMatrix blaze_second_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_first_vector_batch(i, j) = first_vector_batch_unchecked(i, j); + blaze_second_vector_batch(i, j) = second_vector_batch_unchecked(i, j); + } + } + + auto blaze_output = elastica::batch_dot(blaze_first_vector_batch, blaze_second_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + output_arr_unchecked(j) = blaze_output[j]; + } + + return output_arr; +} + +py::array_t batch_norm(py::array_t vector_batch) +{ + const auto v0 = vector_batch.shape(0); + const auto num_elems = vector_batch.shape(1); + + assert(v0 == 3UL); + + auto output_arr = py::array_t{num_elems}; + + auto vector_batch_unchecked = vector_batch.unchecked<2>(); + auto output_arr_unchecked = output_arr.mutable_unchecked<1>(); + + DynamicMatrix blaze_vector_batch(v0, num_elems); + for (py::ssize_t j = 0; j < num_elems; j++) + { + for (py::ssize_t i = 0; i < 3; i++) + { + blaze_vector_batch(i, j) = vector_batch_unchecked(i, j); + } + } + + auto blaze_output = elastica::batch_norm(blaze_vector_batch); + for (py::ssize_t j = 0; j < num_elems; j++) + { + output_arr_unchecked(j) = blaze_output[j]; + } + + return output_arr; +} + +PYBIND11_MODULE(_linalg_numpy, m) +{ + m.doc() = R"pbdoc( + elasticapp _linalg_numpy + --------------- + + .. currentmodule:: _linalg_numpy + + .. autosummary:: + :toctree: _generate + + batch_matmul_naive + batch_matmul_blaze + difference_kernel + batch_matvec + batch_matmul + batch_cross + batch_dot + batch_norm + )pbdoc"; + + m.def("batch_matmul_naive", &batch_matmul_naive, R"pbdoc( + This is batch matrix matrix multiplication function. Only batch + of 3x3 matrices can be multiplied. + )pbdoc"); + + m.def("batch_matmul_blaze", &batch_matmul_blaze, R"pbdoc( + This is batch matrix matrix multiplication function. Only batch + of 3x3 matrices can be multiplied. + )pbdoc"); + + m.def("difference_kernel", &difference_kernel, R"pbdoc( + Vector Difference + )pbdoc"); + + m.def("batch_matvec", &batch_matvec, R"pbdoc( + Computes a batchwise matrix-vector product given in indical notation: + matvec_batch{ik} = matrix_batch{ijk} * vector_batch{jk} + )pbdoc"); + + m.def("batch_matmul", &batch_matmul, R"pbdoc( + Computes a batchwise matrix-matrix product given in indical notation: + matmul_batch{ilk} = first_matrix_batch{ijk} * second_matrix_batch{jlk} + )pbdoc"); + + m.def("batch_cross", &batch_cross, R"pbdoc( + Computes a batchwise vector-vector cross product given in indical notation: + cross_batch{il} = LCT{ijk} * first_vector_batch{jl} * second_vector_batch{kl} + where LCT is the Levi-Civita Tensor + )pbdoc"); + + m.def("batch_dot", &batch_dot, R"pbdoc( + Computes a batchwise vector-vector dot product given in indical notation: + dot_batch{j} = first_vector_batch{ij} * second_vector_batch{ij} + )pbdoc"); + + m.def("batch_norm", &batch_norm, R"pbdoc( + Computes a batchwise vector L2 norm given in indical notation: + norm_batch{j} = (vector_batch{ij} * vector_batch{ij})^0.5 + )pbdoc"); + +#ifdef VERSION_INFO + m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); +#else + m.attr("__version__") = "dev"; +#endif +} diff --git a/backend/src/meson.build b/backend/src/meson.build index 623d13960..42b030ee4 100644 --- a/backend/src/meson.build +++ b/backend/src/meson.build @@ -16,6 +16,15 @@ _linalg = py.extension_module( subdir: package, ) +_linalg_numpy = py.extension_module( + '_linalg_numpy', + sources: ['_linalg_numpy.cpp'], + dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], + link_language: 'cpp', + install: true, + subdir: package, +) + _PyArrays = py.extension_module( '_PyArrays', sources: [ diff --git a/backend/tests/test_linalg.py b/backend/tests/test_linalg.py index 1fba0332b..da92cb370 100644 --- a/backend/tests/test_linalg.py +++ b/backend/tests/test_linalg.py @@ -6,6 +6,7 @@ import numpy as np import pytest from numpy.testing import assert_allclose +from elasticapp._PyArrays import Matrix, Tensor from elasticapp._linalg import ( difference_kernel, batch_matvec, @@ -19,7 +20,9 @@ @pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) def test_difference_kernel(blocksize: int): input_vector_collection = np.random.randn(3, blocksize) - output_vector_collection = difference_kernel(input_vector_collection) + output_vector_collection = np.asarray( + difference_kernel(Matrix(input_vector_collection)) + ) correct_vector_collection = ( input_vector_collection[:, 1:] - input_vector_collection[:, :-1] @@ -33,8 +36,8 @@ def test_batch_matvec(blocksize: int): input_matrix_collection = np.random.randn(3, 3, blocksize) input_vector_collection = np.random.randn(3, blocksize) - test_vector_collection = batch_matvec( - input_matrix_collection, input_vector_collection + test_vector_collection = np.asarray( + batch_matvec(Tensor(input_matrix_collection), Matrix(input_vector_collection)) ) correct_vector_collection = [ @@ -51,8 +54,11 @@ def test_batch_matmul(blocksize: int): input_first_matrix_collection = np.random.randn(3, 3, blocksize) input_second_matrix_collection = np.random.randn(3, 3, blocksize) - test_matrix_collection = batch_matmul( - input_first_matrix_collection, input_second_matrix_collection + test_matrix_collection = np.asarray( + batch_matmul( + Tensor(input_first_matrix_collection), + Tensor(input_second_matrix_collection), + ) ) correct_matrix_collection = np.empty((3, 3, blocksize)) @@ -72,8 +78,11 @@ def test_batch_cross(dim, blocksize: int): input_first_vector_collection = np.random.randn(dim, blocksize) input_second_vector_collection = np.random.randn(dim, blocksize) - test_vector_collection = batch_cross( - input_first_vector_collection, input_second_vector_collection + test_vector_collection = np.asarray( + batch_cross( + Matrix(input_first_vector_collection), + Matrix(input_second_vector_collection), + ) ) correct_vector_collection = np.cross( input_first_vector_collection, input_second_vector_collection, axisa=0, axisb=0 @@ -87,8 +96,11 @@ def test_batch_dot(blocksize: int): input_first_vector_collection = np.random.randn(3, blocksize) input_second_vector_collection = np.random.randn(3, blocksize) - test_vector_collection = batch_dot( - input_first_vector_collection, input_second_vector_collection + test_vector_collection = np.asarray( + batch_dot( + Matrix(input_first_vector_collection), + Matrix(input_second_vector_collection), + ) ) correct_vector_collection = np.einsum( @@ -102,7 +114,9 @@ def test_batch_dot(blocksize: int): def test_batch_norm(blocksize: int): input_first_vector_collection = np.random.randn(3, blocksize) - test_vector_collection = batch_norm(input_first_vector_collection) + test_vector_collection = np.asarray( + batch_norm(Matrix(input_first_vector_collection)) + ) correct_vector_collection = np.sqrt( np.einsum( diff --git a/backend/tests/test_linalg_numpy.py b/backend/tests/test_linalg_numpy.py new file mode 100644 index 000000000..fe4dc0abc --- /dev/null +++ b/backend/tests/test_linalg_numpy.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 + +# This file is based on pyelastica tests/test_math/test_linalg.py + +# System imports +import numpy as np +import pytest +from numpy.testing import assert_allclose +from elasticapp._linalg_numpy import ( + difference_kernel, + batch_matvec, + batch_matmul, + batch_cross, + batch_dot, + batch_norm, +) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_difference_kernel(blocksize: int): + input_vector_collection = np.random.randn(3, blocksize) + output_vector_collection = difference_kernel(input_vector_collection) + + correct_vector_collection = ( + input_vector_collection[:, 1:] - input_vector_collection[:, :-1] + ) + + assert_allclose(output_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_matvec(blocksize: int): + input_matrix_collection = np.random.randn(3, 3, blocksize) + input_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_matvec( + input_matrix_collection, input_vector_collection + ) + + correct_vector_collection = [ + np.dot(input_matrix_collection[..., i], input_vector_collection[..., i]) + for i in range(blocksize) + ] + correct_vector_collection = np.array(correct_vector_collection).T + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_matmul(blocksize: int): + input_first_matrix_collection = np.random.randn(3, 3, blocksize) + input_second_matrix_collection = np.random.randn(3, 3, blocksize) + + test_matrix_collection = batch_matmul( + input_first_matrix_collection, input_second_matrix_collection + ) + + correct_matrix_collection = np.empty((3, 3, blocksize)) + for i in range(blocksize): + correct_matrix_collection[..., i] = np.dot( + input_first_matrix_collection[..., i], + input_second_matrix_collection[..., i], + ) + + assert_allclose(test_matrix_collection, correct_matrix_collection) + + +# TODO : Generalize to two dimensions +@pytest.mark.parametrize("dim", [3]) +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_cross(dim, blocksize: int): + input_first_vector_collection = np.random.randn(dim, blocksize) + input_second_vector_collection = np.random.randn(dim, blocksize) + + test_vector_collection = batch_cross( + input_first_vector_collection, input_second_vector_collection + ) + correct_vector_collection = np.cross( + input_first_vector_collection, input_second_vector_collection, axisa=0, axisb=0 + ).T + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_dot(blocksize: int): + input_first_vector_collection = np.random.randn(3, blocksize) + input_second_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_dot( + input_first_vector_collection, input_second_vector_collection + ) + + correct_vector_collection = np.einsum( + "ij,ij->j", input_first_vector_collection, input_second_vector_collection + ) + + assert_allclose(test_vector_collection, correct_vector_collection) + + +@pytest.mark.parametrize("blocksize", [1, 2, 8, 32]) +def test_batch_norm(blocksize: int): + input_first_vector_collection = np.random.randn(3, blocksize) + + test_vector_collection = batch_norm(input_first_vector_collection) + + correct_vector_collection = np.sqrt( + np.einsum( + "ij,ij->j", input_first_vector_collection, input_first_vector_collection + ) + ) + + assert_allclose(test_vector_collection, correct_vector_collection) From f9d4d001fd6a715bbe371803b629a1d88428a176 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 9 Jul 2024 22:29:12 -0500 Subject: [PATCH 13/62] Update README.md --- backend/README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/backend/README.md b/backend/README.md index 60b6c39ec..33bbee46f 100644 --- a/backend/README.md +++ b/backend/README.md @@ -10,6 +10,8 @@ In the root of the PyElastica source tree, run the following command pip install ./backend ``` +> Make sure you install the package from _PyElastica source tree_. + This command will take care of installation of all build-time dependencies, compilation of C++ source files and install the a python package called `elasticapp`. ## Testing From 7d45e74747344d91009b3d7bb525296d8ea6f3b3 Mon Sep 17 00:00:00 2001 From: Ankith Date: Wed, 10 Jul 2024 17:45:20 +0530 Subject: [PATCH 14/62] Add batchcross benchmarking --- backend/benchmarking/batchcross.py | 50 ++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 backend/benchmarking/batchcross.py diff --git a/backend/benchmarking/batchcross.py b/backend/benchmarking/batchcross.py new file mode 100644 index 000000000..eb9cc78b0 --- /dev/null +++ b/backend/benchmarking/batchcross.py @@ -0,0 +1,50 @@ +from elasticapp._linalg_numpy import batch_cross +from elasticapp._linalg import batch_cross as batch_cross_final +from elasticapp._PyArrays import Matrix +from elastica._linalg import _batch_cross +import numpy +import time + +# warm up jit for fair comparison +random_1 = numpy.random.random((3, 1)) +random_2 = numpy.random.random((3, 1)) +out1 = _batch_cross(random_1, random_2) + + +def benchmark_batchsize(funcs: list, batches: list[int], num_iterations: int = 1000): + ret: dict = {} + for batch_size in batches: + random_a = numpy.random.random((3, batch_size)) + random_b = numpy.random.random((3, batch_size)) + + ret[batch_size] = {} + for func_name, func, func_wrap in funcs: + random_a_w = func_wrap(random_a) if func_wrap else random_a + random_b_w = func_wrap(random_b) if func_wrap else random_b + + start = time.perf_counter() + for _ in range(num_iterations): + func(random_a_w, random_b_w) + + ret[batch_size][func_name] = (time.perf_counter() - start) / num_iterations + + return ret + + +results = benchmark_batchsize( + [ + ("pyelastica", _batch_cross, None), + ("elasticapp_blaze_copy", batch_cross, None), + ("elasticapp_blaze_final", batch_cross_final, Matrix), + ], + [2**i for i in range(14)], +) +for size, data in results.items(): + pyelastica = data["pyelastica"] + elasticapp_blaze_copy = data["elasticapp_blaze_copy"] + elasticapp_blaze_final = data["elasticapp_blaze_final"] + print(f"{size = }") + print(f"{pyelastica = }") + print(f"{elasticapp_blaze_copy = }, ratio: {elasticapp_blaze_copy / pyelastica}") + print(f"{elasticapp_blaze_final = }, ratio: {elasticapp_blaze_final / pyelastica}") + print() From f1eaef8a8da7aa1299b271147017acfecd2b7f0e Mon Sep 17 00:00:00 2001 From: Ankith Date: Wed, 10 Jul 2024 17:48:50 +0530 Subject: [PATCH 15/62] Setup optimization blaze defines globally --- backend/meson.build | 37 +++++++++++++++++++++++++++++++++++++ backend/src/meson.build | 6 +++--- 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/backend/meson.build b/backend/meson.build index 4c5cce8f2..dc97db5f4 100644 --- a/backend/meson.build +++ b/backend/meson.build @@ -12,6 +12,40 @@ cc = meson.get_compiler('cpp') py = import('python').find_installation(pure: false) py_dep = py.dependency() +# set up global defines to optimize blaze usages +add_global_arguments( + # - Enable external BLAS kernels + '-DBLAZE_BLAS_MODE=0', + # - Set default matrix storage order to column-major, since many of our + # functions are implemented for column-major layout. This default reduces + # conversions. + '-DBLAZE_DEFAULT_STORAGE_ORDER=blaze::rowMajor', + # - Disable SMP parallelization. This disables SMP parallelization for all + # possible backends (OpenMP, C++11 threads, Boost, HPX): + # https://bitbucket.org/blaze-lib/blaze/wiki/Serial%20Execution#!option-3-deactivation-of-parallel-execution + '-DBLAZE_USE_SHARED_MEMORY_PARALLELIZATION=0', + # - Disable MPI parallelization + '-DBLAZE_MPI_PARALLEL_MODE=0', + # - Using the default cache size, which may have been configured automatically + # by the Blaze CMake configuration for the machine we are running on. We + # could override it here explicitly to tune performance. + # BLAZE_CACHE_SIZE + '-DBLAZE_USE_PADDING=1', + # - Always enable non-temporal stores for cache optimization of large data + # structures: https://bitbucket.org/blaze-lib/blaze/wiki/Configuration%20Files#!streaming-non-temporal-stores + '-DBLAZE_USE_STREAMING=1', + # - Initializing default-constructed structures for fundamental types + '-DBLAZE_USE_DEFAULT_INITIALIZATON=1', + # Use Sleef for vectorization of more math functions + '-DBLAZE_USE_SLEEF=1', + # Set inlining settings + '-DBLAZE_USE_STRONG_INLINE=1', + '-DBLAZE_USE_ALWAYS_INLINE=1', + # Set vectorization (leave to 1, else there is no use for blaze) + '-DBLAZE_USE_VECTORIZATION=1', + language: 'cpp' +) + # find dependencies and create dep objects pybind11_dep = declare_dependency( include_directories: run_command( @@ -53,6 +87,9 @@ if not sleef_dep.found() ) endif +# blaze and deps commonly used together +blaze_deps = [blaze_dep, blaze_tensor_dep, sleef_dep] + brigand_dep = dependency('brigand', fallback : 'brigand') cxxopts_dep = dependency('cxxopts', fallback : 'cxxopts') diff --git a/backend/src/meson.build b/backend/src/meson.build index 42b030ee4..fcbda76b8 100644 --- a/backend/src/meson.build +++ b/backend/src/meson.build @@ -10,7 +10,7 @@ example_1 = py.extension_module( _linalg = py.extension_module( '_linalg', sources: ['_linalg.cpp'], - dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], + dependencies: [py_dep, pybind11_dep] + blaze_deps, link_language: 'cpp', install: true, subdir: package, @@ -19,7 +19,7 @@ _linalg = py.extension_module( _linalg_numpy = py.extension_module( '_linalg_numpy', sources: ['_linalg_numpy.cpp'], - dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], + dependencies: [py_dep, pybind11_dep] + blaze_deps, link_language: 'cpp', install: true, subdir: package, @@ -33,7 +33,7 @@ _PyArrays = py.extension_module( 'Utilities/Math/Python/BlazeMatrix.cpp', 'Utilities/Math/Python/BlazeTensor.cpp', ], - dependencies: [py_dep, pybind11_dep, blaze_dep, blaze_tensor_dep], + dependencies: [py_dep, pybind11_dep] + blaze_deps, link_language: 'cpp', install: true, subdir: package, From a61a4450cde064774bdf169ccfa1393fbd898597 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 10 Jul 2024 21:49:04 -0500 Subject: [PATCH 16/62] revert: __str__ for blaze objects. Other print-related methods are added. --- .../src/Utilities/Math/Python/BlazeMatrix.cpp | 14 ++-- .../src/Utilities/Math/Python/BlazeTensor.cpp | 12 +-- .../src/Utilities/Math/Python/BlazeVector.cpp | 20 +++-- backend/src/Utilities/PrintHelpers.hpp | 72 ++++++++++++++++++ backend/src/Utilities/Requires.hpp | 70 +++++++++++++++++ .../src/Utilities/StlStreamDeclarations.hpp | 76 +++++++++++++++++++ 6 files changed, 246 insertions(+), 18 deletions(-) create mode 100644 backend/src/Utilities/PrintHelpers.hpp create mode 100644 backend/src/Utilities/Requires.hpp create mode 100644 backend/src/Utilities/StlStreamDeclarations.hpp diff --git a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp index 56f55a110..4ebe13deb 100644 --- a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp +++ b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp @@ -6,7 +6,7 @@ // #include "PythonBindings/BoundChecks.hpp" // #include "Utilities/DefineTypes.h" -// #include "Utilities/MakeString.hpp" +#include "Utilities/MakeString.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" // @@ -81,7 +81,7 @@ namespace py_bindings { "__getitem__", +[](const type& self, const std::tuple& x) { - // matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); + matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); return self(std::get<0>(x), std::get<1>(x)); }) .def( @@ -90,16 +90,16 @@ namespace py_bindings { return array_slice(t, std::move(slice)); }) // Need __str__ for converting to string/printing - // .def( - // "__str__", - // +[](const type& self) { return std::string(MakeString{} << self); }) + .def( + "__str__", + +[](const type& self) { return std::string(MakeString{} << self); }) .def( "__setitem__", +[](type& self, const std::tuple& x, const Real val) { - // matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); + matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); self(std::get<0>(x), std::get<1>(x)) = val; - }); + }); } //**************************************************************************** diff --git a/backend/src/Utilities/Math/Python/BlazeTensor.cpp b/backend/src/Utilities/Math/Python/BlazeTensor.cpp index 156f7f891..9ed4d01c4 100644 --- a/backend/src/Utilities/Math/Python/BlazeTensor.cpp +++ b/backend/src/Utilities/Math/Python/BlazeTensor.cpp @@ -6,7 +6,7 @@ // #include "PythonBindings/BoundChecks.hpp" // #include "Utilities/DefineTypes.h" -// #include "Utilities/MakeString.hpp" +#include "Utilities/MakeString.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" // @@ -99,12 +99,12 @@ namespace py_bindings { // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), // std::get<2>(x)); self(std::get<0>(x), std::get<1>(x), std::get<2>(x)) = val; - }); + }) // Need __str__ for converting to string/printing - // .def( - // "__str__", +[](const type& self) { - // return std::string(MakeString{} << self); - // }); + .def( + "__str__", +[](const type& self) { + return std::string(MakeString{} << self); + }); } //**************************************************************************** diff --git a/backend/src/Utilities/Math/Python/BlazeVector.cpp b/backend/src/Utilities/Math/Python/BlazeVector.cpp index 86239102f..eed3afa66 100644 --- a/backend/src/Utilities/Math/Python/BlazeVector.cpp +++ b/backend/src/Utilities/Math/Python/BlazeVector.cpp @@ -5,6 +5,8 @@ // #include "PythonBindings/BoundChecks.hpp" // #include "Utilities/DefineTypes.h" +// #include "Utilities/PrettyType.hpp" +#include "Utilities/PrintHelpers.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" // @@ -96,7 +98,7 @@ namespace py_bindings { .def( "__getitem__", +[](const type& t, const size_t i) { - // bounds_check(t, i); + bounds_check(t, i); return t[i]; }) .def( @@ -106,17 +108,25 @@ namespace py_bindings { }) .def( "__setitem__", +[](type& t, const size_t i, const Real v) { - // bounds_check(t, i); + bounds_check(t, i); t[i] = v; }); + static const auto printer = [](type const &t) { + // Blaze's default printing adds extra lines and spaces which is + // not what we want + std::ostringstream os; + sequence_print_helper(os, t.begin(), t.end()); + return os.str(); + }; + // Need __str__ for converting to string py_vector - // .def("__str__") + .def("__str__", +printer) // repr allows you to output the object in an interactive python // terminal using obj to get the "string REPResenting the object". - // .def("__repr__") - // .def(py::self += py::self) + .def("__repr__", +printer) + .def(py::self += py::self) // Need to do math explicitly converting to DataVector because we don't // want to represent all the possible expression template types .def( diff --git a/backend/src/Utilities/PrintHelpers.hpp b/backend/src/Utilities/PrintHelpers.hpp new file mode 100644 index 000000000..b2134f5a2 --- /dev/null +++ b/backend/src/Utilities/PrintHelpers.hpp @@ -0,0 +1,72 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include + +#include "Utilities/StlStreamDeclarations.hpp" + +template +void sequence_print_helper(std::ostream& out, ForwardIt&& begin, + ForwardIt&& end, Func f) noexcept { + out << "("; + if (begin != end) { + while (true) { + f(out, begin++); + if (begin == end) { + break; + } + out << ","; + } + } + out << ")"; +} + +/*! + * \ingroup UtilitiesGroup + * \brief Prints all the items as a comma separated list surrounded by parens. + */ +template +void sequence_print_helper(std::ostream& out, ForwardIt&& begin, + ForwardIt&& end) noexcept { + sequence_print_helper(out, std::forward(begin), + std::forward(end), + [](std::ostream& os, const ForwardIt& it) noexcept { + using ::operator<<; + os << *it; + }); +} + +//@{ +/*! + * \ingroup UtilitiesGroup + * Like sequence_print_helper, but sorts the string representations. + */ +template +void unordered_print_helper(std::ostream& out, ForwardIt&& begin, + ForwardIt&& end, Func f) noexcept { + std::vector entries; + while (begin != end) { + std::ostringstream ss; + f(ss, begin++); + entries.push_back(ss.str()); + } + std::sort(entries.begin(), entries.end()); + sequence_print_helper(out, entries.begin(), entries.end()); +} + +template +void unordered_print_helper(std::ostream& out, ForwardIt&& begin, + ForwardIt&& end) noexcept { + unordered_print_helper(out, std::forward(begin), + std::forward(end), + [](std::ostream& os, const ForwardIt& it) noexcept { + using ::operator<<; + os << *it; + }); +} diff --git a/backend/src/Utilities/Requires.hpp b/backend/src/Utilities/Requires.hpp new file mode 100644 index 000000000..4a74ce9b1 --- /dev/null +++ b/backend/src/Utilities/Requires.hpp @@ -0,0 +1,70 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +// Source : +// https://raw.githubusercontent.com/sxs-collaboration/spectre/develop/src/Utilities/Requires.hpp + +#pragma once + +/// \file +/// Defines the type alias Requires + +#include + +namespace Requires_detail { + template + struct requires_impl { + using template_error_type_failed_to_meet_requirements_on_template_parameters = + std::nullptr_t; + }; + + template <> + struct requires_impl {}; +} // namespace Requires_detail + +/*! + * \ingroup UtilitiesGroup + * \brief Express requirements on the template parameters of a function or + * class, replaces `std::enable_if_t` + * + * Replacement for `std::enable_if_t` and Concepts for expressing requirements + * on template parameters. This does not require merging of the Concepts + * TS (whose merit is debatable) and provides an "error message" if substitution + * of a template parameter failed. Specifically, the compiler error will contain + * "template_error_type_failed_to_meet_requirements_on_template_parameters", + * aiding the user of a function or class in tracking down the list of + * requirements on the deduced type. + * + * For example, if a function `foo` is defined as: + * \snippet Utilities/Test_Requires.cpp foo_definition + * then calling the function with a list, `foo(std::list{});` results in + * the following compilation error from clang: + * \code + * ./tests/Unit/Utilities/Test_Requires.cpp:29:3: error: no matching function + * for call to 'foo' + * foo(std::list{}); + * ^~~ + * ./tests/Unit/Utilities/Test_Requires.cpp:15:13: note: candidate + * template ignored: substitution failure [with T = std::__1::list >]: no type named + * 'template_error_type_failed_to_meet_requirements_on_template_parameters' + * in 'Requires_detail::requires_impl' + * std::string foo(const T&) { + * ^ + * 1 error generated. + * \endcode + * + * Here is an example of how write function overloads using `Requires` or to + * express constraints on the template parameters: + * \snippet Utilities/Test_Requires.cpp function_definitions + * + * \note + * Using `Requires` is safer than using `std::enable_if_t` because the + * nested type alias is of type `std::nullptr_t` and so usage is always: + * \code + * template = nullptr> + * \endcode + */ +template +using Requires = typename Requires_detail::requires_impl< + B>::template_error_type_failed_to_meet_requirements_on_template_parameters; diff --git a/backend/src/Utilities/StlStreamDeclarations.hpp b/backend/src/Utilities/StlStreamDeclarations.hpp new file mode 100644 index 000000000..3df51f5bb --- /dev/null +++ b/backend/src/Utilities/StlStreamDeclarations.hpp @@ -0,0 +1,76 @@ +// Reused from SpECTRE : https://spectre-code.org/ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +// Suppress warnings from this header since they can occur from not including +// StdHelpers later. This header is used only in TypeTraits.hpp so that +// is_streamable works correctly for STL types that have stream operators +// defined. +#ifdef __GNUC__ +#pragma GCC system_header +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Utilities/Requires.hpp" + +namespace tt { + template + struct is_streamable; +} // namespace tt + +template +std::ostream& operator<<(std::ostream& os, const std::list& v) noexcept; + +template +std::ostream& operator<<(std::ostream& os, const std::vector& v) noexcept; + +template +std::ostream& operator<<(std::ostream& os, const std::deque& v) noexcept; + +template +std::ostream& operator<<(std::ostream& os, const std::array& a) noexcept; + +template +std::ostream& operator<<(std::ostream& os, + const std::tuple& t) noexcept; + +template +inline std::ostream& operator<<(std::ostream& os, + const std::unordered_map& m) noexcept; + +template +inline std::ostream& operator<<(std::ostream& os, + const std::map& m) noexcept; + +template +std::ostream& operator<<(std::ostream& os, + const std::unordered_set& v) noexcept; + +template +inline std::ostream& operator<<(std::ostream& os, + const std::set& v) noexcept; + +template ::value> = nullptr> +std::ostream& operator<<(std::ostream& os, + const std::unique_ptr& t) noexcept; + +template ::value> = nullptr> +std::ostream& operator<<(std::ostream& os, + const std::shared_ptr& t) noexcept; + +template +std::ostream& operator<<(std::ostream& os, const std::pair& t) noexcept; From 2e34cdd2030fbb5a28e6bfc8d9d75967d25f66f1 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 10 Jul 2024 21:52:33 -0500 Subject: [PATCH 17/62] feat: Add bound check during blaze setter --- backend/README.md | 9 +++ .../src/Utilities/Math/Python/BlazeMatrix.cpp | 6 +- .../src/Utilities/Math/Python/BlazeTensor.cpp | 11 ++- .../src/Utilities/Math/Python/BlazeVector.cpp | 3 +- .../src/Utilities/Math/Python/BoundChecks.hpp | 67 +++++++++++++++++++ 5 files changed, 84 insertions(+), 12 deletions(-) create mode 100644 backend/src/Utilities/Math/Python/BoundChecks.hpp diff --git a/backend/README.md b/backend/README.md index 33bbee46f..506dfeb2f 100644 --- a/backend/README.md +++ b/backend/README.md @@ -33,3 +33,12 @@ For benchmarking various `matmul` implementations, run ``` python3 backend/benchmarking/matmul.py ``` + +## Contributed By + +- Tejaswin Parthasarathy (Teja) +- [Seung Hyun Kim](https://github.com/skim0119) +- Ankith Pai +- [Yashraj Bhosale](https://github.com/bhosale2) +- Arman Tekinalp +- Songyuan Cui diff --git a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp index 4ebe13deb..2834471ab 100644 --- a/backend/src/Utilities/Math/Python/BlazeMatrix.cpp +++ b/backend/src/Utilities/Math/Python/BlazeMatrix.cpp @@ -2,13 +2,11 @@ // Includes //****************************************************************************** -// -// #include "PythonBindings/BoundChecks.hpp" -// #include "Utilities/DefineTypes.h" #include "Utilities/MakeString.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" +#include "Utilities/Math/Python/BoundChecks.hpp" // #include #include @@ -99,7 +97,7 @@ namespace py_bindings { const Real val) { matrix_bounds_check(self, std::get<0>(x), std::get<1>(x)); self(std::get<0>(x), std::get<1>(x)) = val; - }); + }); } //**************************************************************************** diff --git a/backend/src/Utilities/Math/Python/BlazeTensor.cpp b/backend/src/Utilities/Math/Python/BlazeTensor.cpp index 9ed4d01c4..779952b77 100644 --- a/backend/src/Utilities/Math/Python/BlazeTensor.cpp +++ b/backend/src/Utilities/Math/Python/BlazeTensor.cpp @@ -3,12 +3,11 @@ //****************************************************************************** -// #include "PythonBindings/BoundChecks.hpp" -// #include "Utilities/DefineTypes.h" #include "Utilities/MakeString.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" +#include "Utilities/Math/Python/BoundChecks.hpp" // #include #include @@ -82,8 +81,8 @@ namespace py_bindings { "__getitem__", +[](const type& self, const std::tuple& x) { - // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), - // std::get<2>(x)); + tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), + std::get<2>(x)); return self(std::get<0>(x), std::get<1>(x), std::get<2>(x)); }) .def( @@ -96,8 +95,8 @@ namespace py_bindings { +[](type& self, const std::tuple& x, const Real val) { - // tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), - // std::get<2>(x)); + tensor_bounds_check(self, std::get<0>(x), std::get<1>(x), + std::get<2>(x)); self(std::get<0>(x), std::get<1>(x), std::get<2>(x)) = val; }) // Need __str__ for converting to string/printing diff --git a/backend/src/Utilities/Math/Python/BlazeVector.cpp b/backend/src/Utilities/Math/Python/BlazeVector.cpp index eed3afa66..1d9494ac1 100644 --- a/backend/src/Utilities/Math/Python/BlazeVector.cpp +++ b/backend/src/Utilities/Math/Python/BlazeVector.cpp @@ -2,13 +2,12 @@ // Includes //****************************************************************************** -// #include "PythonBindings/BoundChecks.hpp" // #include "Utilities/DefineTypes.h" -// #include "Utilities/PrettyType.hpp" #include "Utilities/PrintHelpers.hpp" // #include "Utilities/Math/Python/SliceHelpers.hpp" +#include "Utilities/Math/Python/BoundChecks.hpp" // #include #include diff --git a/backend/src/Utilities/Math/Python/BoundChecks.hpp b/backend/src/Utilities/Math/Python/BoundChecks.hpp new file mode 100644 index 000000000..84f912457 --- /dev/null +++ b/backend/src/Utilities/Math/Python/BoundChecks.hpp @@ -0,0 +1,67 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +#include +#include + +#include "Utilities/PrettyType.hpp" + +namespace py_bindings { + + //**************************************************************************** + /*! \brief Check if a vector-like object access is in bounds. Throws + * std::runtime_error if it is not. + * \ingroup python_bindings + */ + template + void bounds_check(const T& t, const std::size_t i) { + if (i >= t.size()) { + throw std::runtime_error{"Out of bounds access (" + std::to_string(i) + + ") into " + pretty_type::name() + + " of size " + std::to_string(t.size())}; + } + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Check if a matrix-like object access is in bounds. Throws + * std::runtime_error if it is not. + * \ingroup python_bindings + */ + template + void matrix_bounds_check(const T& matrix, const std::size_t row, + const std::size_t column) { + if (row >= matrix.rows() or column >= matrix.columns()) { + throw std::runtime_error{"Out of bounds access (" + std::to_string(row) + + ", " + std::to_string(column) + + ") into Matrix of size (" + + std::to_string(matrix.rows()) + ", " + + std::to_string(matrix.columns()) + ")"}; + } + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Check if a tensor-like object access is in bounds. Throws + * std::runtime_error if it is not. + * \ingroup python_bindings + */ + template + void tensor_bounds_check(const T& tensor, const std::size_t page, + const std::size_t row, const std::size_t column) { + if (page >= tensor.pages() or row >= tensor.rows() or + column >= tensor.columns()) { + throw std::runtime_error{ + "Out of bounds access (" + std::to_string(page) + ", " + + std::to_string(row) + ", " + std::to_string(column) + + ") into Tensor of size (" + std::to_string(tensor.pages()) + ", " + + std::to_string(tensor.rows()) + ", " + + std::to_string(tensor.columns()) + ")"}; + } + } + //**************************************************************************** + +} // namespace py_bindings From 99ff10f83cb7a4025e6d2bff3e79404c752e3d0c Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 10 Jul 2024 22:03:23 -0500 Subject: [PATCH 18/62] Add default simulator frames definition --- backend/src/Simulator/Frames.hpp | 95 ++++++++++++++++ .../src/Simulator/Frames/EulerianFrame.hpp | 90 +++++++++++++++ .../src/Simulator/Frames/LagrangianFrame.hpp | 103 ++++++++++++++++++ .../Simulator/Frames/RotationConvention.hpp | 27 +++++ backend/src/Utilities/Math/Types.hpp | 13 --- backend/src/Utilities/Math/Vec3.hpp | 13 --- backend/src/Utilities/NonCreatable.hpp | 77 +++++++++++++ 7 files changed, 392 insertions(+), 26 deletions(-) create mode 100644 backend/src/Simulator/Frames.hpp create mode 100644 backend/src/Simulator/Frames/EulerianFrame.hpp create mode 100644 backend/src/Simulator/Frames/LagrangianFrame.hpp create mode 100644 backend/src/Simulator/Frames/RotationConvention.hpp create mode 100644 backend/src/Utilities/NonCreatable.hpp diff --git a/backend/src/Simulator/Frames.hpp b/backend/src/Simulator/Frames.hpp new file mode 100644 index 000000000..2e46a78f9 --- /dev/null +++ b/backend/src/Simulator/Frames.hpp @@ -0,0 +1,95 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include "Simulator/Frames/EulerianFrame.hpp" +#include "Simulator/Frames/LagrangianFrame.hpp" +#include "Utilities/Math/Vec3.hpp" + +namespace elastica { + + // ? + using EulerianFrame = detail::EulerianFrame; + using detail::cast_along; + using LagrangianFrame = detail::LagrangianFrame; + // ? + + //**************************************************************************** + /*!\brief Gets a unit vector along an Eulerian direction + * \ingroup simulator + * + * \details + * Gets the unit vector along a direction `Dir` of an Eulerian Frame + * + * \param dir Direction along which to get a unit vector + */ + inline constexpr auto get_unit_vector_along( + EulerianFrame::DirectionType dir) noexcept -> Vec3 { + return { + cast_along(dir) == 0 ? 1.0 : 0.0, + cast_along(dir) == 1 ? 1.0 : 0.0, + cast_along(dir) == 2 ? 1.0 : 0.0, + }; + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Gets a unit vector along an Eulerian direction + * \ingroup simulator + * + * \details + * Gets the unit vector along a direction `Dir` of an Eulerian Frame + * + * \tparam Dir Direction along which to get a unit vector + */ + template + inline constexpr auto get_unit_vector_along() noexcept -> Vec3 { + return get_unit_vector_along(Dir); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Gets a unit vector along a Lagrangian direction + * \ingroup simulator + * + * \details + * Gets the unit vector along a direction `Dir` of an Lagrangian Frame + * + * \example + * \snippet Test_LagrangianFrame.cpp vector_along_eg + * + * \param dir Direction along which to get a unit vector + */ + inline constexpr auto get_unit_vector_along( + LagrangianFrame::DirectionType dir) noexcept -> Vec3 { + return { + cast_along(dir) == 0 ? 1.0 : 0.0, + cast_along(dir) == 1 ? 1.0 : 0.0, + cast_along(dir) == 2 ? 1.0 : 0.0, + }; + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Gets a unit vector along a Lagrangian direction + * \ingroup simulator + * + * \details + * Gets the unit vector along a direction `Dir` of an Lagrangian Frame + * + * \tparam Dir Direction along which to get a unit vector + */ + template + inline constexpr auto get_unit_vector_along() noexcept -> Vec3 { + return get_unit_vector_along(Dir); + } + //**************************************************************************** + + struct Frames { + static_assert(LagrangianFrame::Dimension == EulerianFrame::Dimension, + "Wrong dimension configuration!"); + static constexpr auto Dimension = LagrangianFrame::Dimension; + }; + +} // namespace elastica diff --git a/backend/src/Simulator/Frames/EulerianFrame.hpp b/backend/src/Simulator/Frames/EulerianFrame.hpp new file mode 100644 index 000000000..c48a35011 --- /dev/null +++ b/backend/src/Simulator/Frames/EulerianFrame.hpp @@ -0,0 +1,90 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include + +#include "Utilities/NonCreatable.hpp" + +namespace elastica { + + namespace detail { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*!\brief Trait defining unique IDs for directions in Eulerian frame + * \ingroup simulator + * + * EulerianFrame is a trait defining the type and number of directions in + * 3D. It serves as a strong type for all directions within the \elastica + * library. + */ + struct EulerianFrame : public NonCreatable { + public: + //**Type definitions****************************************************** + //! The type of direction ID + enum class DirectionType : std::uint8_t { x = 0, y, z, Count }; + //************************************************************************ + + //**Static members******************************************************** + //! unique ID for X direction + static constexpr DirectionType X = DirectionType::x; + //! unique ID for Y direction + static constexpr DirectionType Y = DirectionType::y; + //! unique ID for Z direction + static constexpr DirectionType Z = DirectionType::z; + //! Dimension of simulation + static constexpr auto Dimension = + static_cast(DirectionType::Count); + //************************************************************************ + }; + //************************************************************************** + + //************************************************************************** + /*!\brief Casts a Eulerian frame direction into an index value + * \ingroup simulator + * + * \details + * Converts a EulerianFrame::DirectionType into an index value in + * associative arrays. + * + * \example + * \snippet Test_EulerianFrame.cpp cast_along_eg + * + * \param dir A Direction to cast as index to + */ + inline constexpr auto cast_along(EulerianFrame::DirectionType dir) + -> std::size_t { + return static_cast(dir); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Casts a Eulerian frame direction into an index value + * \ingroup simulator + * + * \details + * Converts a EulerianFrame::DirectionType into an index value in + * associative arrays. + * + * \example + * \snippet Test_EulerianFrame.cpp cast_along_template_eg + * + * \tparam Dir A Direction to cast as index to + */ + template + inline constexpr auto cast_along() -> std::size_t { + return cast_along(Dir); + } + //************************************************************************** + + } // namespace detail + +} // namespace elastica diff --git a/backend/src/Simulator/Frames/LagrangianFrame.hpp b/backend/src/Simulator/Frames/LagrangianFrame.hpp new file mode 100644 index 000000000..d6976364a --- /dev/null +++ b/backend/src/Simulator/Frames/LagrangianFrame.hpp @@ -0,0 +1,103 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +#include + +#include "Simulator/Frames/RotationConvention.hpp" + +#include "Utilities/NonCreatable.hpp" + +namespace elastica { + + namespace detail { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Trait defining unique IDs for different directions + * \ingroup domain + * + * LagrangianFrame is a trait defining the type and number of directions in + * 3D, in a convected frame. It serves as a strong type for all Lagrangian + * directions within the \elastica library. + */ + struct LagrangianFrame : public NonCreatable { + public: + //**Type definitions****************************************************** + //! The type of direction ID + enum class DirectionType : std::uint8_t { d1 = 0, d2, d3, Count }; + //************************************************************************ + + //**Static members******************************************************** + //! unique ID for D1 direction + static constexpr DirectionType D1 = DirectionType::d1; + //! unique ID for D2 direction + static constexpr DirectionType D2 = DirectionType::d2; + //! unique ID for D3 direction + static constexpr DirectionType D3 = DirectionType::d3; + //! unique ID for the normal direction + static constexpr DirectionType Normal = + DirectionType(RotationConvention::Normal); + //! unique ID for the binormal direction + static constexpr DirectionType Binormal = + DirectionType(RotationConvention::Binormal); + //! unique ID for the tangent direction + static constexpr DirectionType Tangent = + DirectionType(RotationConvention::Tangent); + //! Dimension of simulation + static constexpr auto Dimension = + static_cast(DirectionType::Count); + //************************************************************************ + }; + /*! \endcond */ + //************************************************************************** + + //************************************************************************** + /*!\brief Casts a Lagrangian frame direction into an index value + * \ingroup simulator + * + * \details + * Converts a LagrangianFrame::DirectionType into an index value for use in + * generic code or associative arrays. + * + * \example + * \snippet Test_LagrangianFrame.cpp cast_along_eg + * + * \param dir A Direction to cast as index to + */ + inline constexpr auto cast_along(LagrangianFrame::DirectionType dir) + -> std::size_t { + return static_cast(dir); + } + //************************************************************************** + + //************************************************************************** + /*!\brief Casts a Lagrangian frame direction into an index value + * \ingroup simulator + * + * \details + * Converts a LagrangianFrame::DirectionType into an index value for use in + * generic code or associative arrays. + * + * \example + * \snippet Test_LagrangianFrame.cpp cast_along_template_eg + * + * \tparam Dir A Direction to cast as index to + */ + template + inline constexpr auto cast_along() -> std::size_t { + return cast_along(Dir); + } + //************************************************************************** + + } // namespace detail + +} // namespace elastica diff --git a/backend/src/Simulator/Frames/RotationConvention.hpp b/backend/src/Simulator/Frames/RotationConvention.hpp new file mode 100644 index 000000000..724dfb3f9 --- /dev/null +++ b/backend/src/Simulator/Frames/RotationConvention.hpp @@ -0,0 +1,27 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include + +namespace elastica { + + //************************************************************************** + /*!\brief The rotation convention followed by \elastica + * \ingroup utils + * + * \details + * TODO + */ + struct RotationConvention { + //! unique ID for the normal direction + static constexpr std::uint8_t Normal = 0U; + //! unique ID for the binormal direction + static constexpr std::uint8_t Binormal = 1U; + //! unique ID for the tangent direction + static constexpr std::uint8_t Tangent = 2U; + }; + //**************************************************************************** + +} // namespace elastica diff --git a/backend/src/Utilities/Math/Types.hpp b/backend/src/Utilities/Math/Types.hpp index 6c84669e7..3357c1774 100644 --- a/backend/src/Utilities/Math/Types.hpp +++ b/backend/src/Utilities/Math/Types.hpp @@ -1,16 +1,3 @@ -//============================================================================== -/*! -// \file -// \brief Types belonging to the Math module -// -// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved -// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved -// -// Distributed under the MIT License. -// See LICENSE.txt for details. -*/ -//============================================================================== - #pragma once //****************************************************************************** diff --git a/backend/src/Utilities/Math/Vec3.hpp b/backend/src/Utilities/Math/Vec3.hpp index f73a7d5be..f1bf13b31 100644 --- a/backend/src/Utilities/Math/Vec3.hpp +++ b/backend/src/Utilities/Math/Vec3.hpp @@ -1,16 +1,3 @@ -//============================================================================== -/*! -// \file -// \brief Header for a 3D vector for use in \elastica interface -// -// Copyright (C) 2020-2020 Tejaswin Parthasarathy - All Rights Reserved -// Copyright (C) 2020-2020 MattiaLab - All Rights Reserved -// -// Distributed under the MIT License. -// See LICENSE.txt for details. -*/ -//============================================================================== - #pragma once //****************************************************************************** diff --git a/backend/src/Utilities/NonCreatable.hpp b/backend/src/Utilities/NonCreatable.hpp new file mode 100644 index 000000000..ca2caa549 --- /dev/null +++ b/backend/src/Utilities/NonCreatable.hpp @@ -0,0 +1,77 @@ +//============================================================================== +/*! + * From: file pe/util/NonCreatable.h + * \brief Base class for for non-creatable (static) classes + * + * Copyright (C) 2009 Klaus Iglberger + * + * This file is part of pe. + * + * pe is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * pe is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along with + * pe. If not, see . + */ +//============================================================================== + +#pragma once + +namespace elastica { + + //============================================================================ + // + // CLASS DEFINITION + // + //============================================================================ + + //**************************************************************************** + /*!\brief Base class for non-creatable (static) classes. + // \ingroup util + // + // The NonCreatable class is intended to work as a base class for + // non-creatable classes, i.e. classes that cannot be instantiated and + // exclusively offer static functions/data. Both the standard as well as the + // copy constructor and the copy assignment operator are deleted + // and left undefined in order to prohibit the instantiation of objects of + // derived classes. + // + // \note + // It is not necessary to publicly derive from this class. It is sufficient to + // derive privately to prevent the instantiation of the derived class. + // + // \example + \code + class A : private NonCreatable + { ... }; + \endcode + // + // \see NonCopyable + */ + class NonCreatable { + public: + //**Constructors and copy assignment operator******************************* + /*!\name Constructors and copy assignment operator */ + //@{ + //! Constructor (private & deleted) + NonCreatable() = delete; + //! Copy constructor (private & deleted) + NonCreatable(const NonCreatable&) = delete; + //! Move constructor (private & deleted) + NonCreatable(NonCreatable&&) = delete; + //! Copy assignment operator (private & deleted) + NonCreatable& operator=(const NonCreatable&) = delete; + //! Move assignment operator (private & deleted) + NonCreatable& operator=(NonCreatable&&) = delete; + //@} + //************************************************************************** + }; + //**************************************************************************** + +} // namespace elastica From d03eb976e1516744c44fed5d74166bf0d3bfdb65 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 10 Jul 2024 23:13:09 -0500 Subject: [PATCH 19/62] binding: elastica tags --- backend/src/Systems/Block/TypeTraits.hpp | 26 ++ backend/src/Systems/Block/Types.hpp | 10 + backend/src/Systems/Python/BindTags.cpp | 27 ++ backend/src/Systems/Python/BindTags.hpp | 46 +++ backend/src/Systems/Python/Bindings.cpp | 17 + backend/src/Systems/Systems.hpp | 60 +++ backend/src/Systems/Types.hpp | 10 + backend/src/Utilities/Math/Zero.hpp | 50 +++ backend/src/Utilities/TMPL.hpp | 444 +++++++++++++++++++++++ 9 files changed, 690 insertions(+) create mode 100644 backend/src/Systems/Block/TypeTraits.hpp create mode 100644 backend/src/Systems/Block/Types.hpp create mode 100644 backend/src/Systems/Python/BindTags.cpp create mode 100644 backend/src/Systems/Python/BindTags.hpp create mode 100644 backend/src/Systems/Python/Bindings.cpp create mode 100644 backend/src/Systems/Systems.hpp create mode 100644 backend/src/Systems/Types.hpp create mode 100644 backend/src/Utilities/Math/Zero.hpp create mode 100644 backend/src/Utilities/TMPL.hpp diff --git a/backend/src/Systems/Block/TypeTraits.hpp b/backend/src/Systems/Block/TypeTraits.hpp new file mode 100644 index 000000000..62da94b0d --- /dev/null +++ b/backend/src/Systems/Block/TypeTraits.hpp @@ -0,0 +1,26 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +/// +#include "Types.hpp" +/// +#include "Systems/Block/Block/Aliases.hpp" +#include "Systems/Block/Block/TypeTraits.hpp" +#include "Systems/Block/BlockVariables/Aliases.hpp" +#include "Systems/Block/BlockVariables/TypeTraits.hpp" + +//============================================================================== +// +// DOXYGEN DOCUMENTATION +// +//============================================================================== + +//****************************************************************************** +/*!\defgroup block_tt Block Type Traits + * \ingroup blocks + * \brief Type traits for blocks + */ +//****************************************************************************** diff --git a/backend/src/Systems/Block/Types.hpp b/backend/src/Systems/Block/Types.hpp new file mode 100644 index 000000000..baa7f944e --- /dev/null +++ b/backend/src/Systems/Block/Types.hpp @@ -0,0 +1,10 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** + +/// +#include "Systems/Block/Block/Types.hpp" +#include "Systems/Block/BlockVariables/Types.hpp" +/// diff --git a/backend/src/Systems/Python/BindTags.cpp b/backend/src/Systems/Python/BindTags.cpp new file mode 100644 index 000000000..3a56bbeda --- /dev/null +++ b/backend/src/Systems/Python/BindTags.cpp @@ -0,0 +1,27 @@ +//****************************************************************************** +// Includes +//****************************************************************************** + +// +#include +// +#include "Systems/Systems.hpp" +// +#include "Systems/Python/BindTags.hpp" +// +#include "Utilities/TMPL.hpp" + +namespace py = pybind11; + +namespace py_bindings { + + //**************************************************************************** + /*!\brief Helps bind tags for all physical systems in \elastica + * \ingroup python_bindings + */ + void bind_tags(py::module& m) { // NOLINT + bind_tags<::elastica::PhysicalSystemPlugins>(m); + } + //**************************************************************************** + +} // namespace py_bindings diff --git a/backend/src/Systems/Python/BindTags.hpp b/backend/src/Systems/Python/BindTags.hpp new file mode 100644 index 000000000..43ebca155 --- /dev/null +++ b/backend/src/Systems/Python/BindTags.hpp @@ -0,0 +1,46 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +// +#include "Systems/Block/TypeTraits.hpp" +// +// #include "Utilities/PrettyType.hpp" +// +#include +// +#include + +namespace py_bindings { + + //**************************************************************************** + /*!\brief Helps bind (unique) tags of Plugins in \elastica + * \ingroup python_bindings + * + * \details + * + * + * \param m A python module + */ + template + void bind_tags(pybind11::module& m) { + using AllVariables = tmpl::flatten< + tmpl::transform>>; + using UniqueTags = tmpl::remove_duplicates>>; + + // first define variables as a read only property for named evaluation + tmpl::for_each([&](auto v) { + namespace py = pybind11; + using Tag = tmpl::type_from; + + // FIXME: avoiding pretty_type for now. + // std::string tag_name = pretty_type::short_name(); + std::string tag_name = typeid(v).name(); + py::class_(m, ("_" + tag_name).c_str(), + ("Symbol corresponding to C++ tag " + tag_name).c_str()); + }); + } + +} // namespace py_bindings diff --git a/backend/src/Systems/Python/Bindings.cpp b/backend/src/Systems/Python/Bindings.cpp new file mode 100644 index 000000000..de7d08ac4 --- /dev/null +++ b/backend/src/Systems/Python/Bindings.cpp @@ -0,0 +1,17 @@ +//****************************************************************************** +// Includes +//****************************************************************************** +#include + +namespace py = pybind11; + +namespace py_bindings { + void bind_tags(py::module& m); // NOLINT +} // namespace py_bindings + +PYBIND11_MODULE(_PyTags, m) { // NOLINT + m.doc() = R"pbdoc( + Bindings for Elastica++ tag types + )pbdoc"; + py_bindings::bind_tags(m); +} diff --git a/backend/src/Systems/Systems.hpp b/backend/src/Systems/Systems.hpp new file mode 100644 index 000000000..c2691208b --- /dev/null +++ b/backend/src/Systems/Systems.hpp @@ -0,0 +1,60 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +// +#include "Systems/Types.hpp" +// +#include "Systems/Protocols.hpp" +#include "Systems/Tags.hpp" +// +#include "Systems/Block.hpp" +#include "Systems/CosseratRods.hpp" +// #include "Systems/RigidBody.hpp" +#include "Systems/States/States.hpp" +// +#include "Utilities/TMPL.hpp" + +//============================================================================== +// +// DOXYGEN DOCUMENTATION +// +//============================================================================== + +//****************************************************************************** +/*!\defgroup systems Physical systems + * \brief Physical systems (that occupy space and evolve in time) in \elastica + */ +//****************************************************************************** + +namespace elastica { + + //**************************************************************************** + /*!\brief All implemented CosseratRod plugins in \elastica + * \ingroup systems + */ + using CosseratRodPlugins = + tmpl::list<::elastica::cosserat_rod::CosseratRod, + ::elastica::cosserat_rod::CosseratRodWithoutDamping>; + //**************************************************************************** + + //**************************************************************************** + /*!\brief All implemented RigidBody plugins in \elastica + * \ingroup systems + */ + // using RigidBodyPlugins = tmpl::list<::elastica::rigid_body::Sphere>; + //**************************************************************************** + + //**************************************************************************** + /*!\brief All implemented physical system plugins in \elastica + * \ingroup systems + * + * \see elastica::CosseratRodPlugins, elastica::RigidBodyPlugins + */ + using PhysicalSystemPlugins = + tmpl::append; + // tmpl::append; + //**************************************************************************** + +} // namespace elastica diff --git a/backend/src/Systems/Types.hpp b/backend/src/Systems/Types.hpp new file mode 100644 index 000000000..35695934a --- /dev/null +++ b/backend/src/Systems/Types.hpp @@ -0,0 +1,10 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include "Systems/Block/Types.hpp" +#include "Systems/CosseratRods/Types.hpp" +// #include "Systems/RigidBody/Types.hpp" +#include "Systems/States/Types.hpp" +#include "Systems/common/Types.hpp" diff --git a/backend/src/Utilities/Math/Zero.hpp b/backend/src/Utilities/Math/Zero.hpp new file mode 100644 index 000000000..3aa3eea19 --- /dev/null +++ b/backend/src/Utilities/Math/Zero.hpp @@ -0,0 +1,50 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +#include + +#include "Utilities/Math/Vec3.hpp" + +namespace elastica { + + //**************************************************************************** + /*!\brief Checks for a zero number + */ + inline auto is_zero(real_t real_number) noexcept -> bool { + return ::blaze::isZero(real_number); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Checks for a zero vector + */ + inline auto is_zero(const Vec3& vec) noexcept -> bool { + return ::blaze::isZero(vec); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Checks for a zero length vector + * + * Checks if the length of a vector is zero or as close to zero that it + * can not be distinguished form zero + */ + inline auto is_zero_length(const Vec3& vec) noexcept -> bool { + // return vec.sqrLength() < Limits::fpuAccuracy(); + return ::blaze::sqrLength(vec) < ::blaze::accuracy; + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Returns a zero vector + */ + inline auto zero3() noexcept -> Vec3 { + return Vec3{static_cast(0.0), static_cast(0.0), + static_cast(0.0)}; + } + //**************************************************************************** + +} // namespace elastica diff --git a/backend/src/Utilities/TMPL.hpp b/backend/src/Utilities/TMPL.hpp new file mode 100644 index 000000000..8898752c1 --- /dev/null +++ b/backend/src/Utilities/TMPL.hpp @@ -0,0 +1,444 @@ +#pragma once + +// Since this header only wraps brigand and several additions to it we mark +// it as a system header file so that clang-tidy ignores it. +#ifdef __GNUC__ +#pragma GCC system_header +#endif + +#define BRIGAND_NO_BOOST_SUPPORT +#include +#include + +#include "Requires.hpp" +//#include "TypeTraits.hpp" + +namespace tmpl = brigand; + +// Use spectre's TMPL + +namespace brigand { + /// Check if a typelist contains an item. + template + using list_contains = + tmpl::found>>; + + template + constexpr const bool list_contains_v = list_contains::value; + + /// Obtain the elements of `Sequence1` that are not in `Sequence2`. + template + using list_difference = + fold>; + +} // namespace brigand + +namespace brigand { + namespace detail { + template + struct remove_duplicates_helper { + using type = typename std::conditional< + std::is_same, no_such_type_>::value, push_back, + S>::type; + }; + } // namespace detail + + template + using remove_duplicates = + fold, detail::remove_duplicates_helper<_state, _element>>; +} // namespace brigand + +namespace brigand { + template + using cartesian_product = reverse_fold< + list, list>, + lazy::join, + defer>>>>>>>>>; + +} + +namespace brigand { + template + struct conditional; + + template <> + struct conditional { + template + using type = T; + }; + + template <> + struct conditional { + template + using type = F; + }; + + template + using conditional_t = typename conditional::template type; +} // namespace brigand + +/*! + * \ingroup UtilitiesGroup + * \brief Allows zero-cost unordered expansion of a parameter + * + * \details + * Expands a parameter pack, typically useful for runtime evaluation via a + * Callable such as a lambda, function, or function object. For example, + * an unordered transform of a std::tuple can be implemented as: + * \snippet Utilities/Test_TMPL.cpp expand_pack_example + * + * \see tuple_fold tuple_counted_fold tuple_transform std::tuple + * EXPAND_PACK_LEFT_TO_RIGHT + */ +template +constexpr void expand_pack(Ts&&...) noexcept {} + +/*! + * \ingroup UtilitiesGroup + * \brief Expand a parameter pack evaluating the terms from left to right. + * + * The parameter pack inside the argument to the macro must not be expanded + * since the macro will do the expansion correctly for you. In the below example + * a parameter pack of `std::integral_constant` is passed to the + * function. The closure `lambda` is used to sum up the values of all the `Ts`. + * Note that the `Ts` passed to `EXPAND_PACK_LEFT_TO_RIGHT` is not expanded. + * + * \snippet Utilities/Test_TMPL.cpp expand_pack_left_to_right + * + * \see tuple_fold tuple_counted_fold tuple_transform std::tuple expand_pack + */ +#define EXPAND_PACK_LEFT_TO_RIGHT(...) \ + (void)std::initializer_list { ((void)(__VA_ARGS__), '0')... } + +/*! + * \ingroup UtilitiesGroup + * \brief Returns the first argument of a parameter pack + */ +template +constexpr decltype(auto) get_first_argument(T&& t, Ts&&... /*rest*/) noexcept { + return t; +} + +// namespace tmpl { +// // Types list +// template +// struct list {}; +// +// // Values list +// template +// using integral_list = list...>; +// +// // Empty list +// using empty_sequence = tmpl::list<>; +// +// // Predefined placeholders +// struct _1 {}; +// struct _2 {}; +// +// // Utility to expand a typelist into its variadic constituents +// // and apply a functor/expression on the pack +// namespace detail { +// template class B> +// struct wrap; +// +// template