Investigating multi precision libraries

[jenshnielsen]

Just some notes about using multi precision.

Eigen has unsupported MPFR C++ support http://eigen.tuxfamily.org/dox/unsupported/group__MPRealSupport__Module.html

This using MPFR C++ as a c++ wrapper to
MPFR multi precision floating point library.
Which in turn depends on a multi precision integer library. Normally
GMP However, on windows we likely need to use MPIR
see https://github.com/wbhart/mpir/tree/master/build.vc14 for some visual studio support.

[jenshnielsen]

A good starting point is probably to get the Eigen example running on windows.

Works on my Mac with:

clang++ test.cpp -I/usr/local/include/eigen3/ -I. -lmpfr -lgmp

[jenshnielsen]

It strikes me that we should be able to improve the precision in round overflow protection by changing it from dividing by a power of 10 (currently 10^50) to dividing by a power of 2 which can be done exactly for floating point numbers.

The conversion back when doing log10 can naturally not be done exactly then but that is only one operation and can be done as

math.log10(a/(2**b)) + b * math.log10(2)

[jenshnielsen]

Possibly relevant from CGAL http://www.cgal.org/exact.html thou somewhat lacking in any real details

[jenshnielsen]

An example of something that diverges from fitGlyR4ke.py

With Gcc

iteration # 7169; log-lik = 233149.204258
[  4.39168065e+05   4.55322454e+06   2.75936209e+04   6.19107695e+04
   8.07015209e+03   6.30184502e+04   3.58607958e+04   1.18087961e+03
   4.00347448e+03   5.76489604e+03   1.77272876e+04   4.69242928e+02
   1.53484232e+04   6.08467481e+03   1.00839783e+01   5.30306673e+02
   3.34845112e+02   1.18528813e+04   7.83808560e+03   2.54316451e+07
   2.74799567e+03   5.65415757e+05   1.19073627e-01   3.33573106e-02
   6.67684178e+03   5.50803653e+04   6.66499465e+02]
iteration # 7170; log-lik = 233149.204258
[  4.39168065e+05   4.55322454e+06   2.75936209e+04   6.19107695e+04
   8.07015209e+03   6.30184502e+04   3.58607958e+04   1.18087961e+03
   4.00347448e+03   5.76489604e+03   1.77272876e+04   4.69242928e+02
   1.53484232e+04   6.08467481e+03   1.00839783e+01   5.30306673e+02
   3.34845112e+02   1.18528813e+04   7.83808560e+03   2.54316451e+07
   2.74799567e+03   5.65415757e+05   1.19073627e-01   3.33573106e-02
   6.67684178e+03   5.50803653e+04   6.66499465e+02]
iteration # 7171; log-lik = 233149.204792
[  4.40493424e+05   4.54128124e+06   2.75978784e+04   6.18953993e+04
   8.07929151e+03   6.30566727e+04   3.58492216e+04   1.18239595e+03
   4.00785231e+03   5.76743417e+03   1.77141563e+04   4.67914563e+02
   1.52776836e+04   6.08899781e+03   1.01180468e+01   5.30782378e+02
   3.34399842e+02   1.18454027e+04   7.84042543e+03   2.53133933e+07
   2.75072352e+03   5.63404190e+05   1.19096829e-01   3.33501525e-02
   6.67720048e+03   5.50575109e+04   6.66274946e+02]
iteration # 7172; log-lik = 233149.204792
[  4.40493424e+05   4.54128124e+06   2.75978784e+04   6.18953993e+04
   8.07929151e+03   6.30566727e+04   3.58492216e+04   1.18239595e+03
   4.00785231e+03   5.76743417e+03   1.77141563e+04   4.67914563e+02
   1.52776836e+04   6.08899781e+03   1.01180468e+01   5.30782378e+02
   3.34399842e+02   1.18454027e+04   7.84042543e+03   2.53133933e+07
   2.75072352e+03   5.63404190e+05   1.19096829e-01   3.33501525e-02
   6.67720048e+03   5.50575109e+04   6.66274946e+02]

With Clang

