eigen_extensions.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//eigen_extensions.hpp is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//eigen_extensions.hpp 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 General Public License
//along with eigen_extensions.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) 2015 - 2021 by Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license,
// as well as COPYING.  Bertini2 is provided with permitted
// additional terms in the b2/licenses/ directory.
 
// individual authors of this file include:
// silviana amethyst, University of Wisconsin Eau Claire

 
#ifndef BERTINI_EIGEN_EXTENSIONS_HPP
#define BERTINI_EIGEN_EXTENSIONS_HPP
 
 
#include "bertini2/num_traits.hpp"
 
#include <boost/serialization/complex.hpp>
#include <boost/serialization/array.hpp>
#include <boost/serialization/split_member.hpp>
 
#define EIGEN_DENSEBASE_PLUGIN "bertini2/eigen_serialization_addon.hpp"
 
#include <Eigen/Core>
 
namespace {
using mpfr_real = bertini::mpfr_float;
using mpfr_complex = bertini::mpfr_complex;
}
 
namespace Eigen {
 
    template<> struct NumTraits<mpfr_real> : GenericNumTraits<mpfr_real> // permits to get the epsilon, dummy_precision, lowest, highest functions
    {
 
        typedef mpfr_real Real;
        typedef mpfr_real NonInteger;
        typedef mpfr_real Nested;
        enum {
            IsComplex = 0,
            IsInteger = 0,
            IsSigned = 1,
            RequireInitialization = 1, // yes, require initialization, otherwise get crashes
            ReadCost = 20,
            AddCost = 30,
            MulCost = 40
        };
 
 
        inline static Real highest() {
 
            return (mpfr_real(1) - epsilon()) * pow(mpfr_real(2),mpfr_get_emax()-1);//);//DefaultPrecision());
        }
 
        inline static Real lowest() {
            return -highest();
        }
 
        inline static Real dummy_precision()
        {
            using bertini::DefaultPrecision;
            return pow( mpfr_real(10),-int(DefaultPrecision()-3));
        }
 
        inline static Real epsilon()
        {
            using bertini::DefaultPrecision;
            return pow(mpfr_real(10),-int(DefaultPrecision()));
        }
 
        static inline int digits10()
        {
            return bertini::DefaultPrecision();
            // return internal::default_digits10_impl<T>::run();
        }
        //http://www.manpagez.com/info/mpfr/mpfr-2.3.2/mpfr_31.php
    };
 
 
    template<typename Expr1,typename Expr2, typename Expr3>
    struct NumTraits<
        boost::multiprecision::detail::expression<Expr1,
        Expr2, Expr3, void, void>
            > : NumTraits<mpfr_real> // permits to get the epsilon, dummy_precision, lowest, highest functions
    {
 
        typedef mpfr_real Real;
        typedef mpfr_real NonInteger;
        typedef mpfr_real Nested;
 
        //http://www.manpagez.com/info/mpfr/mpfr-2.3.2/mpfr_31.php
    };
 
 

    template<> struct NumTraits<mpfr_complex> : NumTraits<mpfr_real>
    {
        typedef mpfr_real Real;
        typedef mpfr_real NonInteger;
        typedef mpfr_complex Nested;// Nested;
        enum {
            IsComplex = 1,
            IsInteger = 0,
            IsSigned = 1,
            RequireInitialization = 1, // yes, require initialization, otherwise get crashes
            ReadCost = 2 * NumTraits<Real>::ReadCost,
            AddCost = 2 * NumTraits<Real>::AddCost,
            MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
        };
 
    };
 
    namespace internal {
        template<>
        struct abs2_impl<mpfr_complex>
        {
            static inline mpfr_real run(const mpfr_complex& x)
            {
                return real(x)*real(x) + imag(x)*imag(x);
            }
        };
 
 
        template<> inline mpfr_complex random<mpfr_complex>()
        {
            return bertini::multiprecision::rand();
        }
 
        template<> inline mpfr_complex random<mpfr_complex>(const mpfr_complex& a, const mpfr_complex& b)
        {
            return a + (b-a) * random<mpfr_complex>();
        }
 
        template<>
        struct conj_helper<mpfr_complex, mpfr_complex, false, true>
        {
            typedef mpfr_complex Scalar;
            EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
            { return c + pmul(x,y); }
 
            EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
            { return Scalar(real(x)*real(y) + imag(x)*imag(y), imag(x)*real(y) - real(x)*imag(y)); }
        };
 
        template<>
        struct conj_helper<mpfr_complex, mpfr_complex, true, false>
        {
            typedef mpfr_complex Scalar;
            EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
            { return c + pmul(x,y); }
 
            EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
            { return Scalar(real(x)*real(y) + imag(x)*imag(y), real(x)*imag(y) - imag(x)*real(y)); }
        };
 
        // see https://forum.kde.org/viewtopic.php?f=74&t=111176
 
        //int
        template<>
        struct scalar_product_traits<int,mpfr_complex>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        template<>
        struct scalar_product_traits<mpfr_complex, int>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        //long
        template<>
        struct scalar_product_traits<long,mpfr_complex>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        template<>
        struct scalar_product_traits<mpfr_complex, long>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        //long long
        template<>
        struct scalar_product_traits<long long,mpfr_complex>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        template<>
        struct scalar_product_traits<mpfr_complex, long long>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
 
        //mpfr_real
        template<>
        struct scalar_product_traits<mpfr_real,mpfr_complex>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        template<>
        struct scalar_product_traits<mpfr_complex, mpfr_real>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        //mpz_int
        template<>
        struct scalar_product_traits<bertini::mpz_int,mpfr_complex>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
        template<>
        struct scalar_product_traits<mpfr_complex, bertini::mpz_int>
        {
            enum { Defined = 1 };
            typedef mpfr_complex ReturnType;
        };
 
    } // re: namespace internal
} // re: namespace Eigen
 
 
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Eigenvalues>
#include <Eigen/SVD>
#include <Eigen/Dense>
 
namespace bertini {
 
    template<typename NumType> using Vec = Eigen::Matrix<NumType, Eigen::Dynamic, 1>;
    template<typename NumType> using Mat = Eigen::Matrix<NumType, Eigen::Dynamic, Eigen::Dynamic>;
 

    template<typename Derived>
    inline
    bool IsEmpty(Eigen::MatrixBase<Derived> const & v)
    {
        return (v.rows()==0) || (v.cols()==0);
    }
 

    template<typename Derived>
    inline
    unsigned Precision(Eigen::MatrixBase<Derived> const & v)
    {
        if (IsEmpty(v))
            throw std::runtime_error("getting precision of empty object");
 
        return Precision(v(0,0));
    }
 

    template<typename Derived>
    unsigned PrecisionRequireUniform(Eigen::MatrixBase<Derived> const & v)
    {
        const auto a = Precision(v(0,0));
        for (int ii=0; ii<v.cols(); ++ii)
            for (int jj=0; jj<v.rows(); ++jj)
            {
                if (a!=Precision(v(jj,ii)))
                {
                    std::stringstream ss;
                    ss << "non-uniform precision in object! (" << a << "!=" << Precision(v(jj,ii)) <<   ") at position " << jj << "," << ii;
                    throw std::runtime_error(ss.str());
                }
            }
        return a;
    }
 

    template<typename Derived>
    void Precision(Eigen::MatrixBase<Derived> & v, unsigned prec)
    {
        using bertini::Precision;
        for (int ii=0; ii<v.cols(); ++ii)
            for (int jj=0; jj<v.rows(); ++jj)
                Precision(v(jj,ii),prec);
    }
 
 

    template<typename T>
    inline
    bool IsSmallValue(T const& testme)
    {
        using std::abs;
        return abs(testme) <= Eigen::NumTraits<T>::epsilon()*100;
    }

    template<typename T>
    inline
    bool IsSymmRelDiffSmall(T const& a, T const& b, typename Eigen::NumTraits<T>::Real const& e)
    {
        using std::abs;
        if (a==b)
            return true;
 
        typename Eigen::NumTraits<T>::Real c = abs(a-b);
        return (c/abs(a) <= e) || (c/abs(b) <= e) ;
    }

    template<typename T>
    inline
    bool IsLargeChange(T const& numerator, T const& denomenator)
    {
        static_assert(!Eigen::NumTraits<T>::IsInteger, "IsLargeChange cannot be used safely on non-integral types");
        using std::abs;
        return abs(numerator/denomenator) >= 1/Eigen::NumTraits<T>::dummy_precision();
    }
 
    enum class MatrixSuccessCode
    {
        Success,
        LargeChange,
        SmallValue
    };

