Оптимизация алгоритма Верхоффа в R

Я написал следующую функцию для вычисления контрольной цифры в R.

verhoeffCheck <- function(x)
{
## calculates check digit based on Verhoeff algorithm
## note that due to the way strsplit works, to call for vector x, use sapply(x,verhoeffCheck)

## check for string since leading zeros with numbers will be lost
if (class(x)!="character"){stop("Must enter a string")}

#split and convert to numbers
digs <- strsplit(x,"")[[1]]
digs <- as.numeric(digs)

digs <- rev(digs)   ## right to left algorithm

## tables required for D_5 group

d5_mult <- matrix(c(
                 0:9,
                 c(1:4,0,6:9,5),
                 c(2:4,0:1,7:9,5:6),
                 c(3:4,0:2,8:9,5:7),
                 c(4,0:3,9,5:8),
                 c(5,9:6,0,4:1),
                 c(6:5,9:7,1:0,4:2),
                 c(7:5,9:8,2:0,4:3),
                 c(8:5,9,3:0,4),
                 9:0
                 ),10,10,byrow=T)

d5_perm <- matrix(c(
                 0:9,
                 c(1,5,7,6,2,8,3,0,9,4),
                 c(5,8,0,3,7,9,6,1,4,2),
                 c(8,9,1,6,0,4,3,5,2,7),
                 c(9,4,5,3,1,2,6,8,7,0),
                 c(4,2,8,6,5,7,3,9,0,1),
                 c(2,7,9,3,8,0,6,4,1,5),
                 c(7,0,4,6,9,1,3,2,5,8)
                 ),8,10,byrow=T)

d5_inv <- c(0,4:1,5:9)

## apply algoritm - note 1-based indexing in R
d <- 0

for (i in 1:length(digs)){
    d <- d5_mult[d+1,(d5_perm[(i%%8)+1,digs[i]+1])+1]
    }

d5_inv[d+1]
}

Чтобы работать с вектором строк, необходимо использовать sapply. Отчасти это связано с использованием strsplit, возвращающего список векторов. Это влияет на производительность даже для входных данных среднего размера.

Как эту функцию можно векторизовать?

Я также знаю, что некоторая производительность теряется из-за необходимости создавать таблицы на каждой итерации. Будет ли хранение их в новой среде лучшим решением?


person James    schedule 13.08.2010    source источник


Ответы (3)


Если ваши входные строки могут содержать разное количество символов, то я не вижу никакого способа обойти вызовы lapply (или эквивалент plyr). Хитрость заключается в том, чтобы переместить их внутрь функции, чтобы verhoeffCheck могла принимать векторные входные данные. Таким образом, вам нужно создать матрицы только один раз.

verhoeffCheckNew <- function(x)
{
## calculates check digit based on Verhoeff algorithm

## check for string since leading zeros with numbers will be lost
  if (!is.character(x)) stop("Must enter a string")

  #split and convert to numbers
  digs <- strsplit(x, "")
  digs <- lapply(digs, function(x) rev(as.numeric(x)))

  ## tables required for D_5 group
  d5_mult <- matrix(c(
                   0:9,
                   c(1:4,0,6:9,5),
                   c(2:4,0:1,7:9,5:6),
                   c(3:4,0:2,8:9,5:7),
                   c(4,0:3,9,5:8),
                   c(5,9:6,0,4:1),
                   c(6:5,9:7,1:0,4:2),
                   c(7:5,9:8,2:0,4:3),
                   c(8:5,9,3:0,4),
                   9:0
                   ),10,10,byrow=T)

  d5_perm <- matrix(c(
                   0:9,
                   c(1,5,7,6,2,8,3,0,9,4),
                   c(5,8,0,3,7,9,6,1,4,2),
                   c(8,9,1,6,0,4,3,5,2,7),
                   c(9,4,5,3,1,2,6,8,7,0),
                   c(4,2,8,6,5,7,3,9,0,1),
                   c(2,7,9,3,8,0,6,4,1,5),
                   c(7,0,4,6,9,1,3,2,5,8)
                   ),8,10,byrow=T)

  d5_inv <- c(0,4:1,5:9)

  ## apply algorithm - note 1-based indexing in R      
  sapply(digs, function(x)
  {
    d <- 0  
    for (i in 1:length(x)){
        d <- d5_mult[d + 1, (d5_perm[(i %% 8) + 1, x[i] + 1]) + 1]
        }  
    d5_inv[d+1]
  })
}