iteration # 7169; log-lik = 233149.204258
[  4.39168065e+05   4.55322454e+06   2.75936209e+04   6.19107695e+04
   8.07015209e+03   6.30184502e+04   3.58607958e+04   1.18087961e+03
   4.00347448e+03   5.76489604e+03   1.77272876e+04   4.69242928e+02
   1.53484232e+04   6.08467481e+03   1.00839783e+01   5.30306673e+02
   3.34845112e+02   1.18528813e+04   7.83808560e+03   2.54316451e+07
   2.74799567e+03   5.65415757e+05   1.19073627e-01   3.33573106e-02
   6.67684178e+03   5.50803653e+04   6.66499465e+02]
iteration # 7170; log-lik = 233149.204258
[  4.39168065e+05   4.55322454e+06   2.75936209e+04   6.19107695e+04
   8.07015209e+03   6.30184502e+04   3.58607958e+04   1.18087961e+03
   4.00347448e+03   5.76489604e+03   1.77272876e+04   4.69242928e+02
   1.53484232e+04   6.08467481e+03   1.00839783e+01   5.30306673e+02
   3.34845112e+02   1.18528813e+04   7.83808560e+03   2.54316451e+07
   2.74799567e+03   5.65415757e+05   1.19073627e-01   3.33573106e-02
   6.67684178e+03   5.50803653e+04   6.66499465e+02]
iteration # 7171; log-lik = 233149.204636
[  4.41052761e+05   4.55003103e+06   2.75925932e+04   6.18974291e+04
   8.07811319e+03   6.30488833e+04   3.58334764e+04   1.18077263e+03
   4.00829165e+03   5.76705755e+03   1.77202964e+04   4.67890671e+02
   1.52908901e+04   6.09078592e+03   1.01109456e+01   5.30755317e+02
   3.34429074e+02   1.18485003e+04   7.84223692e+03   2.53147142e+07
   2.75048901e+03   5.63541642e+05   1.19088442e-01   3.33508712e-02
   6.67669155e+03   5.50741422e+04   6.66502946e+02]
iteration # 7172; log-lik = 233149.204636
[  4.41052761e+05   4.55003103e+06   2.75925932e+04   6.18974291e+04
   8.07811319e+03   6.30488833e+04   3.58334764e+04   1.18077263e+03
   4.00829165e+03   5.76705755e+03   1.77202964e+04   4.67890671e+02
   1.52908901e+04   6.09078592e+03   1.01109456e+01   5.30755317e+02
   3.34429074e+02   1.18485003e+04   7.84223692e+03   2.53147142e+07
   2.75048901e+03   5.63541642e+05   1.19088442e-01   3.33508712e-02
   6.67669155e+03   5.50741422e+04   6.66502946e+02]

[jenshnielsen]

Trying to bulk replace the real data type with mpreal fails first because mpfr ++ is not constexpr safe.

In file included from /Users/jhn/ucl/dcprogs/HJCFIT/likelihood/determinant_equation.h:26:0,
                 from /Users/jhn/ucl/dcprogs/HJCFIT/likelihood/determinant_equation.cc:25:
/Users/jhn/ucl/dcprogs/HJCFIT/likelihood/laplace_survivor.h:102:40: error: the type 'const t_real {aka const mpfr::mpreal}' of constexpr variable 'DCProgs::LaplaceSurvivor::ZERO' is not literal
         constexpr static t_real ZERO = 1e-12;
                                        ^~~~~
In file included from /usr/local/include/eigen3/unsupported/Eigen/MPRealSupport:16:0,
                 from /Users/jhn/ucl/dcprogs/HJCFIT/build/DCProgsConfig.h:32,
                 from /Users/jhn/ucl/dcprogs/HJCFIT/likelihood/determinant_equation.cc:21:
