Гипотеза Коллатца в R

Я до сих пор преподаю немного R в основном себе (и своим ученикам).

Вот реализация последовательности Коллатца в R:

f <- function(n)
{
    # construct the entire Collatz path starting from n
    if (n==1) return(1)
    if (n %% 2 == 0) return(c(n, f(n/2)))
    return(c(n, f(3*n + 1)))
}

При вызове f(13) я получаю 13, 40, 20, 10, 5, 16, 8, 4, 2, 1

Однако обратите внимание, что вектор здесь динамически увеличивается в размере. Такие шаги, как правило, являются рецептом неэффективного кода. Есть ли более производительная версия?

В Python я бы использовал

def collatz(n):
    assert isinstance(n, int)
    assert n >= 1

    def __colla(n):

        while n > 1:
            yield n

            if n % 2 == 0:
                n = int(n / 2)
            else:
                n = int(3 * n + 1)

        yield 1

    return list([x for x in __colla(n)])

Я нашел способ записывать в векторы без указания их размерности априори. Поэтому решение может быть

collatz <-function(n)
{
  stopifnot(n >= 1)  
  # define a vector without specifying the length
  x = c()

  i = 1
  while (n > 1)
  {
    x[i] = n
    i = i + 1
    n = ifelse(n %% 2, 3*n + 1, n/2)
  }
  x[i] = 1
  # now "cut" the vector
  dim(x) = c(i)
  return(x)
}

person tschm    schedule 03.10.2018    source источник
comment
a-look-ar- векторизация через гипотезу Коллатца   -  person Patrick Artner    schedule 03.10.2018
comment
r-and-the-collatz- предположение-часть-2.html   -  person Patrick Artner    schedule 03.10.2018
comment
google R collatz efficient - первые 2 попадания - должны помочь :) понятия не имею о R   -  person Patrick Artner    schedule 03.10.2018
comment
К сожалению, не совсем. Я уже видел эти ссылки. Они даже не пытаются явно построить последовательность чисел в гипотезе. Они просто обновляют счетчик, чтобы получить длину последовательности (ей). Они бегут за вектором аргументов. Тем не менее, очень интересно. Спасибо...   -  person tschm    schedule 04.10.2018


Ответы (1)


Мне было любопытно посмотреть, как реализация C++ через Rcpp будет сравниваться с двумя вашими базовыми подходами R. Вот мои результаты.

Сначала давайте определим функцию collatz_Rcpp, которая возвращает последовательность Hailstone для заданного целого числа n. (Нерекурсивная) реализация была адаптирована из Rosetta Code.

library(Rcpp)
cppFunction("
    std::vector<int> collatz_Rcpp(int i) {
        std::vector<int> v;
        while(true) {
            v.push_back(i);
            if (i == 1) break;
            i = (i % 2) ? (3 * i + 1) : (i / 2);
        }
        return v;
    }
")

Теперь мы запускаем анализ microbenchmark, используя как вашу базовую R, так и реализацию Rcpp. Мы вычисляем последовательности града для первых 10000 целых чисел

# base R implementation
collatz_R <- function(n) {
    # construct the entire Collatz path starting from n
    if (n==1) return(1)
    if (n %% 2 == 0) return(c(n, collatz(n/2)))
    return(c(n, collatz(3*n + 1)))
}

# "updated" base R implementation
collatz_R_updated <-function(n) {
  stopifnot(n >= 1)
  # define a vector without specifying the length
  x = c()
  i = 1
  while (n > 1) {
    x[i] = n
    i = i + 1
    n = ifelse(n %% 2, 3*n + 1, n/2)
  }
  x[i] = 1
  # now "cut" the vector
  dim(x) = c(i)
  return(x)
}

library(microbenchmark)
n <- 10000
res <- microbenchmark(
    baseR = sapply(1:n, collatz_R),
    baseR_updated = sapply(1:n, collatz_R_updated),
    Rcpp = sapply(1:n, collatz_Rcpp))

res
#         expr        min         lq       mean     median         uq       max
#        baseR   65.68623   73.56471   81.42989   77.46592   83.87024  193.2609
#baseR_updated 3861.99336 3997.45091 4240.30315 4122.88577 4348.97153 5463.7787
#         Rcpp   36.52132   46.06178   51.61129   49.27667   53.10080  168.9824

library(ggplot2)
autoplot(res)

введите здесь описание изображения

(Нерекурсивная) реализация Rcpp кажется примерно на 30% быстрее, чем исходная (рекурсивная) базовая реализация R. «Обновленная» (нерекурсивная) базовая реализация R значительно медленнее, чем исходный (рекурсивный) базовый подход R (на моем MacBook Air microbenchmark требуется около 10 минут из-за baseR_updated).

person Maurits Evers    schedule 03.10.2018
comment
Вау, впечатляющие результаты. Большое спасибо. Теперь я стараюсь держаться подальше от C++, так как он добавляет еще один уровень сложности. Я очень удивлен различиями между baseR и baseR_updated. Я предположил, что baseR_updated будет быстрее, чем baseR. Таким образом, baseR выигрывает как по эффективности, так и по элегантности (по сравнению с baseR_updated...) - person tschm; 04.10.2018
comment
Спасибо за интересный вопрос @tschm. Я предполагаю, что вывод состоит в том, что рекурсивный алгоритм (baseR) очень быстр, несмотря на динамическое увеличение вектора в R. Подход Rcpp не является рекурсивным и значительно быстрее, чем эквивалентный базовый подход R baseR_updated. Рекурсивная Rcpp версия, вероятно, была бы победителем ;-) - person Maurits Evers; 07.10.2018
comment
Я предполагаю, что следующей остановкой будет кэширование, в частности, если вы вычислите 10000 последовательностей, как вы это делаете.... - person tschm; 07.10.2018
comment
Я предполагаю, что кеширование функции в R не является тривиальным. Я мог бы сделать m_f‹-memoise(f), но в этом случае рекурсивные вызовы f в f используют кэширование... - person tschm; 07.10.2018