Поскольку d зависит от того, что было раньше, не так просто векторизовать цикл for.

Моя версия работает примерно вдвое быстрее для строк 1e5.

rand_string <- function(n = 12) 
{
  paste(sample(as.character(0:9), sample(n), replace = TRUE), collapse = "")
}
big_test <- replicate(1e5, rand_string())

tic()
res1 <- unname(sapply(big_test, verhoeffCheck))
toc()

tic()
res2 <- verhoeffCheckNew(big_test)
toc()

identical(res1, res2) #hopefully TRUE!

См. этот вопрос для tic и toc.

Дальнейшие мысли:

Вам может потребоваться дополнительная проверка ввода для "" и других строк, которые возвращают NA при преобразовании в числовой формат.

Поскольку вы имеете дело исключительно с целыми числами, вы можете получить небольшой выигрыш в производительности от их использования, а не от удвоений. (Используйте as.integer вместо as.numeric и добавьте L к значениям в ваших матрицах.)

person Richie Cotton    schedule 13.08.2010
comment
Очень хорошо! Я нахожу подобное ускорение. Обтекание последнего lapply в as.numeric гарантирует, что будет возвращен вектор, а не список. - person James; 13.08.2010
comment
@James: использование sapply вместо lapply сделает это за вас (без необходимости as.numeric). - person Richie Cotton; 13.08.2010

Начнем с определения матриц поиска. Я разместил их таким образом, чтобы их было легче сверить со ссылкой, например. http://en.wikipedia.org/wiki/Verhoeff_algorithm.

d5_mult <- matrix(as.integer(c(
  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
  1, 2, 3, 4, 0, 6, 7, 8, 9, 5,
  2, 3, 4, 0, 1, 7, 8, 9, 5, 6,
  3, 4, 0, 1, 2, 8, 9, 5, 6, 7,
  4, 0, 1, 2, 3, 9, 5, 6, 7, 8,
  5, 9, 8, 7, 6, 0, 4, 3, 2, 1,
  6, 5, 9, 8, 7, 1, 0, 4, 3, 2,
  7, 6, 5, 9, 8, 2, 1, 0, 4, 3,
  8, 7, 6, 5, 9, 3, 2, 1, 0, 4,
  9, 8, 7, 6, 5, 4, 3, 2, 1, 0
)), ncol = 10, byrow = TRUE)

d5_perm <- matrix(as.integer(c(
  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
  1, 5, 7, 6, 2, 8, 3, 0, 9, 4,
  5, 8, 0, 3, 7, 9, 6, 1, 4, 2,
  8, 9, 1, 6, 0, 4, 3, 5, 2, 7,
  9, 4, 5, 3, 1, 2, 6, 8, 7, 0,
  4, 2, 8, 6, 5, 7, 3, 9, 0, 1,
  2, 7, 9, 3, 8, 0, 6, 4, 1, 5,
  7, 0, 4, 6, 9, 1, 3, 2, 5, 8
)), ncol = 10, byrow = TRUE)

d5_inv <- as.integer(c(0, 4, 3, 2, 1, 5, 6, 7, 8, 9))

Далее мы определим функцию проверки и попробуем ее с тестовым вводом. Я максимально внимательно следил за выводом в Википедии.

p <- function(i, n_i) {
  d5_perm[(i %% 8) + 1, n_i + 1] + 1
}
d <- function(c, p) {
  d5_mult[c + 1, p]
}

verhoeff <- function(x) {
  #split and convert to numbers
  digs <- strsplit(as.character(x), "")[[1]]
  digs <- as.numeric(digs)
  digs <- rev(digs)   ## right to left algorithm

  ## apply algoritm - note 1-based indexing in R
  c <- 0
  for (i in 1:length(digs)) {
    c <- d(c, p(i, digs[i]))
  }

  d5_inv[c + 1]
}
verhoeff(142857)

## [1] 0

Эта функция принципиально итеративна, поскольку каждая итерация зависит от значения предыдущей. Это означает, что мы вряд ли сможем векторизовать в R, поэтому, если мы хотим векторизовать, нам нужно использовать Rcpp.

Однако, прежде чем мы обратимся к этому, стоит выяснить, можем ли мы сделать начальное разделение быстрее. Сначала мы делаем небольшой микробенчмарк, чтобы увидеть, стоит ли заморачиваться:

library(microbenchmark)
digits <- function(x) {
  digs <- strsplit(as.character(x), "")[[1]]
  digs <- as.numeric(digs)
  rev(digs)
}

microbenchmark(
  digits(142857),
  verhoeff(142857)
)

## Unit: microseconds
##              expr   min    lq median    uq   max neval
##    digits(142857) 11.30 12.01  12.43 12.85 28.79   100
##  verhoeff(142857) 32.24 33.81  34.66 35.47 95.85   100

Похоже на это! На моем компьютере verhoeff_prepare() составляет около 50% времени выполнения. Небольшой поиск в stackoverflow показывает другой подход к преобразованию числа в цифры:

digits2 <- function(x) {
   n <- floor(log10(x))
   x %/% 10^(0:n) %% 10
}
digits2(12345)

## [1] 5 4 3 2 1

microbenchmark(
  digits(142857),
  digits2(142857)
)

## Unit: microseconds
##             expr    min     lq median     uq   max neval
##   digits(142857) 11.495 12.102 12.468 12.834 79.60   100
##  digits2(142857)  2.322  2.784  3.358  3.561 13.69   100

digits2() намного быстрее, чем digits(), но имеет ограниченное влияние на всю среду выполнения.

verhoeff2 <- function(x) {
  digs <- digits2(x)

  c <- 0
  for (i in 1:length(digs)) {
    c <- d(c, p(i, digs[i]))
  }

  d5_inv[c + 1]
}
verhoeff2(142857)

## [1] 0

microbenchmark(
  verhoeff(142857),
  verhoeff2(142857)
)

## Unit: microseconds
##               expr   min    lq median    uq   max neval
##   verhoeff(142857) 33.06 34.49  35.19 35.92 73.38   100
##  verhoeff2(142857) 20.98 22.58  24.05 25.28 48.69   100

Чтобы сделать это еще быстрее, мы могли бы попробовать C++.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int verhoeff3_c(IntegerVector digits, IntegerMatrix mult, IntegerMatrix perm,
                IntegerVector inv) {
  int n = digits.size();
  int c = 0;

  for(int i = 0; i < n; ++i) {
    int p = perm(i % 8, digits[i]);
    c = mult(c, p);
  }

  return inv[c];
}

verhoeff3 <- function(x) {
  verhoeff3_c(digits(x), d5_mult, d5_perm, d5_inv)
}
verhoeff3(142857)

## [1] 3

microbenchmark(
  verhoeff2(142857),
  verhoeff3(142857)
)

## Unit: microseconds
##               expr   min    lq median    uq   max neval
##  verhoeff2(142857) 21.00 22.85  25.53 27.11 63.71   100
##  verhoeff3(142857) 16.75 17.99  18.87 19.64 79.54   100

Это не дает большого улучшения. Возможно, мы сможем добиться большего успеха, если передадим число в C++ и обработаем цифры в цикле:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int verhoeff4_c(int number, IntegerMatrix mult, IntegerMatrix perm,
                IntegerVector inv) {
  int c = 0;
  int i = 0;

  for (int i = 0; number > 0; ++i, number /= 10) {
    int p = perm(i % 8, number % 10);
    c = mult(c, p);
  }

  return inv[c];
}

