From 4fe224b5c079b61a8037fe1b7aeb7bf70b20036c Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Tue, 16 Jan 2024 12:36:19 -0600 Subject: [PATCH 1/4] Remove Eigen support, add xtensor We are planning to replace Eigen with xtensor. This is a first step. It builds successfully, but the tests will not pass, because the inverse distortion function needs to be re-written to use xtensor instead of Eigen. Signed-off-by: Patrick Avery --- conda.recipe/meta.yaml | 2 +- .../cpp_sublibrary/src/inverse_distortion.cpp | 19 ++------- scripts/install/install_build_dependencies.py | 42 ++++++++++++++++++- setup.py | 5 ++- 4 files changed, 48 insertions(+), 20 deletions(-) diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index 94bb3fe9a..211273269 100644 --- a/conda.recipe/meta.yaml +++ b/conda.recipe/meta.yaml @@ -30,7 +30,7 @@ requirements: - numba - pybind11 - xsimd - - eigen + - xtensor run: - appdirs diff --git a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp index 21a0001b5..2a0ab7b9f 100644 --- a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp +++ b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp @@ -1,26 +1,13 @@ #include -#include -#include +#include namespace py = pybind11; constexpr double FOUR_THIRDS_PI = 4.1887902; constexpr double N_THREE_HALVES_SQRT_3 = -2.59807621; constexpr double TWO_OVER_SQRT_THREE = 1.154700538; -Eigen::ArrayXXd ge_41rt_inverse_distortion(const Eigen::ArrayXXd& inputs, const double rhoMax, const Eigen::ArrayXd& params) { - Eigen::ArrayXd radii = inputs.matrix().rowwise().norm(); - if (radii.maxCoeff() < 1e-10) { - return Eigen::ArrayXXd::Zero(inputs.rows(), inputs.cols()); - } - Eigen::ArrayXd inverted_radii = radii.cwiseInverse(); - Eigen::ArrayXd cosines = inputs.col(0) * inverted_radii; - Eigen::ArrayXd cosine_double_angles = 2*cosines.square() - 1; - Eigen::ArrayXd cosine_quadruple_angles = 2*cosine_double_angles.square() - 1; - Eigen::ArrayXd sqrt_p_is = rhoMax / (-params[0]*cosine_double_angles - params[1]*cosine_quadruple_angles - params[2]).sqrt(); - Eigen::ArrayXd solutions = TWO_OVER_SQRT_THREE*sqrt_p_is*(acos(N_THREE_HALVES_SQRT_3*radii/sqrt_p_is)/3 + FOUR_THIRDS_PI).cos(); - Eigen::ArrayXXd results = solutions.rowwise().replicate(inputs.cols()).array() * inputs * inverted_radii.rowwise().replicate(inputs.cols()).array(); - - return results; +void ge_41rt_inverse_distortion() { + return; } PYBIND11_MODULE(inverse_distortion, m) diff --git a/scripts/install/install_build_dependencies.py b/scripts/install/install_build_dependencies.py index 2e35e2de4..451efc189 100755 --- a/scripts/install/install_build_dependencies.py +++ b/scripts/install/install_build_dependencies.py @@ -36,6 +36,44 @@ def download_and_extract_tgz(url, md5sum, path): temp_file.unlink() +def download_xtensor(path): + url = 'https://github.com/xtensor-stack/xtensor/archive/refs/tags/0.24.7.tar.gz' # noqa + md5sum = 'e21a14d679db71e92a703bccd3c5866a' + out_dir_name = 'xtensor-0.24.7' + + with tempfile.TemporaryDirectory() as temp_dir: + download_and_extract_tgz(url, md5sum, temp_dir) + + target_path = Path(path) / 'xtensor' + if target_path.exists(): + shutil.rmtree(target_path) + + os.makedirs(path, exist_ok=True) + shutil.move(str(Path(temp_dir) / out_dir_name / 'include/xtensor'), + str(Path(path) / 'xtensor/xtensor')) + + return str(target_path) + + +def download_xtl(path): + url = 'https://github.com/xtensor-stack/xtl/archive/refs/tags/0.7.7.tar.gz' # noqa + md5sum = '6df56ae8bc30471f6773b3f18642c8ab' + out_dir_name = 'xtl-0.7.7' + + with tempfile.TemporaryDirectory() as temp_dir: + download_and_extract_tgz(url, md5sum, temp_dir) + + target_path = Path(path) / 'xtl' + if target_path.exists(): + shutil.rmtree(target_path) + + os.makedirs(path, exist_ok=True) + shutil.move(str(Path(temp_dir) / out_dir_name / 'include/xtl'), + str(Path(path) / 'xtl/xtl')) + + return str(target_path) + + def download_xsimd(path): url = 'https://github.com/xtensor-stack/xsimd/archive/refs/tags/7.6.0.tar.gz' # noqa md5sum = '6e52c2af8b3cb4688993a0e70825f4e8' @@ -93,8 +131,10 @@ def download_pybind11(path): INSTALL_FUNCTIONS = { 'eigen3': download_eigen3, - 'xsimd': download_xsimd, 'pybind11': download_pybind11, + 'xsimd': download_xsimd, + 'xtensor': download_xtensor, + 'xtl': download_xtl, } diff --git a/setup.py b/setup.py index addf3d567..efefc4c4e 100644 --- a/setup.py +++ b/setup.py @@ -115,15 +115,16 @@ def get_cpp_extensions(): cpp_transform_pkgdir = Path('hexrd') / 'transforms/cpp_sublibrary' src_files = [str(cpp_transform_pkgdir / 'src/inverse_distortion.cpp')] - extra_compile_args = ['-O3', '-Wall', '-shared', '-std=c++11', + extra_compile_args = ['-O3', '-Wall', '-shared', '-std=c++14', '-funroll-loops'] if not sys.platform.startswith('win'): extra_compile_args.append('-fPIC') # Define include directories include_dirs = [ - get_include_path('eigen3'), get_include_path('xsimd'), + get_include_path('xtensor'), + get_include_path('xtl'), get_pybind11_include_path(), numpy.get_include(), ] From 7f56e54e60bb07b2748ee8feaacc8c78e56396ff Mon Sep 17 00:00:00 2001 From: Zack Singer Date: Wed, 17 Jan 2024 22:16:27 -0500 Subject: [PATCH 2/4] Convert Eigen -> xTensor. Array dims incorrect sometimes? --- hexrd/transforms/cpp_sublibrary/Makefile | 4 +-- .../cpp_sublibrary/src/inverse_distortion.cpp | 28 ++++++++++++++++--- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/hexrd/transforms/cpp_sublibrary/Makefile b/hexrd/transforms/cpp_sublibrary/Makefile index 7ce236b58..fb38bdcd7 100644 --- a/hexrd/transforms/cpp_sublibrary/Makefile +++ b/hexrd/transforms/cpp_sublibrary/Makefile @@ -1,6 +1,6 @@ CXX = c++ -CXXFLAGS = -O3 -Wall -shared -fPIC -funroll-loops -Wno-maybe-uninitialized -fopenmp -INCLUDES = -I/${CONDA_PREFIX}/include/eigen3/ -I/${CONDA_PREFIX}/include/xsimd -I/${CONDA_PREFIX}/include -I/${CONDA_PREFIX}/include/python3.11/ +CXXFLAGS = -O3 -Wall -shared -fPIC -funroll-loops -Wno-maybe-uninitialized -Wno-attributes -fopenmp +INCLUDES = -I/${CONDA_PREFIX}/include -I/${CONDA_PREFIX}/include/python3.11/ -I/${CONDA_PREFIX}/lib/python3.11/site-packages/numpy/core/include SRCS_INVERSE_DISTORTION = src/inverse_distortion.cpp TARGET_DIR = ../../extensions/ TARGET_NAME_INVERSE_DISTORTION = inverse_distortion diff --git a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp index 2a0ab7b9f..bebc439c4 100644 --- a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp +++ b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp @@ -1,18 +1,38 @@ +#define FORCE_IMPORT_ARRAY + #include +#include #include +#include +#include namespace py = pybind11; constexpr double FOUR_THIRDS_PI = 4.1887902; constexpr double N_THREE_HALVES_SQRT_3 = -2.59807621; constexpr double TWO_OVER_SQRT_THREE = 1.154700538; -void ge_41rt_inverse_distortion() { - return; +xt::pyarray ge_41rt_inverse_distortion(const xt::pyarray& inputs, const double rhoMax, const xt::pyarray& params) { + auto radii = xt::sqrt(xt::sum(xt::square(inputs), {1})); + + if (xt::amax(radii)() < 1e-10) { + return xt::zeros({inputs.shape()[0], inputs.shape()[1]}); + } + + auto inverted_radii = xt::pow(radii, -1); + xt::pyarray cosines = inputs * xt::view(inverted_radii, xt::newaxis()); + auto cosine_double_angles = 2 * xt::square(cosines) - 1; + auto cosine_quadruple_angles = 2 * xt::square(cosine_double_angles) - 1; + auto sqrt_p_is = rhoMax / xt::sqrt(-params[0] * cosine_double_angles - params[1] * cosine_quadruple_angles - params[2]); + auto solutions = TWO_OVER_SQRT_THREE * sqrt_p_is * xt::cos(xt::acos(N_THREE_HALVES_SQRT_3 * radii / sqrt_p_is) / 3 + FOUR_THIRDS_PI); + xt::pyarray results = solutions * inputs * xt::view(inverted_radii, xt::newaxis()); + + return results; } PYBIND11_MODULE(inverse_distortion, m) { - m.doc() = "HEXRD inverse distribution plugin"; + xt::import_numpy(); + m.doc() = "HEXRD inverse distribution plugin"; - m.def("ge_41rt_inverse_distortion", &ge_41rt_inverse_distortion, "Inverse distortion for ge_41rt"); + m.def("ge_41rt_inverse_distortion", &ge_41rt_inverse_distortion, "Inverse distortion for ge_41rt"); } From 9caa1eea14ae23721717300e9cc8c0a0420ef536 Mon Sep 17 00:00:00 2001 From: Zack Singer Date: Fri, 19 Jan 2024 11:39:03 -0500 Subject: [PATCH 3/4] Resolve matrix output dimensionality issues --- hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp index bebc439c4..799f9dfff 100644 --- a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp +++ b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp @@ -13,18 +13,18 @@ constexpr double TWO_OVER_SQRT_THREE = 1.154700538; xt::pyarray ge_41rt_inverse_distortion(const xt::pyarray& inputs, const double rhoMax, const xt::pyarray& params) { auto radii = xt::sqrt(xt::sum(xt::square(inputs), {1})); - + if (xt::amax(radii)() < 1e-10) { return xt::zeros({inputs.shape()[0], inputs.shape()[1]}); } auto inverted_radii = xt::pow(radii, -1); - xt::pyarray cosines = inputs * xt::view(inverted_radii, xt::newaxis()); + xt::pyarray cosines = xt::view(inputs, xt::all(), 0) * inverted_radii; auto cosine_double_angles = 2 * xt::square(cosines) - 1; auto cosine_quadruple_angles = 2 * xt::square(cosine_double_angles) - 1; auto sqrt_p_is = rhoMax / xt::sqrt(-params[0] * cosine_double_angles - params[1] * cosine_quadruple_angles - params[2]); auto solutions = TWO_OVER_SQRT_THREE * sqrt_p_is * xt::cos(xt::acos(N_THREE_HALVES_SQRT_3 * radii / sqrt_p_is) / 3 + FOUR_THIRDS_PI); - xt::pyarray results = solutions * inputs * xt::view(inverted_radii, xt::newaxis()); + xt::pyarray results = xt::view(solutions, xt::all(), xt::newaxis()) * inputs * xt::view(inverted_radii, xt::all(), xt::newaxis()); return results; } From 475987fd0718aac3acab43e1adfee6750cdf9ff0 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Fri, 19 Jan 2024 11:32:31 -0600 Subject: [PATCH 4/4] Install xtensor-python and upgrade xsimd This was required for the compilation with xsimd to succeed. Signed-off-by: Patrick Avery --- conda.recipe/meta.yaml | 1 + .../cpp_sublibrary/src/inverse_distortion.cpp | 1 + scripts/install/install_build_dependencies.py | 28 ++++++++++++++++--- setup.py | 1 + 4 files changed, 27 insertions(+), 4 deletions(-) diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index 211273269..9a03d7c34 100644 --- a/conda.recipe/meta.yaml +++ b/conda.recipe/meta.yaml @@ -31,6 +31,7 @@ requirements: - pybind11 - xsimd - xtensor + - xtensor-python run: - appdirs diff --git a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp index 799f9dfff..196c31394 100644 --- a/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp +++ b/hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp @@ -1,4 +1,5 @@ #define FORCE_IMPORT_ARRAY +#define XTENSOR_USE_XSIMD #include #include diff --git a/scripts/install/install_build_dependencies.py b/scripts/install/install_build_dependencies.py index 451efc189..abc5299f0 100755 --- a/scripts/install/install_build_dependencies.py +++ b/scripts/install/install_build_dependencies.py @@ -55,6 +55,26 @@ def download_xtensor(path): return str(target_path) +def download_xtensor_python(path): + url = 'https://github.com/xtensor-stack/xtensor-python/archive/refs/tags/0.26.1.tar.gz' # noqa + md5sum = '5d05edf71ac948dc620968229320c791' + out_dir_name = 'xtensor-python-0.26.1' + + with tempfile.TemporaryDirectory() as temp_dir: + download_and_extract_tgz(url, md5sum, temp_dir) + + target_path = Path(path) / 'xtensor-python' + if target_path.exists(): + shutil.rmtree(target_path) + + os.makedirs(path, exist_ok=True) + shutil.move( + str(Path(temp_dir) / out_dir_name / 'include/xtensor-python'), + str(Path(path) / 'xtensor-python/xtensor-python')) + + return str(target_path) + + def download_xtl(path): url = 'https://github.com/xtensor-stack/xtl/archive/refs/tags/0.7.7.tar.gz' # noqa md5sum = '6df56ae8bc30471f6773b3f18642c8ab' @@ -75,9 +95,9 @@ def download_xtl(path): def download_xsimd(path): - url = 'https://github.com/xtensor-stack/xsimd/archive/refs/tags/7.6.0.tar.gz' # noqa - md5sum = '6e52c2af8b3cb4688993a0e70825f4e8' - out_dir_name = 'xsimd-7.6.0' + url = 'https://github.com/xtensor-stack/xsimd/archive/refs/tags/12.1.1.tar.gz' # noqa + md5sum = 'e8887de343bd6036bdfa1f4a4752dc64' + out_dir_name = 'xsimd-12.1.1' with tempfile.TemporaryDirectory() as temp_dir: download_and_extract_tgz(url, md5sum, temp_dir) @@ -125,7 +145,6 @@ def download_pybind11(path): shutil.move(str(Path(temp_dir) / out_dir_name / 'include/pybind11'), str(Path(path) / 'pybind11/pybind11')) - return str(target_path) @@ -134,6 +153,7 @@ def download_pybind11(path): 'pybind11': download_pybind11, 'xsimd': download_xsimd, 'xtensor': download_xtensor, + 'xtensor-python': download_xtensor_python, 'xtl': download_xtl, } diff --git a/setup.py b/setup.py index efefc4c4e..5707bf151 100644 --- a/setup.py +++ b/setup.py @@ -124,6 +124,7 @@ def get_cpp_extensions(): include_dirs = [ get_include_path('xsimd'), get_include_path('xtensor'), + get_include_path('xtensor-python'), get_include_path('xtl'), get_pybind11_include_path(), numpy.get_include(),