/Users/jhn/ucl/dcprogs/HJCFIT/mpreal.h:145:7: note: 'mpfr::mpreal' is not literal because:
 class mpreal {
       ^~~~~~
/Users/jhn/ucl/dcprogs/HJCFIT/mpreal.h:145:7: note:   'mpfr::mpreal' has a non-trivial destructor
In file included from /Users/jhn/ucl/dcprogs/HJCFIT/likelihood/determinant_equation.h:26:0,
                 from /Users/jhn/ucl/dcprogs/HJCFIT/likelihood/determinant_equation.cc:25:
/Users/jhn/ucl/dcprogs/HJCFIT/likelihood/laplace_survivor.h:102:40: error: in-class initialization of static data member 'const t_real DCProgs::LaplaceSurvivor::ZERO' of non-literal type
         constexpr static t_real ZERO = 1e-12;
                                        ^~~~~
/Users/jhn/ucl/dcprogs/HJCFIT/likelihood/laplace_survivor.h:102:40: error: call to non-constexpr function 'mpfr::mpreal::mpreal(double, mpfr_prec_t, mpfr_rnd_t)'
make[2]: *** [likelihood/CMakeFiles/likelihood.dir/determinant_equation.cc.o] Error 1

[jenshnielsen]

Hacking HJCFIT to disable constexpr gets the compilation a bit further:

In file included from /usr/local/include/eigen3/unsupported/Eigen/MatrixFunctions:58:
/usr/local/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:186:3: error: no matching member function for call to 'computeUV'
  computeUV(RealScalar());
  ^~~~~~~~~
/usr/local/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:422:10: note: in instantiation of function template specialization 'Eigen::MatrixExponential<Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1> >::compute<Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1> >' requested here
      me.compute(result);
         ^
/usr/local/include/eigen3/Eigen/src/Core/ReturnByValue.h:61:42: note: in instantiation of function template specialization 'Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > >::evalTo<Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1> >' requested here
    { static_cast<const Derived*>(this)->evalTo(dst); }
                                         ^
/usr/local/include/eigen3/Eigen/src/Core/Matrix.h:311:13: note: in instantiation of function template specialization 'Eigen::ReturnByValue<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > >::evalTo<Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1> >' requested here
      other.evalTo(*this);
            ^
/usr/local/include/eigen3/Eigen/src/Core/ProductBase.h:98:9: note: in instantiation of function template specialization 'Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>::Matrix<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > >' requested here
      : m_lhs(a_lhs), m_rhs(a_rhs)
        ^
/usr/local/include/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:393:54: note: in instantiation of member function 'Eigen::ProductBase<Eigen::GeneralProduct<Eigen::ReturnByValue<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > >, Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false>, 5>, Eigen::ReturnByValue<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > >, Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> >::ProductBase' requested here
    GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
                                                     ^
/usr/local/include/eigen3/Eigen/src/Core/GeneralProduct.h:598:10: note: in instantiation of member function 'Eigen::GeneralProduct<Eigen::ReturnByValue<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > >, Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false>, 5>::GeneralProduct' requested here
  return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
         ^
/Users/jhn/ucl/dcprogs/HJCFIT/likelihood/idealG.h:76:41: note: in instantiation of function template specialization 'Eigen::MatrixBase<Eigen::ReturnByValue<Eigen::MatrixExponentialReturnValue<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<mpfr::mpreal>, const Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> > > > >::operator*<Eigen::Block<const Eigen::Matrix<mpfr::mpreal, -1, -1, 0, -1, -1>, -1, -1, false> >' requested here
        { return (t*QMatrix::ff()).exp()*QMatrix::fa(); }
                                        ^
/usr/local/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:119:10: note: candidate function not viable: no known conversion from 'RealScalar' (aka 'mpfr::mpreal') to 'double' for 1st argument
    void computeUV(double);
         ^
/usr/local/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:125:10: note: candidate function not viable: no known conversion from 'RealScalar' (aka 'mpfr::mpreal') to 'float' for 1st argument
    void computeUV(float);
         ^
/usr/local/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:131:10: note: candidate function not viable: no known conversion from 'RealScalar' (aka 'mpfr::mpreal') to 'long double' for 1st argument
    void computeUV(long double);
         ^
1 error generated.

[jenshnielsen]

However, it's probably not worth pursuing this further since a total replacement with mpreal would be extremely slow

[jenshnielsen]

cmake -DCMAKE_BUILD_TYPE=Release .. -DGMP_INCLUDES=C:\Users\Jens\Documents\GitHub\mpir\lib\x64\Release -DGMP_LIBRARIES=C:\Users\Jens\Documents\GitHub\mpir\lib\x64\Release\mpir.lib -DMPFR_INCLUDES=C:\Users\Jens\Documents\GitHub\mpfr\lib\x64\Release -DMPFR_LIBRARIES=C:\Users\Jens\Documents\GitHub\mpfr\lib\x64\Release\mpfr.lib -G "NMake Makefiles"