verhoeff4 <- function(x) {
  verhoeff4_c(x, d5_mult, d5_perm, d5_inv)
}
verhoeff4(142857)

## [1] 3

microbenchmark(
  verhoeff2(142857),
  verhoeff3(142857),
  verhoeff4(142857)
)

## Unit: microseconds
##               expr    min     lq median     uq   max neval
##  verhoeff2(142857) 21.808 24.910 26.838 27.797 64.22   100
##  verhoeff3(142857) 17.699 18.742 19.599 20.764 81.67   100
##  verhoeff4(142857)  3.143  3.797  4.095  4.396 13.21   100

И получаем отдачу: verhoeff4() примерно в 5 раз быстрее, чем verhoeff2().

person hadley    schedule 19.03.2014
comment
Спасибо, Хэдли, хорошая работа! Единственная проблема заключается в том, что, как правило, начальные нули важны, поэтому ввод нельзя преобразовать в числовой, как это делают digits2, verhoeff3 и verhoeff4 с некоторым дополнением векторов разделенных цифр. Кроме того, использование целых чисел наложило бы ограничения на длину ввода: для 32-битных целых чисел было бы безопасно только 8 цифр. Штрих-коды, хотя и используют другую схему, состоят из 12 цифр плюс контрольная цифра. - person James; 25.03.2014
comment
@James, стратегия по-прежнему будет работать с числовыми векторами (для увеличения диапазона), или я ожидаю, что С++ будет почти так же быстро перебирать chars из String - person hadley; 25.03.2014

Richie C хорошо ответил на вопрос о векторизации; что касается создания таблиц только один раз без загромождения глобального пространства имен, одним быстрым решением, не требующим пакета, является

verhoeffCheck <- local(function(x)
{
## calculates check digit based on Verhoeff algorithm
## note that due to the way strsplit works, to call for vector x, use sapply(x,verhoeffCheck)
## check for string since leading zeros with numbers will be lost
if (class(x)!="character"){stop("Must enter a string")}
#split and convert to numbers
digs <- strsplit(x,"")[[1]]
digs <- as.numeric(digs)
digs <- rev(digs)   ## right to left algorithm
## apply algoritm - note 1-based indexing in R
d <- 0
for (i in 1:length(digs)){
    d <- d5_mult[d+1,(d5_perm[(i%%8)+1,digs[i]+1])+1]
    }
d5_inv[d+1]
})

assign("d5_mult", matrix(c(
    0:9, c(1:4,0,6:9,5), c(2:4,0:1,7:9,5:6), c(3:4,0:2,8:9,5:7),
    c(4,0:3,9,5:8), c(5,9:6,0,4:1), c(6:5,9:7,1:0,4:2), c(7:5,9:8,2:0,4:3),
    c(8:5,9,3:0,4), 9:0), 10, 10, byrow = TRUE), 
    envir = environment(verhoeffCheck))

assign("d5_perm", matrix(c(
    0:9, c(1,5,7,6,2,8,3,0,9,4), c(5,8,0,3,7,9,6,1,4,2),
    c(8,9,1,6,0,4,3,5,2,7), c(9,4,5,3,1,2,6,8,7,0), c(4,2,8,6,5,7,3,9,0,1),
    c(2,7,9,3,8,0,6,4,1,5), c(7,0,4,6,9,1,3,2,5,8)), 8, 10, byrow = TRUE),
    envir = environment(verhoeffCheck))

assign("d5_inv", c(0,4:1,5:9), envir = environment(verhoeffCheck))
## Now just use the function

который хранит данные в среде функции. Вы можете рассчитать время, чтобы увидеть, насколько это быстрее.

Надеюсь это поможет.

Аллан

person Allan Engelhardt    schedule 13.08.2010