Использование повышенной точности с RcppEigen - PullRequest
0 голосов
/ 18 октября 2018

Я написал библиотеку, используя RcppArmadillo.Проблема, с которой я столкнулся, заключается в том, что для некоторых параметров функция arma :: solve не дает мне точного решения, поскольку у матрицы "A" r-секунда близка к 0. Если arma :: solve может точно решить линейное уравнение, тоне будет проблемой.Но это дает мне приблизительное решение, которое мне не подходит.

Затем я подумал об использовании RcppEigen и использовании переменных повышенной точности.Если я правильно понимаю Eigen, решатель Eigen даст мне решение в мультиточности, и вполне вероятно, что это решение будет достаточно для меня (даже с float128).

Но когда я пытаюсь реализовать этопланирую, у меня ошибка, и я не знаю что делать.Например, работает следующее:

#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <cmath>
#include <cstdint>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/eigen.hpp>

// Correctly setup the build environment 
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(BH)]]
using namespace Rcpp;
using namespace arma;
using namespace Eigen;
using namespace boost;
using namespace boost::multiprecision;
using Eigen::Map;
using Eigen::MatrixXd;                  // variable size matrix, double precision
using Eigen::VectorXd;
using Eigen::Matrix;
using Eigen::Dynamic;
namespace mp = boost::multiprecision;


    // [[Rcpp::export]]
    Eigen::MatrixXd onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
      //Definitions & pre-computations:
      unsigned int tdim;
      double ttc, powttc, wlog, sinwlog, coswlog;
      tdim = t2 - t1 + 1;
      ttc = tc - t1 + 1;
      Eigen::MatrixXd output(tdim, 4); //first dimension is time, second is onesfgh
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = pow(ttc, m);
        wlog = w * log(ttc);
        sinwlog = sin(wlog);
        coswlog = cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

Но:

// [[Rcpp::export]]
Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
  //Definitions & pre-computations:
  unsigned int tdim;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> tc128(tc), m128(m), w128(w);
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> ttc, powttc, wlog, sinwlog, coswlog;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> t1128(double(t1));
  tdim = t2 - t1 + 1;
  ttc = tc128 - t1128 + 1.0;
  Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> output(tdim, 4); //first dimension is time, second is onesfgh
  //Main calculations:
  for (unsigned int i = 0; i < tdim; i++) {
    output(i, 0) = 1.0;
    ttc--;
    powttc = boost::multiprecision::pow(ttc, m128);
    wlog = w128 * boost::multiprecision::log(ttc);
    sinwlog = boost::multiprecision::sin(wlog);
    coswlog = boost::multiprecision::cos(wlog);
    //sincos(wlog, &sinwlog, &coswlog);
    output(i, 1) = powttc;
    output(i, 2) = coswlog * powttc;
    output(i, 3) = sinwlog * powttc;
  }
  return output;
}

не работает.Я получаю сообщение об ошибке:

RcppExports.cpp:11:15: error: 'boost' was not declared in this scope

В частности, я попробовал совет, приведенный в https://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/high_precision/using_test.html, но без использования typedef, поскольку в другом посте Дирк Эддельбюттель предлагает, что typedefs не легкоиметь дело (и они должны быть помещены в другой файл .h).

Есть ли какие-либо предложения, как действовать?

РЕДАКТИРОВАТЬ:

Iизменили функцию, не экспортируя.Теперь код:

  #include <RcppArmadillo.h>
    #include <RcppEigen.h>
    #include <cmath>
    #include <cstdint>
    #include <boost/multiprecision/float128.hpp>
    #include <boost/multiprecision/eigen.hpp>

    // Correctly setup the build environment 
    // [[Rcpp::depends(RcppArmadillo)]]
    // [[Rcpp::depends(RcppEigen)]]
    // [[Rcpp::depends(BH)]]
    using namespace Rcpp;
    using namespace arma;
    using namespace Eigen;
    using namespace boost;
    using namespace boost::multiprecision;
    using Eigen::Map;
    using Eigen::MatrixXd;                  
    using Eigen::VectorXd;
    using Eigen::Matrix;
    using Eigen::Dynamic;
    namespace mp = boost::multiprecision;

    Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int `t2, int t1, double tc, double m, double w) {`
      //Definitions & pre-computations:
      unsigned int tdim;
      mp::float128 tc128(tc), m128(m), w128(w);
      mp::float128 ttc, powttc, wlog, sinwlog, coswlog;
      mp::float128 t1128(double(t1));
      tdim = t2 - t1 + 1;
      ttc = tc128 - t1128 + 1.0;
      Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> output(tdim, 4); 
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = mp::pow(ttc, m128);
        wlog = w128 * mp::log(ttc);
        sinwlog = mp::sin(wlog);
        coswlog = mp::cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

И теперь я получаю следующее сообщение об ошибке:

babel_RcppArmadillo.cpp:6:42: fatal error: boost/multiprecision/eigen.hpp: 
No such file or directory
#include <boost/multiprecision/eigen.hpp>                                          
compilation terminated.

Что странно, так как я предполагаю, что BH включает подпапки boost.

1 Ответ

0 голосов
/ 19 октября 2018

Обновление: Это уже не нужно, поскольку пакет BH включает в себя требуемый заголовок с января 2019 года.


Заголовок boost/multiprecision/eigen.hppбыл добавлен в версии 1.68 , в то время как пакет BH в настоящее время обеспечивает повышение 1.66.Вы должны установить обновленную версию отдельно.Следующее должно работать, но не проверено:

  1. Загрузите boost 1.68 и распакуйте его в какой-нибудь подходящий каталог.В Linux и других Unix-подобных системах я бы, вероятно, использовал /usr/local/include.В противном случае я бы использовал любой путь без пробелов в имени.
  2. Не зависит от пакета BH, т.е. удалите // [[Rcpp::depends(BH)]] для sourceCpp и др.или Imports: BH для использования пакета.
  3. Если boost был установлен в нестандартном месте на шаге 1., вы должны указать компилятору, чтобы он там смотрел.Это можно сделать с помощью

    PKG_CPPFLAGS=-I<path-to-boost>
    

    либо через Sys.setenv для sourceCpp и др.или в пределах src/Makevars(.win) для использования пакета.

...