From 0e24abc5c12bd74f1359cefd62405c3af175bf80 Mon Sep 17 00:00:00 2001 From: Antony Chan Date: Wed, 27 Nov 2024 14:28:38 -0800 Subject: [PATCH] Option to use Eigen versus Lapack to solve equations On Windows system, it is often tedious to install Lapack/Atlas/Openblas packages so as to simulate PSFs. Configure the project: ``` meson setup -Duse_eigen=true -Darmadillo-code:lapack=none build-windows/ ``` In order to bypass LAPACK depdendencies. Download Eigen3 the header-only library. Implement the Bidiagonal divide and conquer SVD solver. Write a wrapper around the Eigen logic to avoid C++ namespace pollusion. --- .github/workflows/validate-build.yml | 19 ++++--- .../packagefiles/armadillo-code/meson.build | 5 +- .../packagefiles/armadillo-code/meson.options | 1 + besselj_armadillo_support/meson.build | 3 +- examples/generate-psf.cpp | 9 ++++ examples/meson.build | 5 +- meson.options | 11 +++- microsc-psf/inc/linsolver.h | 23 ++++++++ microsc-psf/inc/make_psf.h | 10 +++- microsc-psf/meson.build | 35 ++++++++++++ microsc-psf/src/linsolver.cpp | 53 +++++++++++++++++++ microsc-psf/src/main.cpp | 39 ++++++++++++-- 12 files changed, 194 insertions(+), 19 deletions(-) create mode 100644 microsc-psf/inc/linsolver.h create mode 100644 microsc-psf/src/linsolver.cpp diff --git a/.github/workflows/validate-build.yml b/.github/workflows/validate-build.yml index aa31738..545712d 100644 --- a/.github/workflows/validate-build.yml +++ b/.github/workflows/validate-build.yml @@ -16,7 +16,7 @@ jobs: strategy: matrix: #os: [ubuntu-20.04, ubuntu-latest, windows-latest, macos-latest] - os: [ubuntu-20.04, ubuntu-latest] + os: [ubuntu-20.04, ubuntu-latest, windows-latest] use_boost: ["false", "true"] steps: @@ -54,14 +54,17 @@ jobs: build/ if: startsWith(matrix.os, 'window') == false - #- name: Resolve C++ build dependencies (msvc toolchain) - # run: meson setup --vsenv build/ - # if: startsWith(matrix.os, 'window') == true + - name: Resolve C++ build dependencies (msvc toolchain) + run: >- + meson setup + --vsenv + -Darmadillo-code:lapack=none + -Duse_eigen=true + build/ + if: startsWith(matrix.os, 'window') == true - name: Build everything - #run: meson compile -C build/ - run: ninja -C build/ all + run: meson compile -C build/ - name: Test everything - #run: meson test -C build/ - run: ninja -C build/ test + run: meson test -C build/ diff --git a/3rdparty/packagefiles/armadillo-code/meson.build b/3rdparty/packagefiles/armadillo-code/meson.build index 46c8b42..4b6f415 100644 --- a/3rdparty/packagefiles/armadillo-code/meson.build +++ b/3rdparty/packagefiles/armadillo-code/meson.build @@ -5,7 +5,10 @@ project('armadillo', 'cpp', # Header-only library armadillo_inc = include_directories('include') -if get_option('lapack') == 'openblas64' +if get_option('lapack') == 'none' + warning('Disabled Lapack dependencies. SVD and LU algorithms will be undefined.') + lapack_dep = dependency('', required: false) +elif get_option('lapack') == 'openblas64' lapack_dep = dependency('openblas64') elif get_option('lapack') == 'openblas' lapack_dep = dependency('openblas') diff --git a/3rdparty/packagefiles/armadillo-code/meson.options b/3rdparty/packagefiles/armadillo-code/meson.options index 017520c..36b741a 100644 --- a/3rdparty/packagefiles/armadillo-code/meson.options +++ b/3rdparty/packagefiles/armadillo-code/meson.options @@ -6,6 +6,7 @@ option('lapack', 'openblas', 'lapack64', 'lapack', + 'none', ], value: 'openblas', ) diff --git a/besselj_armadillo_support/meson.build b/besselj_armadillo_support/meson.build index e72a0cc..480d30c 100644 --- a/besselj_armadillo_support/meson.build +++ b/besselj_armadillo_support/meson.build @@ -1,5 +1,6 @@ armadillo_besselj_support_inc = include_directories('inc') -armadillo_dep = subproject('armadillo-code').get_variable('armadillo_dep') +armadillo_proj = subproject('armadillo-code') +armadillo_dep = armadillo_proj.get_variable('armadillo_dep') # TODO(Antony): use "cxx.has_function" when it supports std namespace. if get_option('use_boost') diff --git a/examples/generate-psf.cpp b/examples/generate-psf.cpp index c697dcd..c75f6db 100644 --- a/examples/generate-psf.cpp +++ b/examples/generate-psf.cpp @@ -2,6 +2,14 @@ #include "make_psf.h" +namespace { +#ifdef LSTSQ_USE_LAPACK +constexpr bool has_lapack = true; +#else +constexpr bool has_lapack = false; +#endif +} // namespace + int main() { using arma::hdf5_name; using microsc_psf::makePSF; @@ -20,6 +28,7 @@ int main() { precision_li2017_t precision{}; precision.num_basis = 153; precision.rho_samples = 1000; + precision.solver = has_lapack ? microsc_psf::SandersonAndCurtin2020 : microsc_psf::EigenBdcSVD; const auto psf = makePSF(params, {0.1_um, 0.25_um}, {256, 128}, 0.610_um, precision); diff --git a/examples/meson.build b/examples/meson.build index bb39fe2..0d471aa 100644 --- a/examples/meson.build +++ b/examples/meson.build @@ -12,7 +12,10 @@ endif generate_psf_exe = executable('generate-psf', sources: 'generate-psf.cpp', - cpp_args: generate_psf_compile_args, + cpp_args: [ + generate_psf_compile_args, + use_lapack_args, + ], dependencies: [ microsc_psf_dep, hdf5_dep, diff --git a/meson.options b/meson.options index e944008..9142b21 100644 --- a/meson.options +++ b/meson.options @@ -1,2 +1,11 @@ option('install_examples', type: 'boolean', value: false) -option('use_boost', type: 'boolean', value: false) +option('use_boost', + description: 'Compute bessel function of the first kind with Boost::Math', + type: 'boolean', + value: false, +) +option('use_eigen', + description: 'Solve linear least squares with Eigen3 the matrix library', + type: 'boolean', + value: false, +) diff --git a/microsc-psf/inc/linsolver.h b/microsc-psf/inc/linsolver.h new file mode 100644 index 0000000..476f5a3 --- /dev/null +++ b/microsc-psf/inc/linsolver.h @@ -0,0 +1,23 @@ +#pragma once + +#include + +namespace microsc_psf { +namespace internal { + +/** Solve min || A * x - b ||^2 for x , using Eigen3 solver. + * + * Typical ussage: simulate PSF on machines not having LAPACK support, e.g. + * Windows on ARM64 CPU. + * + * Note: This may seem contrived to have two near-identical C++ libraries + * serving the linear algebra code, but Armadillo has a nice Matlab-like syntax + * that simplifies the code review process a lot. + * + * @tparam tranpose_b True if the signal b needs a Hermitian transpose ahead of + * the least squares solver. + */ +template +arma::cx_mat solveWithEigen(const arma::mat& A, const arma::cx_mat& b); +} // namespace internal +} // namespace microsc_psf \ No newline at end of file diff --git a/microsc-psf/inc/make_psf.h b/microsc-psf/inc/make_psf.h index 20a606e..591f0e1 100644 --- a/microsc-psf/inc/make_psf.h +++ b/microsc-psf/inc/make_psf.h @@ -30,8 +30,14 @@ enum linear_solver_option_t : uint8_t { // https://arma.sourceforge.net/armadillo_solver_2020.pdf SandersonAndCurtin2020, - // Penrose pseudo inverse. - PenroseInverse + /** Penrose pseudo inverse. */ + PenroseInverse, + + /** Eigen's SVD solver. + * + * Reference: https://eigen.tuxfamily.org/dox/classEigen_1_1BDCSVD.html + */ + EigenBdcSVD, }; struct precision_li2017_t { diff --git a/microsc-psf/meson.build b/microsc-psf/meson.build index afcf89c..7ad8750 100644 --- a/microsc-psf/meson.build +++ b/microsc-psf/meson.build @@ -1,10 +1,44 @@ microsc_psf_inc = include_directories('inc') + +if get_option('use_eigen') + warning('Solving linear least squares with Eigen3, not Lapack') + + # Eigen3 is a direct competitor of Armadillo library. This may seem + # contrived to require two nearly identicial libraries to compute PSF, but + # Armadillo has a Matlab-like syntax that simplifies the code review process + # a lot. + eigen_solver_lib = static_library('eigen_solver', + sources: 'src/linsolver.cpp', + include_directories: include_directories('inc'), + dependencies: [ + dependency('eigen', fallback: ['eigen', 'eigen_dep']), + armadillo_dep, + ], + ) + + eigen_solver_dep = declare_dependency( + link_with: eigen_solver_lib, + include_directories: include_directories('inc'), + compile_args: '-DLSTSQ_USE_EIGEN', + ) +else + eigen_solver_dep = [] +endif + +lapack_dep = armadillo_proj.get_variable('lapack_dep') +if lapack_dep.found() + use_lapack_args = '-DLSTSQ_USE_LAPACK' +else + use_lapack_args = [] +endif + microsc_psf_lib = static_library('microsc-psf', sources: [ 'src/main.cpp', ], cpp_args: [ '-D' + bessel_source, + use_lapack_args, ], include_directories: [ microsc_psf_inc, @@ -13,6 +47,7 @@ microsc_psf_lib = static_library('microsc-psf', dependencies: [ armadillo_dep, boost_math_dep, + eigen_solver_dep, ], ) diff --git a/microsc-psf/src/linsolver.cpp b/microsc-psf/src/linsolver.cpp new file mode 100644 index 0000000..4224a73 --- /dev/null +++ b/microsc-psf/src/linsolver.cpp @@ -0,0 +1,53 @@ +#include "linsolver.h" + +#include +#include + +namespace microsc_psf { + +namespace internal { + +template +arma::cx_mat +solveWithEigen(const arma::mat& A_buffer, const arma::cx_mat& b_buffer) { + using arma::cx_double; + using Eigen::ComputeThinU; + using Eigen::ComputeThinV; + using Eigen::Map; + using Eigen::MatrixX; + using Eigen::MatrixXd; + + const auto m = A_buffer.n_rows; + const auto n = A_buffer.n_cols; + const auto k = (transpose_b ? b_buffer.n_rows : b_buffer.n_cols); + + if constexpr (transpose_b) { + assert(b_buffer.n_cols == m); + } else { + assert(b_buffer.n_rows == m); + } + + // First, map to Eigen's primary data structure. + const Map A(A_buffer.memptr(), m, n); + const Map> b(b_buffer.memptr(), b_buffer.n_rows, b_buffer.n_cols); + + // Allocate the result buffer + arma::cx_mat xopt_buffer(n, k); + + // Next, solve with Eigen's SVD solver + Map> xopt(xopt_buffer.memptr(), n, k); + if constexpr (transpose_b) { + // Eigen may have an efficient Hermitian transposed solver. Move the transpose step into + // Eigen. + xopt = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b.transpose()); + } else { + xopt = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b); + } + + return xopt_buffer; +} + +template arma::cx_mat solveWithEigen(const arma::mat&, const arma::cx_mat&); +template arma::cx_mat solveWithEigen(const arma::mat&, const arma::cx_mat&); +} // namespace internal +} // namespace microsc_psf diff --git a/microsc-psf/src/main.cpp b/microsc-psf/src/main.cpp index 967fdde..60387ea 100644 --- a/microsc-psf/src/main.cpp +++ b/microsc-psf/src/main.cpp @@ -4,6 +4,10 @@ // This must come after #include "make_psf.h" +#ifdef LSTSQ_USE_EIGEN +#include "linsolver.h" +#endif + namespace { // Return the range [0, N), excluding N. @@ -137,7 +141,7 @@ makePSF(microscope_params_t params, scale_t voxel, scale_t vol { // Define the basis of Bessel functions. // Shape: number of basis function by number of rho samples. - auto&& J = besselj<0>(scaling_factor * Rho); +#define J besselj<0>(scaling_factor * Rho) // Compute the approximation to the sampled pupil phase by finding the least squares // solution to the complex coefficients of the Fourier-Bessel expansion. @@ -146,10 +150,35 @@ makePSF(microscope_params_t params, scale_t voxel, scale_t vol // // Note: Armadillo does not have solver for real-valued matrix and complex-valued vector. // Reference: https://arma.sourceforge.net/armadillo_solver_2020.pdf - cx_mat&& Ci = (precision.solver == SandersonAndCurtin2020) - ? solve(conv_to::from(J.t()), phase.t()) - : cx_mat(pinv(J.t()) * phase.t()); - + cx_mat Ci; + switch (precision.solver) { +#ifdef LSTSQ_USE_LAPACK + case SandersonAndCurtin2020: + Ci = solve(conv_to::from(J.t()), phase.t()); + break; + case PenroseInverse: + Ci = pinv(J.t()) * phase.t(); + break; +#else + case SandersonAndCurtin2020: + case PenroseInverse: + throw std::invalid_argument(R"(LAPACK not found. +Re-configure the C++ project with +"meson configure -Darmadillo-code:lapack=..." and try again.)"); + break; +#endif + + case EigenBdcSVD: { +#ifndef LSTSQ_USE_EIGEN + throw std::invalid_argument(R"(Eigen::BdcSVD solver not implemented. +Re-configure the C++ project with +"meson configure -Duse_eigen=true" and try again.)"); +#else + constexpr bool always_transpose_phase = true; + Ci = microsc_psf::internal::solveWithEigen(J.t(), phase); +#endif + } + } const cx_mat ciEle = Ele * Ci; PSF0 = real(ciEle % conj(ciEle)); }