search

DCPROGS / HJCFIT

Full maximum likelihood fitting of a mechanism directly to the entire sequence of open and shut times, with exact missed events correction.
GNU General Public License v3.0
7 stars 4 forks source link

Investigating multi precision libraries #129

Open jenshnielsen opened 8 years ago

jenshnielsen commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

jenshnielsen commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

jenshnielsen commented 7 years ago 
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"