eigen.h

Go to the documentation of this file.

//--------------------------------------------------------------------------------------------------
//
//  UnitConversion: A compile-time c++23 unit conversion library with no dependencies
//
//--------------------------------------------------------------------------------------------------
//
// The MIT License (MIT)
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
// BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
// DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//
//--------------------------------------------------------------------------------------------------
//
// Copyright (c) 2016 Nic Holthaus
//
//--------------------------------------------------------------------------------------------------
//
//
//--------------------------------------------------------------------------------------------------
 
#pragma once
 
#ifndef units_eigen_h_
#define units_eigen_h_
 
#include <units/core.h>
 
#if defined __has_include
#if __has_include(<Eigen/Core>)
#include <Eigen/Core>
 
#include <cmath>
 
namespace Eigen
{
    template<units::UnitType U>
    struct NumTraits<U> : NumTraits<typename units::traits::unit_traits<U>::underlying_type>
    {
        using T          = typename units::traits::unit_traits<U>::underlying_type; 
        using Real       = U;
        using NonInteger = U;
        using Nested     = U;
        using Literal    = T;
 
        enum
        {
            IsComplex             = 0,
            IsInteger             = NumTraits<T>::IsInteger,
            IsSigned              = NumTraits<T>::IsSigned,
            RequireInitialization = 1,
            ReadCost              = 1,
            AddCost               = 1,
            MulCost               = 1
        };
    };

    template<units::UnitType U, class X>
    struct ScalarBinaryOpTraits<U, X, internal::scalar_product_op<U, X>>
    {
        using ReturnType = U; 
    };

    template<class X, units::UnitType U>
    struct ScalarBinaryOpTraits<X, U, internal::scalar_product_op<X, U>>
    {
        using ReturnType = U; 
    };

    template<units::UnitType U, class X>
    struct ScalarBinaryOpTraits<U, X, internal::scalar_quotient_op<U, X>>
    {
        using ReturnType = U; 
    };
} // namespace Eigen
 
namespace units
{
    template<class DerivedA, class DerivedB>
    auto unit_dot(const Eigen::MatrixBase<DerivedA>& lhs, const Eigen::MatrixBase<DerivedB>& rhs)
    {
        using UnitA   = typename DerivedA::Scalar;
        using UnitB   = typename DerivedB::Scalar;
        using Product = decltype(std::declval<UnitA>() * std::declval<UnitB>());
 
        Product accumulator = Product(0);
        for (Eigen::Index i = 0; i < lhs.size(); ++i)
            accumulator += lhs(i) * rhs(i);
        return accumulator;
    }

    template<class Derived>
    auto unit_squared_norm(const Eigen::MatrixBase<Derived>& v)
    {
        return unit_dot(v, v);
    }

    template<class Derived>
    auto unit_norm(const Eigen::MatrixBase<Derived>& v)
    {
        return units::sqrt(unit_squared_norm(v));
    }

    template<class Derived>
    auto unit_normalized(const Eigen::MatrixBase<Derived>& v)
    {
        using Unit       = typename Derived::Scalar;
        using Underlying = detail::floating_point_promotion_t<typename traits::unit_traits<Unit>::underlying_type>;
 
        const auto magnitude = unit_norm(v);
 
        Eigen::Matrix<Underlying, Derived::RowsAtCompileTime, 1> direction(v.size());
        for (Eigen::Index i = 0; i < v.size(); ++i)
            direction(i) = (v(i) / magnitude).template to<Underlying>();
        return direction;
    }

    template<class DerivedA, class DerivedB>
    auto unit_cross(const Eigen::MatrixBase<DerivedA>& lhs, const Eigen::MatrixBase<DerivedB>& rhs)
    {
        using UnitA   = typename DerivedA::Scalar;
        using UnitB   = typename DerivedB::Scalar;
        using Product = decltype(std::declval<UnitA>() * std::declval<UnitB>());
 
        Eigen::Matrix<Product, 3, 1> result;
        result(0) = lhs(1) * rhs(2) - lhs(2) * rhs(1);
        result(1) = lhs(2) * rhs(0) - lhs(0) * rhs(2);
        result(2) = lhs(0) * rhs(1) - lhs(1) * rhs(0);
        return result;
    }

    template<class MatrixDerived, class VectorDerived>
    auto unit_transform(const Eigen::MatrixBase<MatrixDerived>& matrix, const Eigen::MatrixBase<VectorDerived>& vector)
    {
        using Unit       = typename VectorDerived::Scalar;
        using Underlying = typename traits::unit_traits<Unit>::underlying_type;
 
        // Operate on each component's own-scale value (raw), so a ratio-scaled dimensionless input transforms in the
        // scale the caller sees and the result reconstructs on that same scale. Every intermediate is sized to its
        // runtime length so a dynamically-sized transform does not index a zero-length matrix.
        Eigen::Matrix<Underlying, VectorDerived::RowsAtCompileTime, 1> raw(vector.size());
        for (Eigen::Index i = 0; i < vector.size(); ++i)
            raw(i) = vector(i).raw();
 
        const Eigen::Matrix<Underlying, MatrixDerived::RowsAtCompileTime, 1> product = matrix * raw;
 
        Eigen::Matrix<Unit, MatrixDerived::RowsAtCompileTime, 1> result(product.size());
        for (Eigen::Index i = 0; i < product.size(); ++i)
            result(i) = Unit(product(i));
        return result;
    }
} // namespace units
 
#endif // __has_include(<Eigen/Core>)
#endif // defined __has_include
 
#endif // units_eigen_h_