598 lines
16 KiB
C++
598 lines
16 KiB
C++
//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
|
|
|
|
/**
|
|
\file eigen_extensions.hpp
|
|
|
|
\brief Bertini extensions of the Eigen linear algebra library.
|
|
*/
|
|
|
|
#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
|
|
};
|
|
|
|
|
|
|
|
/**
|
|
\brief This templated struct permits us to use the mpfr_complex type in Eigen matrices.
|
|
|
|
Provides methods to get the epsilon, dummy_precision, lowest, highest functions, largely by inheritance from the NumTraits<mpfr_real> contained in mpfr_extensions.
|
|
*/
|
|
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>;
|
|
|
|
|
|
/**
|
|
\brief Checks whether the number of rows and columns are either 0
|
|
*/
|
|
template<typename Derived>
|
|
inline
|
|
bool IsEmpty(Eigen::MatrixBase<Derived> const & v)
|
|
{
|
|
return (v.rows()==0) || (v.cols()==0);
|
|
}
|
|
|
|
|
|
/**
|
|
\brief Get the precision of an Eigen object. If the object is empty, it's the precision of a default-constructed Scalar. If it actually has content, then it's the precision of the first element.
|
|
|
|
If you require that the object be of uniform precision when you check, use PrecisionRequireUniform
|
|
*/
|
|
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));
|
|
}
|
|
|
|
|
|
/**
|
|
\brief Query the precision of an Eigen object, requiring it to be non-empty, and of uniform precision
|
|
|
|
\throw runtime_error if not uniform precision
|
|
Dies in an Eigen assert if empty (assuming you didn't disable these)
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/**
|
|
\brief Set the precision of an Eigen object.
|
|
|
|
\param v the matrix to change
|
|
\param prec the new precision
|
|
*/
|
|
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);
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
Test a numbers being very small.
|
|
|
|
Compares the number against machine epsilon (or software epsilon if a multiple precision type), times 100.
|
|
|
|
\note Machine epsilon for doubles is about 1e-16, for mpfr_real, it's 10^-current precision.
|
|
|
|
\param testme The number to test
|
|
|
|
\return true, if the number is very small. False otherwise.
|
|
*/
|
|
template<typename T>
|
|
inline
|
|
bool IsSmallValue(T const& testme)
|
|
{
|
|
using std::abs;
|
|
return abs(testme) <= Eigen::NumTraits<T>::epsilon()*100;
|
|
}
|
|
|
|
/**
|
|
\brief Check whether two values are very close to each other.
|
|
|
|
\e The tolerance for being close
|
|
See \url http://www.boost.org/doc/libs/1_34_0/libs/test/doc/components/test_tools/floating_point_comparison.html, for example.
|
|
*/
|
|
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) ;
|
|
}
|
|
|
|
/**
|
|
Test two numbers for having large ratio.
|
|
|
|
The basis for comparison is Eigen's fuzzy precision, Eigen::NumTraits<T>::dummy_precision();
|
|
|
|
\note For doubles, the threshold is 1e-12, for mpfr_real is 1e3*current precision.
|
|
|
|
\param numerator The numerator for the ratio.
|
|
\param denomenator The denomenator for the ratio.
|
|
|
|
\return true, if the ratio is very large, false otherwise
|
|
*/
|
|
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
|
|
};
|
|
|
|
/**
|
|
\brief Check the diagonal elements of an LU decomposition for small values and large ratios.
|
|
|
|
\return Success if things are ok. LargeChange or SmallValue if one is found.
|
|
|
|
This function requires a square non-empty matrix.
|
|
|
|
\tparam Derived Matrix type from Eigen.
|
|
*/
|
|
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;
|
|
}
|
|
|
|
/**
|
|
\brief Make a Kahan matrix with a given number type.
|
|
*/
|
|
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;
|
|
}
|
|
|
|
|
|
/**
|
|
\brief Make a random matrix of units (numbers with norm 1).
|
|
|
|
\return The random matrix of units.
|
|
|
|
\param rows The number of rows.
|
|
\param cols The number of columns.
|
|
\tparam NumberType the type of number to fill the matrix with.
|
|
*/
|
|
template <typename NumberType>
|
|
inline
|
|
Mat<NumberType> RandomOfUnits(uint rows, uint cols)
|
|
{
|
|
return Mat<NumberType>(rows,cols).unaryExpr([](NumberType const& x) { return RandomUnit<NumberType>(); });
|
|
}
|
|
|
|
/**
|
|
\brief Make a random vector of units (numbers with norm 1).
|
|
|
|
\return The random vector of units.
|
|
|
|
\param size The length of the vector.
|
|
\tparam NumberType the type of number to fill the vector with.
|
|
*/
|
|
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
|