Homec4science
Diffusion µSpectre (master) Browse

muSpectre/src/libmugrid/eigen_tools.hh
fa590e5a9f93master

Last Change
Raw File

eigen_tools.hh
View Options

/**
* @file eigen_tools.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Sep 2017
*
* @brief small tools to be used with Eigen
*
* Copyright © 2017 Till Junge
*
* µGrid is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µGrid 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 Lesser General Public License
* along with µGrid; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* * Boston, MA 02111-1307, USA.
*
* Additional permission under GNU GPL version 3 section 7
*
* If you modify this Program, or any covered work, by linking or combining it
* with proprietary FFT implementations or numerical libraries, containing parts
* covered by the terms of those libraries' licenses, the licensors of this
* Program grant you additional permission to convey the resulting work.
*/
#ifndef SRC_LIBMUGRID_EIGEN_TOOLS_HH_
#define SRC_LIBMUGRID_EIGEN_TOOLS_HH_
#include "grid_common.hh"
#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>
#include <type_traits>
#include <utility>
namespace muGrid {
/* ---------------------------------------------------------------------- */
namespace internal {
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t order, Dim_t dim, Dim_t... dims>
struct SizesByOrderHelper {
//! type to use
using Sizes =
typename SizesByOrderHelper<order - 1, dim, dim, dims...>::Sizes;
};
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t dim, Dim_t... dims>
struct SizesByOrderHelper<0, dim, dims...> {
//! type to use
using Sizes = Eigen::Sizes<dims...>;
};
} // namespace internal
//! Creates a Eigen::Sizes type for a Tensor defined by an order and dim
template <Dim_t order, Dim_t dim>
struct SizesByOrder {
static_assert(order > 0, "works only for order greater than zero");
//! `Eigen::Sizes`
using Sizes =
typename internal::SizesByOrderHelper<order - 1, dim, dim>::Sizes;
};
/* ---------------------------------------------------------------------- */
namespace internal {
/* ---------------------------------------------------------------------- */
//! Call a passed lambda with the unpacked sizes as arguments
template <Dim_t order, typename Fun_t, Dim_t dim, Dim_t... args>
struct CallSizesHelper {
//! applies the call
static decltype(auto) call(Fun_t && fun) {
static_assert(order > 0, "can't handle empty sizes b)");
return CallSizesHelper<order - 1, Fun_t, dim, dim, args...>::call(fun);
}
};
/* ---------------------------------------------------------------------- */
template <typename Fun_t, Dim_t dim, Dim_t... args>
//! Call a passed lambda with the unpacked sizes as arguments
struct CallSizesHelper<0, Fun_t, dim, args...> {
//! applies the call
static decltype(auto) call(Fun_t && fun) { return fun(args...); }
};
} // namespace internal
/**
* takes a lambda and calls it with the proper `Eigen::Sizes`
* unpacked as arguments. Is used to call constructors of a
* `Eigen::Tensor` or map thereof in a context where the spatial
* dimension is templated
*/
template <Dim_t order, Dim_t dim, typename Fun_t>
inline decltype(auto) call_sizes(Fun_t && fun) {
static_assert(order > 1, "can't handle empty sizes");
return internal::CallSizesHelper<order - 1, Fun_t, dim, dim>::call(
std::forward<Fun_t>(fun));
}
// compile-time square root
static constexpr Dim_t ct_sqrt(Dim_t res, Dim_t l, Dim_t r) {
if (l == r) {
return r;
} else {
const auto mid = (r + l) / 2;
if (mid * mid >= res) {
return ct_sqrt(res, l, mid);
} else {
return ct_sqrt(res, mid + 1, r);
}
}
}
static constexpr Dim_t ct_sqrt(Dim_t res) { return ct_sqrt(res, 1, res); }
namespace EigenCheck {
/**
* Structure to determine whether an expression can be evaluated
* into a `Eigen::Matrix`, `Eigen::Array`, etc. and which helps
* determine compile-time size
*/
template <class Derived>
struct is_matrix {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! evaluated test
constexpr static bool value{
std::is_same<typename Eigen::internal::traits<T>::XprKind,
Eigen::MatrixXpr>::value};
};
/**
* Helper class to check whether an `Eigen::Array` or
* `Eigen::Matrix` is statically sized
*/
template <class Derived>
struct is_fixed {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! evaluated test
constexpr static bool value{T::SizeAtCompileTime != Eigen::Dynamic};
};
/**
* Helper class to check whether an `Eigen::Array` or `Eigen::Matrix` is a
* static-size and square.
*/
template <class Derived>
struct is_square {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
//! true if the object is square and statically sized
constexpr static bool value{
(T::RowsAtCompileTime == T::ColsAtCompileTime) && is_fixed<T>::value};
};
/**
* computes the dimension from a second order tensor represented
* square matrix or array
*/
template <class Derived>
struct tensor_dim {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
static_assert(is_matrix<T>::value,
"The type of t is not understood as an Eigen::Matrix");
static_assert(is_square<T>::value, "t's matrix isn't square");
//! evaluated dimension
constexpr static Dim_t value{T::RowsAtCompileTime};
};
//! computes the dimension from a fourth order tensor represented
//! by a square matrix
template <class Derived>
struct tensor_4_dim {
//! raw type for testing
using T = std::remove_reference_t<Derived>;
static_assert(is_matrix<T>::value,
"The type of t is not understood as an Eigen::Matrix");
static_assert(is_square<T>::value, "t's matrix isn't square");
//! evaluated dimension
constexpr static Dim_t value{ct_sqrt(T::RowsAtCompileTime)};
static_assert(value * value == T::RowsAtCompileTime,
"This is not a fourth-order tensor mapped on a square "
"matrix");
};
namespace internal {
template <Dim_t dim, Dim_t nb_row, Dim_t nb_col>
constexpr inline Dim_t get_rank() {
constexpr bool is_vec{(nb_row == dim) and (nb_col == 1)};
constexpr bool is_mat{(nb_row == dim) and (nb_col == nb_row)};
constexpr bool is_ten{(nb_row == dim * dim) and (nb_col == dim * dim)};
static_assert(is_vec or is_mat or is_ten,
"can't understand the data type as a first-, second-, or "
"fourth-order tensor");
if (is_vec) {
return firstOrder;
} else if (is_mat) {
return secondOrder;
} else if (is_ten) {
return fourthOrder;
}
}
} // namespace internal
/**
* computes the rank of a tensor given the spatial dimension
*/
template <class Derived, Dim_t Dim>
struct tensor_rank {
using T = std::remove_reference_t<Derived>;
static_assert(is_matrix<T>::value,
"The type of t is not understood as an Eigen::Matrix");
static constexpr Dim_t value{internal::get_rank<Dim, T::RowsAtCompileTime,
T::ColsAtCompileTime>()};
};
} // namespace EigenCheck
namespace log_comp {
//! Matrix type used for logarithm evaluation
template <Dim_t dim>
using Mat_t = Eigen::Matrix<Real, dim, dim>;
//! Vector type used for logarithm evaluation
template <Dim_t dim>
using Vec_t = Eigen::Matrix<Real, dim, 1>;
//! This is a static implementation of the explicit determination
//! of log(Tensor) following Jog, C.S. J Elasticity (2008) 93:
//! 141. https://doi.org/10.1007/s10659-008-9169-x
/* ---------------------------------------------------------------------- */
template <Dim_t dim, Dim_t i, Dim_t j = dim - 1>
struct Proj {
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
static_assert(dim > 0, "only works for positive dimensions");
return 1. / (eigs(i) - eigs(j)) *
(T - eigs(j) * Mat_t<dim>::Identity()) *
Proj<dim, i, j - 1>::compute(eigs, T);
}
};
//! catch the case when there's nothing to do
template <Dim_t dim, Dim_t other>
struct Proj<dim, other, other> {
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
static_assert(dim > 0, "only works for positive dimensions");
return Proj<dim, other, other - 1>::compute(eigs, T);
}
};
//! catch the normal tail case
template <Dim_t dim, Dim_t i>
struct Proj<dim, i, 0> {
static constexpr Dim_t j{0}; //!< short-hand
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
static_assert(dim > 0, "only works for positive dimensions");
return 1. / (eigs(i) - eigs(j)) *
(T - eigs(j) * Mat_t<dim>::Identity());
}
};
//! catch the tail case when the last dimension is i
template <Dim_t dim>
struct Proj<dim, 0, 1> {
static constexpr Dim_t i{0}; //!< short-hand
static constexpr Dim_t j{1}; //!< short-hand
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
static_assert(dim > 0, "only works for positive dimensions");
return 1. / (eigs(i) - eigs(j)) *
(T - eigs(j) * Mat_t<dim>::Identity());
}
};
//! catch the general tail case
template <>
struct Proj<1, 0, 0> {
static constexpr Dim_t dim{1}; //!< short-hand
static constexpr Dim_t i{0}; //!< short-hand
static constexpr Dim_t j{0}; //!< short-hand
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & /*eigs*/,
const Mat_t<dim> & /*T*/) {
return Mat_t<dim>::Identity();
}
};
//! Product term
template <Dim_t dim, Dim_t i>
inline decltype(auto) P(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
return Proj<dim, i>::compute(eigs, T);
}
//! sum term
template <Dim_t dim, Dim_t i = dim - 1>
struct Summand {
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
return std::log(eigs(i)) * P<dim, i>(eigs, T) +
Summand<dim, i - 1>::compute(eigs, T);
}
};
//! sum term
template <Dim_t dim>
struct Summand<dim, 0> {
static constexpr Dim_t i{0}; //!< short-hand
//! wrapped function (raison d'être)
static inline decltype(auto) compute(const Vec_t<dim> & eigs,
const Mat_t<dim> & T) {
return std::log(eigs(i)) * P<dim, i>(eigs, T);
}
};
//! sum implementation
template <Dim_t dim>
inline decltype(auto) Sum(const Vec_t<dim> & eigs, const Mat_t<dim> & T) {
return Summand<dim>::compute(eigs, T);
}
} // namespace log_comp
/**
* computes the matrix logarithm efficiently for dim=1, 2, or 3 for
* a diagonizable tensor. For larger tensors, better use the direct
* eigenvalue/vector computation
*/
template <Dim_t dim>
inline decltype(auto) logm(const log_comp::Mat_t<dim> & mat) {
using Mat = log_comp::Mat_t<dim>;
Eigen::SelfAdjointEigenSolver<Mat> Solver{};
Solver.computeDirect(mat, Eigen::EigenvaluesOnly);
return Mat{log_comp::Sum(Solver.eigenvalues(), mat)};
}
/**
* compute the spectral decomposition
*/
template <class Derived>
inline decltype(auto)
spectral_decomposition(const Eigen::MatrixBase<Derived> & mat) {
static_assert(Derived::SizeAtCompileTime != Eigen::Dynamic,
"works only for static matrices");
static_assert(Derived::RowsAtCompileTime == Derived::ColsAtCompileTime,
"works only for square matrices");
constexpr Dim_t dim{Derived::RowsAtCompileTime};
using Mat = log_comp::Mat_t<dim>;
Eigen::SelfAdjointEigenSolver<Mat> Solver{};
Solver.computeDirect(mat, Eigen::ComputeEigenvectors);
return Solver;
}
/**
* compute the matrix log. This may not be the most
* efficient way to do this
*/
template <Dim_t Dim>
using Decomp_t = Eigen::SelfAdjointEigenSolver<Eigen::Matrix<Real, Dim, Dim>>;
template <Dim_t Dim>
inline decltype(auto) logm_alt(const Decomp_t<Dim> & spectral_decomp) {
using Mat = log_comp::Mat_t<Dim>;
Mat retval{Mat::Zero()};
for (Dim_t i = 0; i < Dim; ++i) {
const Real & val = spectral_decomp.eigenvalues()(i);
auto & vec = spectral_decomp.eigenvectors().col(i);
retval += std::log(val) * vec * vec.transpose();
}
return retval;
}
template <class Derived>
inline decltype(auto) logm_alt(const Eigen::MatrixBase<Derived> & mat) {
static_assert(Derived::SizeAtCompileTime != Eigen::Dynamic,
"works only for static matrices");
static_assert(Derived::RowsAtCompileTime == Derived::ColsAtCompileTime,
"works only for square matrices");
constexpr Dim_t dim{Derived::RowsAtCompileTime};
using Mat = log_comp::Mat_t<dim>;
using Decomp_t = Eigen::SelfAdjointEigenSolver<Mat>;
Decomp_t decomp{spectral_decomposition(mat)};
return logm_alt(decomp);
}
/**
* compute the matrix exponential. This may not be the most
* efficient way to do this
*/
template <Dim_t Dim>
inline decltype(auto) expm(const Decomp_t<Dim> & spectral_decomp) {
using Mat = log_comp::Mat_t<Dim>;
Mat retval{Mat::Zero()};
for (Dim_t i = 0; i < Dim; ++i) {
const Real & val = spectral_decomp.eigenvalues()(i);
auto & vec = spectral_decomp.eigenvectors().col(i);
retval += std::exp(val) * vec * vec.transpose();
}
return retval;
}
template <class Derived>
inline decltype(auto) expm(const Eigen::MatrixBase<Derived> & mat) {
static_assert(Derived::SizeAtCompileTime != Eigen::Dynamic,
"works only for static matrices");
static_assert(Derived::RowsAtCompileTime == Derived::ColsAtCompileTime,
"works only for square matrices");
constexpr Dim_t Dim{Derived::RowsAtCompileTime};
using Mat = log_comp::Mat_t<Dim>;
using Decomp_t = Eigen::SelfAdjointEigenSolver<Mat>;
Decomp_t decomp{spectral_decomposition(mat)};
return expm(decomp);
}
} // namespace muGrid
#endif // SRC_LIBMUGRID_EIGEN_TOOLS_HH_
c4science · Help