Я написал рекурсивную функцию на R и использовал мемоиз, чтобы ускорить ее. Я попытался еще больше ускорить его, написав в Rcpp, а затем запомнив функцию Rcpp, но функция R работает быстрее. Почему это так и есть ли способ ускорить это в моем использовании?
require(microbenchmark)
require(Rcpp)
require(memoise)
Функция Rcpp:
cppFunction('
double FunCpp (unsigned int i, double d1, double d2,
double p, double s, unsigned int imax,
double n, double k, double r,
double m, double t) {
if (i == 0) return 0;
if (i == 1) return log2(-1*d1);
if (i == 2) return log2(d2*d1 - p*s);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunCpp(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunCpp(i-2, d1, d2, p, s, imax, n, k, r, m, t);
return x + log2(1 - pow(2,y-x));
}
')
FunCpp = memoise(FunCpp)
Функция R:
FunR = memoise(function(i, d1, d2, p, s, imax, n, k, r, m, t) {
if(i == 0) 0
else if(i == 1) log2(-1*d1)
else if(i == 2) log2(d2*d1 - p*s)
else {
x = log2(abs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)))
x = x + FunR(i-1, d1, d2, p, s, imax, n, k, r, m, t)
y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax))
y = y + FunR(i-2, d1, d2, p, s, imax, n, k, r, m, t)
x + log2(1 - 2^(y-x))
}
})
Это сравнение скорости реалистично для моих целей. Рекурсивная функция используется для диапазона целых чисел, но после этого она больше не будет вызываться с теми же входными данными. Итак, для сравнения скорости здесь я вызываю функции из других функций, и после того, как я закончил вызов рекурсивных функций, я использую Forgot () для сброса кеша.
TestFunR = function() {
x = sapply(1:31, function(i) {
FunR(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
forget(FunR)
}
TestFunCpp = function() {
x = sapply(1:31, function(i) {
FunCpp(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
forget(FunCpp)
}
microbenchmark(TestFunR(), TestFunCpp())
Unit: milliseconds
expr min lq mean median uq max neval cld
TestFunR() 9.979738 10.4910 12.83228 10.91887 11.89264 61.61513 100 a
TestFunCpp() 520.955483 528.6965 547.31103 536.73058 547.66377 729.57631 100 b
Изменить: я получил метод, работающий из книги Дирка, до публикации этого.
includeText = '
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>
class F {
public:
F(unsigned int n = 200, double d1 = 0, double d2 = 0, double p = 0, double s = 0) {
memo.resize(n);
std::fill( memo.begin(), memo.end(), NAN );
memo[0] = 0;
memo[1] = log2(-1*d1);
memo[2] = log2(d2*d1 - p*s);
}
double FunIL(int i, double d1, double d2, double p, double s, double imax,
double n, double k, double r, double m, double t) {
if (i < 0) return((double) NAN);
if (i >= (int) memo.size()) throw std::range_error(\"i too large\");
if (!std::isnan(memo[i])) return(memo[i]);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunIL(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunIL(i-2, d1, d2, p, s, imax, n, k, r, m, t);
memo[i] = x + log2(1 - pow(2,y-x));
return(memo[i]);
}
private:
std::vector< double > memo;
};
'
bodyText = '
int is = Rcpp::as<int>(i);
double d1s = Rcpp::as<double>(d1);
double d2s = Rcpp::as<double>(d2);
double ps = Rcpp::as<double>(p);
double ss = Rcpp::as<double>(s);
double imaxs = Rcpp::as<double>(imax);
double ns = Rcpp::as<double>(n);
double ks = Rcpp::as<double>(k);
double rs = Rcpp::as<double>(r);
double ms = Rcpp::as<double>(m);
double ts = Rcpp::as<double>(t);
F f(ns, d1s, d2s, ps, ss);
return Rcpp::wrap( f.FunIL(is, d1s, d2s, ps, ss, imaxs, ns, ks, rs, ms, ts) );
'
FunInline = cxxfunction(signature(i = "integer", d1 = "numeric", d2 = "numeric", p = "numeric",
s = "numeric", imax = "numeric", n = "numeric", k = "numeric",
r = "numeric", m = "numeric", t = "numeric"),
plugin = "Rcpp",
verbose = T,
incl = includeText,
body = bodyText)
Он работает одинаково хорошо (см. TestFunInline):
microbenchmark(TestFunR(), TestFunCpp(), TestFunCpp_Mem(), TestFunInline())
Unit: microseconds
expr min lq mean median uq max neval cld
TestFunR() 8871.251 9067.758 10301.8003 9287.5725 9593.1310 19270.081 100 b
TestFunCpp() 514415.356 517160.251 522431.2980 519321.6130 523811.7640 584812.731 100 c
TestFunCpp_Mem() 245.474 264.291 284.8908 281.6105 292.0885 526.870 100 a
TestFunInline() 279.686 295.723 378.2134 306.8425 316.0370 6621.364 100 a
Однако мне не удалось заставить его работать с doParallel. Я оптимизирую целевую функцию для каждого процесса, используя пакеты optim и optimx, и когда я использую% do%, он работает, но когда я использую% dopar%, все, что я вижу, это то, что целевую функцию нельзя оценить по начальным параметрам. Я воспользовался советом Дирка из его множества других сообщений и поместил метод Coatless в пакет, но я не уверен, как поместить метод из книги Дирка в пакет. Это просто моя неопытность в C ++.
Изменить 2: Наконец, он щелкнул, как поместить метод Дирка в исходный файл в моем пакете. Я знаю, что есть и другие дискуссии об использовании Rcpp с doParallel, но я помещаю этот код здесь, потому что это еще один хороший способ решить мою проблему, и, добавив этот код в исходный файл в моем пакете, это оказалось намного проще для меня, чтобы заставить это работать в моем параллельном подходе, чем был встроенный.
class F {
public:
F(unsigned int n = 200, double d1 = 0, double d2 = 0, double p = 0, double s = 0) {
memo.resize(n);
std::fill( memo.begin(), memo.end(), NAN );
memo[0] = 0;
memo[1] = log2(-1*d1);
memo[2] = log2(d2*d1 - p*s);
}
double FunIL(int i, double d1, double d2, double p, double s, double imax,
double n, double k, double r, double m, double t) {
if (i < 0) return((double) NAN);
if (i >= (int) memo.size()) throw std::range_error("\"i too large\"");
if (!std::isnan(memo[i])) return(memo[i]);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunIL(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunIL(i-2, d1, d2, p, s, imax, n, k, r, m, t);
memo[i] = x + log2(1 - pow(2,y-x));
return(memo[i]);
}
private:
std::vector< double > memo;
};
// [[Rcpp::export]]
double FunDirk(int i, double d1, double d2, double p, double s,
double imax, double n, double k, double r,
double m, double t) {
F f(n, d1, d2, p, s);
return f.FunIL(i, d1, d2, p, s, imax, n, k, r, m, t);
}
