Commit 66b0a5ef authored by Philip Trettner's avatar Philip Trettner

Merge branch 'feature/eigenvaluedecomposition' into 'develop'

Feature/eigenvaluedecomposition

See merge request !85
parents e18667ab 13203d0e
......@@ -3,6 +3,7 @@
#include <typed-geometry/feature/vector.hh>
#include <typed-geometry/functions/matrix/adjugate.hh>
#include <typed-geometry/functions/matrix/covariance.hh>
#include <typed-geometry/functions/matrix/determinant.hh>
#include <typed-geometry/functions/matrix/diag.hh>
#include <typed-geometry/functions/matrix/eigenvalues.hh>
......
......@@ -36,8 +36,9 @@ template <class T = void, class RangeT, class TransformF, class ReduceF>
using R = same_or<T, element_type<RangeT>>;
auto it = tg::begin(values);
using U = std::decay_t<decltype(f(t(R(*it)), t(R(*it))))>;
auto const e = tg::end(values);
auto r = t(R(*it));
U r = t(R(*it));
it++;
while (it != e)
{
......@@ -126,7 +127,21 @@ template <class T = void, class RangeT = void, class TransformT = identity_fun>
[[nodiscard]] constexpr auto arithmetic_mean(RangeT const& values, TransformT&& transform = {})
{
auto const s = sum<T>(values, transform);
return s / static_cast<decltype(s)>(tg::container_size(values));
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wsign-conversion"
#endif
#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable : 4244) // possible loss of data
#endif
return s / tg::container_size(values);
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#ifdef _MSC_VER
#pragma warning(pop)
#endif
}
template <class T = void, class RangeT = void, class TransformT = identity_fun>
......
#pragma once
#include <typed-geometry/functions/basic/statistics.hh>
#include <typed-geometry/functions/matrix/outer_product.hh>
#include <typed-geometry/types/mat.hh>
namespace tg
{
template <class PosRangeT, class Transform = identity_fun>
auto covariance_matrix(PosRangeT&& r, Transform&& t = {}) -> decltype(self_outer_product(t(*tg::begin(r)) - t(*tg::begin(r))))
{
using MatT = decltype(self_outer_product(t(*tg::begin(r)) - t(*tg::begin(r))));
auto const avg = average(r, t);
MatT cov;
for (auto const& e : r)
cov += self_outer_product(t(e) - avg);
return cov;
}
}
This diff is collapsed.
#pragma once
#include <typed-geometry/types/comp.hh>
#include <typed-geometry/types/array.hh>
#include <typed-geometry/types/mat.hh>
#include <typed-geometry/types/vec.hh>
#include <typed-geometry/functions/basic/scalar_math.hh>
#include <typed-geometry/functions/matrix/transpose.hh>
namespace tg
{
template <class T>
[[nodiscard]] constexpr comp<2, T> eigenvalues(mat<2, 2, T> const& m)
template <int D, class ScalarT>
struct eigen_decomposition_result
{
auto mp_half = (m[0][0] + m[1][1]) / T(2);
auto q = m[0][0] * m[1][1] - m[1][0] * m[0][1];
auto d = sqrt(max(T(0), mp_half * mp_half - q));
return {mp_half + d, mp_half - d};
vec<D, ScalarT> eigenvector;
ScalarT eigenvalue;
};
namespace detail
{
array<eigen_decomposition_result<2, float>, 2> eigen_decomposition_impl(mat<2, 2, float> const& m);
array<eigen_decomposition_result<2, double>, 2> eigen_decomposition_impl(mat<2, 2, double> const& m);
array<eigen_decomposition_result<3, float>, 3> eigen_decomposition_impl(mat<3, 3, float> const& m);
array<eigen_decomposition_result<3, double>, 3> eigen_decomposition_impl(mat<3, 3, double> const& m);
array<eigen_decomposition_result<4, float>, 4> eigen_decomposition_impl(mat<4, 4, float> const& m);
array<eigen_decomposition_result<4, double>, 4> eigen_decomposition_impl(mat<4, 4, double> const& m);
array<float, 2> eigenvalues_impl(mat<2, 2, float> const& m);
array<double, 2> eigenvalues_impl(mat<2, 2, double> const& m);
array<float, 3> eigenvalues_impl(mat<3, 3, float> const& m);
array<double, 3> eigenvalues_impl(mat<3, 3, double> const& m);
array<float, 4> eigenvalues_impl(mat<4, 4, float> const& m);
array<double, 4> eigenvalues_impl(mat<4, 4, double> const& m);
array<vec<2, float>, 2> eigenvectors_impl(mat<2, 2, float> const& m);
array<vec<2, double>, 2> eigenvectors_impl(mat<2, 2, double> const& m);
array<vec<3, float>, 3> eigenvectors_impl(mat<3, 3, float> const& m);
array<vec<3, double>, 3> eigenvectors_impl(mat<3, 3, double> const& m);
array<vec<4, float>, 4> eigenvectors_impl(mat<4, 4, float> const& m);
array<vec<4, double>, 4> eigenvectors_impl(mat<4, 4, double> const& m);
}
template <class ScalarT, int D>
[[nodiscard]] array<eigen_decomposition_result<D, ScalarT>, D> eigen_decomposition_symmetric(mat<D, D, ScalarT> const& m)
{
if constexpr (D == 1)
{
return {m[0], m[0]};
}
else if constexpr (D == 2)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigen_decomposition_impl(m);
}
else if constexpr (D == 3)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigen_decomposition_impl(m);
}
else if constexpr (D == 4)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigen_decomposition_impl(m);
}
else
{
static_assert(tg::always_false<D>, "unsupported dimension");
return {};
}
}
template <class T>
[[nodiscard]] constexpr comp<2, T> singular_values(mat<2, 2, T> const& m)
template <class ScalarT, int D>
[[nodiscard]] array<ScalarT, D> eigenvalues_symmetric(mat<D, D, ScalarT> const& m)
{
auto [s0, s1] = eigenvalues(m * transpose(m));
return {sqrt(s0), sqrt(s1)};
if constexpr (D == 1)
{
return {m[0], m[0]};
}
else if constexpr (D == 2)
{
auto mp_half = (m[0][0] + m[1][1]) / ScalarT(2);
auto q = m[0][0] * m[1][1] - m[1][0] * m[0][1];
auto d = sqrt(max(ScalarT(0), mp_half * mp_half - q));
return {mp_half + d, mp_half - d};
}
else if constexpr (D == 3)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigenvalues_impl(m);
}
else if constexpr (D == 4)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigenvalues_impl(m);
}
else
{
static_assert(tg::always_false<D>, "unsupported dimension");
return {};
}
}
template <class ScalarT, int D>
[[nodiscard]] array<vec<D, ScalarT>, D> eigenvectors_symmetric(mat<D, D, ScalarT> const& m)
{
if constexpr (D == 1)
{
return {m[0], m[0]};
}
else if constexpr (D == 2)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigenvectors_impl(m);
}
else if constexpr (D == 3)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigenvectors_impl(m);
}
else if constexpr (D == 4)
{
static_assert(std::is_same_v<ScalarT, float> || std::is_same_v<ScalarT, double>, "currently only suports float or double");
return detail::eigenvectors_impl(m);
}
else
{
static_assert(tg::always_false<D>, "unsupported dimension");
return {};
}
}
template <class ScalarT, int D>
[[nodiscard]] array<ScalarT, D> singular_values_symmetric(mat<D, D, ScalarT> const& m)
{
if constexpr (D == 1)
{
return {m[0]};
}
else if constexpr (D == 2)
{
auto [s0, s1] = eigenvalues(m * transpose(m));
return {sqrt(s0), sqrt(s1)};
}
else if constexpr (D == 3)
{
static_assert(tg::always_false<D>, "not yet supported");
}
else if constexpr (D == 4)
{
static_assert(tg::always_false<D>, "not yet supported");
}
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment