Память функции Rcpp?

Question

Я написал рекурсивную функцию в R и использовал memoise, чтобы ускорить ее. Я попытался ускорить его, записав его в 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)) 
    } 
}) 

Эта скорость сравнения реалистична для моих целей. Рекурсивная функция используется в диапазоне целых чисел, но впоследствии она не будет вызываться с теми же входами снова. Поэтому для сравнения скорости здесь я вызываю функции из других функций, и после того, как я закончил вызов рекурсивных функций, я использую функцию() для сброса кеша.

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%, все, что я вижу, это то, что целевую функцию нельзя оценить при начальных параметрах. Я взял совет Дирка из его многих других сообщений и поместил метод Коатлеса в пакет, но я не уверен, как положить этот метод из книги Дирка в пакет. Это просто моя неопытность в C++.

Редактировать 2: Он, наконец, нажал, как установить метод Dirk в исходный файл в моем пакете. Я знаю, что есть другие обсуждения об использовании Rcpp с doParallel, но я помещаю этот код здесь, потому что это еще один хороший способ решить мою проблему, и, добавив этот код в исходный файл в моем пакете, оказалось, что это намного проще для меня, чтобы это работало в моем параллельном подходе, чем inline.

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); 

} 

Answer 1

Memoisize меня

Ну, во-первых, давайте думать о цели memoise. Целью memoise является Результаты функции кеша и повторное использование. Таким образом, после одного вычисления ему больше не нужно пересчитывать значение снова для любой другой последовательности в вычислении, он может просто извлечь значение из кеша. Это особенно актуально для создания рекурсивной структуры.

memoise доступа к кэшу по отношению к R и C++

установка для memoisize является кэширование значений функции значения R. В этом случае он кэширует эти значения. Однако код C++ не может получить доступ к кешированным значениям. Таким образом, версия C++ пересчитывает каждое из этих значений.В сущности, вы в основном с помощью:

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) 
    }) 

Big O Algo

Предупреждение: Там должно быть более формальный аргумент следующий, но его уже некоторое время.

Чтобы понять алгоритмы, иногда нам нужно использовать то, что называется Big O notation, что позволяет нам наблюдать, как код выполняет асимптотически. Теперь Big O в этом случае является O (2^N) из-за двух вызовов для вычисления: Fun(i-1) и FunR(i-2). Тем не менее, memoise использует таблицу хэша/таблицу поиска с вероятностью Big O O(n) в худшем случае и O(1) в лучшем виде. По существу, мы имеем постоянные и экспоненциальные асимптотические результаты.

Улучшение Microbenchmarks - W/O Memosizing в C++

Однако, это не обязательно означает, что функция С ++ мусора. Одним из недостатков R-Rcpp и обратного моста является время задержки между передачей значений между двумя доменами. Таким образом, одним из способов, которым мы можем немного снизить время вычисления, является полное размещение цикла в C++.

например.

// [[Rcpp::export]] 
Rcpp::NumericVector FunCpp_loop(unsigned int e, 
           double d1, double d2, 
           double p, double s, unsigned int imax, 
           double n, double k, double r, 
           double m, double t){ 

    Rcpp::NumericVector o(e); 

    for(unsigned int i = 0; i < e; i++){ 

    o(i) = FunCpp(31-(i+1), -152, -147.33, 150, 0.03, 30, 31, 1, 1, 2, 5); 

    } 

    return o; 
} 

Однако ориентиры здесь не реально улучшить ситуацию (даже до создания вектора 1:31)

Unit: milliseconds 
       expr  min   lq  mean  median  uq  max neval 
     TestFunR(tv) 8.467568 9.077262 9.986837 9.449952 10.60555 14.91243 100 
    TestFunCpp(tv) 476.678391 482.489094 487.687811 486.351087 490.25346 579.38161 100 
TestFunCpp_loop() 478.348070 482.588307 488.234200 486.211347 492.33965 521.10918 100 

Memoziation в C++

Мы можем применить те же методы memoziation данные в memoise в C++. Реализация не такая красивая и приятная, но она служит для того, чтобы показать, что одни и те же принципы применимы.

Для начала сделаем хэш-карту.

// Memoization structure to hold the hash map 
struct mem_map{ 

    // Initializer to create the static (presistent) map 
    static std::map<int, double> create_map() 
    { 
    std::map<int, double> m; 
    m.clear(); 
    return m; 
    } 

    // Name of the static map for the class 
    static std::map<int, double> memo; 

}; 

// Actuall instantiate the class in the global scope (I know, bad me...) 
std::map<int, double> mem_map::memo = mem_map::create_map(); 

Теперь мы, вероятно, должны сделать некоторые аксессоры для работы с объектом карты.

// Reset the map 
// [[Rcpp::export]] 
void clear_mem(){ 
    mem_map::memo.clear(); 
} 

// Get the values of the map. 
// [[Rcpp::export]] 
std::map<int, double> get_mem(){ 
    return mem_map::memo; 
} 

И, наконец, давайте изменим некоторые внутренние функции в вашей функции.

// Users function 
// [[Rcpp::export]] 
double FunCpp_Mem (int i, double d1, double d2, 
        double p, double s, unsigned int imax, 
        double n, double k, double r, 
        double m, double t) { 

    // We have already computed the value... 
    if(mem_map::memo.count(i) > 0) 
    return mem_map::memo[i]; 


    // Otherwise, let us get ready to compute it! 
    double res = 0; 

    if (i <= 2){ 
    if (i <= 0) { // i == 1 
     res = 0.0; 
    }else if (i == 1) { 
     res = log2(-1.0*d1); 
    }else { // i == 2 
     res = log2(d2*d1 - p*s); 
    } 

    // Store result in hashmap 
    mem_map::memo[i] = res; 

    return res; 
    } 

    // Calculate if not in special case. 

    double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t))); 
    x = x + FunCpp_Mem(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_Mem(i-2, d1, d2, p, s, imax, n, k, r, m, t); 


    res = x + log2(1 - pow(2,y-x)); 


    // Update the hashmap for uncalculated value 
    mem_map::memo[i] = res; 

    return res; 
} 

Уф Много работы. Давайте проверим это, посмотрим, стоит ли это.

# Benchmark for Rcpp Memoization 
TestFunCpp_mem = function(tv) { 
    x = sapply(tv, function(i) { 
    FunCpp_Mem(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) 
    }) 
    clear_mem() 
} 

TestFunR = function(tv) { 
    x = sapply(tv, 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) 
} 

# Pre-generate vector 
tv = 1:31 

microbenchmark(TestFunR(tv),TestFunCpp_mem(tv)) 

И результаты ....

microbenchmark(TestFunR(tv),TestFunCpp_mem(tv)) 
Unit: microseconds 
       expr  min  lq  mean median  uq  max neval 
     TestFunR(tv) 8246.324 8662.694 9345.6947 9009.868 9797.126 13001.995 100 
TestFunCpp_mem(tv) 203.832 214.939 253.7931 228.898 240.906 1277.325 100 

Функция Cpp с запоминанием о 40.5x быстрее, чем R-версии! Память была определенно стоит!