    template <typename Derived>
    MatrixSuccessCode LUPartialPivotDecompositionSuccessful(Eigen::MatrixBase<Derived> const& LU)
    {
        #ifndef BERTINI_DISABLE_ASSERTS
            assert(LU.rows()==LU.cols() && "non-square matrix in LUPartialPivotDecompositionSuccessful");
            assert(LU.rows()>0 && "empty matrix in LUPartialPivotDecompositionSuccessful");
        #endif
 
            // this loop won't test entry (0,0).  it's tested separately after.
        for (unsigned int ii = LU.rows()-1; ii > 0; ii--)
        {
            if (IsSmallValue(LU(ii,ii)))
            {
                return MatrixSuccessCode::SmallValue;
            }
 
            if (IsLargeChange(LU(ii-1,ii-1),LU(ii,ii)))
            {
                return MatrixSuccessCode::LargeChange;
            }
        }
 
        // this line is the reason for the above assert on non-empty matrix.
        if (IsSmallValue(LU(0,0)))
        {
            return MatrixSuccessCode::SmallValue;
        }
 
        return MatrixSuccessCode::Success;
    }

    template <typename NumberType>
    Eigen::Matrix<NumberType, Eigen::Dynamic, Eigen::Dynamic> KahanMatrix(unsigned int mat_size, NumberType c)
    {
        using std::sqrt;
        NumberType s, scale(1.0);
        s = sqrt( (NumberType(1.0)-c) * (NumberType(1.0)+c) );
 
 
        Eigen::Matrix<NumberType, Eigen::Dynamic, Eigen::Dynamic> A(mat_size,mat_size);
 
 
        for (unsigned int ii=0; ii<mat_size; ii++) {
            for (unsigned int jj=0; jj<ii; jj++) {
                A(ii,jj) = NumberType(0.0);
            }
            for (unsigned int jj=ii; jj<mat_size; jj++) {
                A(ii,jj) = -c * NumberType(1.0);
            }
        }
 
 
        for (unsigned int ii=0; ii<mat_size; ii++) {
            A(ii,ii) += NumberType(1)+c;
        }
 
        for (unsigned int jj=0; jj<mat_size; jj++){
            for (unsigned int kk=0;kk<mat_size;kk++){
                A(jj,kk) *= scale;
            }
            scale *= s;
        }
 
        for (unsigned int jj=0;jj<mat_size;jj++){
            for (unsigned int kk=0;kk<mat_size;kk++){
                A(kk,jj)/= NumberType(jj) + NumberType(1);
            }
        }
        return A;
    }
 

    template <typename NumberType>
    inline
    Mat<NumberType> RandomOfUnits(uint rows, uint cols)
    {
        return Mat<NumberType>(rows,cols).unaryExpr([](NumberType const& x) { return RandomUnit<NumberType>(); });
    }

    template <typename NumberType>
    inline
    Vec<NumberType> RandomOfUnits(uint size)
    {
        return Vec<NumberType>(size).unaryExpr([](NumberType const& x) { return RandomUnit<NumberType>(); });
    }
 
}
 
 
// // the following code comes from
// // https://stackoverflow.com/questions/18382457/eigen-and-boostserialize
// // and adds support for serialization of Eigen types
// //
// // question asked by user Gabriel and answered by Shmuel Levine
// // answer code repeated here verbatim.
// // please update this comment if this code is changed,
// // and post the modifications to the above referenced post on SO.
// namespace boost{
//     namespace serialization{
 
//         template<   class Archive,
//                     class S,
//                     int Rows_,
//                     int Cols_,
//                     int Ops_,
//                     int MaxRows_,
//                     int MaxCols_>
//         inline void save(
//             Archive & ar,
//             const Eigen::Matrix<S, Rows_, Cols_, Ops_, MaxRows_, MaxCols_> & g,
//             const unsigned int version)
//             {
//                 int rows = g.rows();
//                 int cols = g.cols();
 
//                 ar & rows;
//                 ar & cols;
//                 ar & boost::serialization::make_array(g.data(), rows * cols);
//             }
 
//         template<   class Archive,
//                     class S,
//                     int Rows_,
//                     int Cols_,
//                     int Ops_,
//                     int MaxRows_,
//                     int MaxCols_>
//         inline void load(
//             Archive & ar,
//             Eigen::Matrix<S, Rows_, Cols_, Ops_, MaxRows_, MaxCols_> & g,
//             const unsigned int version)
//         {
//             int rows, cols;
//             ar & rows;
//             ar & cols;
//             g.resize(rows, cols);
//             ar & boost::serialization::make_array(g.data(), rows * cols);
//         }
 
//         template<   class Archive,
//                     class S,
//                     int Rows_,
//                     int Cols_,
//                     int Ops_,
//                     int MaxRows_,
//                     int MaxCols_>
//         inline void serialize(
//             Archive & ar,
//             Eigen::Matrix<S, Rows_, Cols_, Ops_, MaxRows_, MaxCols_> & g,
//             const unsigned int version)
//         {
//             split_free(ar, g, version);
//         }
 
 
//     } // namespace serialization
// } // namespace boost
 
 
 
 
#endif