R: Цикл векторизации для создания попарной матрицы

Кажется, что невозможно векторизовать исходный фрейм данных df. Но если вы трансформируете его с помощью reshape2::dcast(), чтобы на каждый магазин была одна строчка:

require(reshape2)

df$Fruit <- as.character(df$Fruit)

by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")

#   Shop apple orange pear
# 1    A     1      2    3
# 2    B    NA      1    2
# 3    C     2     NA    1
# 4    D     2      3    1
# 5    E     1      1    1

..., то вы можете легко векторизовать по крайней мере для каждой комбинации [m, n]:

fruits <- unique(df$Fruit)
outer(fruits, fruits, 
    Vectorize(
        function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE), 
        c("m", "n")
    ), 
    by_shop)
#      [,1] [,2] [,3]
# [1,]    0    0    2
# [2,]    2    0    1
# [3,]    1    2    0

Вероятно, это решение, которое вы хотели использовать с outer. Гораздо более быстрым решением была бы настоящая векторизация всех комбинаций фруктов [m, n], но я думал об этом и не вижу способа сделать это. Поэтому мне пришлось использовать функцию Vectorize, которая, конечно, намного медленнее, чем настоящая векторизация.

Контрольное сравнение с вашей исходной функцией:

Unit: milliseconds
                  expr      min       lq     mean   median       uq      max neval
     loop.function(df) 3.788794 3.926851 4.157606 4.002502 4.090898 9.529923   100
 loop.function.TMS(df) 1.582858 1.625566 1.804140 1.670095 1.756671 8.569813   100

Код функции и теста (также добавлено сохранение dimnames):

require(reshape2)   
loop.function.TMS <- function(df) { 
    df$Fruit <- as.character(df$Fruit)
    by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")
    fruits <- unique(df$Fruit)
    o <- outer(fruits, fruits, Vectorize(function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE), c("m", "n")), by_shop)
    colnames(o) <- rownames(o) <- fruits
    o
}

require(microbenchmark)
microbenchmark(loop.function(df), loop.function.TMS(df))

data.table решение:

library(data.table)
setDT(df)
setkey(df,Shop)
dcast(df[df,on=.(Shop=Shop),allow.cartesian=T][
           ,.(cnt=sum(i.Order<Order&i.Fruit!=Fruit)),by=.(Fruit,i.Fruit)]
      ,Fruit~i.Fruit,value.var='cnt')

    Fruit apple orange pear
1:  apple     0      0    2
2: orange     2      0    1
3:   pear     1      2    0

В этом примере индекс Shop не нужен, но он, вероятно, повысит производительность для большего набора данных.

Поскольку вопрос вызвал много комментариев по поводу производительности, я решил проверить, что может принести Rcpp:

library(Rcpp)
cppFunction('NumericMatrix rcppPair(DataFrame df) {

std::vector<std::string> Shop = Rcpp::as<std::vector<std::string> >(df["Shop"]);
Rcpp::NumericVector Order = df["Order"];
Rcpp::StringVector Fruit = df["Fruit"];
StringVector FruitLevels = sort_unique(Fruit);
IntegerVector FruitInt = match(Fruit, FruitLevels);
int n  = FruitLevels.length();

std::string currentShop = "";
int order, fruit, i, f;

NumericMatrix result(n,n);
NumericVector fruitOrder(n);

for (i=0;i<Fruit.length();i++){
    if (currentShop != Shop[i]) {
       //Init counter for each shop
       currentShop = Shop[i];
       std::fill(fruitOrder.begin(), fruitOrder.end(), 0);
    }
    order = Order[i];
    fruit = FruitInt[i];
    fruitOrder[fruit-1] = order;
    for (f=0;f<n;f++) {
       if (order > fruitOrder[f] & fruitOrder[f]>0 ) { 
         result(fruit-1,f) = result(fruit-1,f)+1; 
    }
  }
}
rownames(result) = FruitLevels;
colnames(result) = FruitLevels;
return(result);
}
')

rcppPair(df)

       apple orange pear
apple      0      0    2
orange     2      0    1
pear       1      2    0

В примере набора данных это выполняется ›в 500 раз быстрее, чем решение data.table, вероятно, потому, что оно не имеет проблемы с декартовым произведением. Это не должно быть устойчивым при неправильном вводе и предполагает, что магазины / заказ находятся в порядке возрастания.

Учитывая несколько минут, потраченных на поиск 3 строк кода для data.table решения, по сравнению с гораздо более длительным Rcpp процессом решения / отладки, я бы не рекомендовал использовать Rcpp здесь, если нет реального узкого места в производительности.

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

Вот подход, который вносит простые изменения, чтобы сделать его в 5 раз быстрее.

loop.function2 <- function(df){

    spl_df = split(df[, c(1L, 3L)], df[[2L]])
    
    mat <- array(0L,
                 dim=c(length(spl_df), length(spl_df)),
                 dimnames = list(names(spl_df), names(spl_df)))
    
    for (m in 1:(length(spl_df) - 1L)) {
        xm = spl_df[[m]]
        mShop = xm$Shop
        for (n in ((1+m):length(spl_df))) {
            xn = spl_df[[n]]
            mm = match(mShop, xn$Shop)
            inds = which(!is.na(mm))
            mOrder = xm[inds, "Order"]
            nOrder = xn[mm[inds], "Order"]

            mat[m, n] <- sum(nOrder < mOrder)
            mat[n, m] <- sum(mOrder < nOrder)
        }
    }
    mat
}

Есть 3 основных понятия:

1. Исходные df[df$Fruits == fruits[m], ] строки были неэффективными, поскольку вы бы проводили такое же сравнение length(Fruits)^2 раз. Вместо этого мы можем использовать split(), что означает, что мы сканируем Фрукты только один раз.
2. Очень часто использовался df$var, который извлекает вектор во время каждого цикла. Здесь мы размещаем присваивание xm за пределами внутреннего цикла и пытаемся минимизировать то, что нам нужно для подмножества / извлечения.
3. Я изменил его, чтобы он был ближе к combn, так как мы можем повторно использовать наше match() условие, выполнив оба sum(xmOrder > xnOrder), а затем переключив его на sum(xmOrder < xnOrder).

Представление:

bench::mark(loop.function(df), loop.function2(df))

# A tibble: 2 x 13
##  expression              min median
##  <bch:expr>         <bch:tm> <bch:>
##1 loop.function(df)    3.57ms 4.34ms
##2 loop.function2(df)  677.2us 858.6us

Я догадываюсь, что для вашего более крупного набора данных от @ Waldi data.table будет быстрее. Но для небольших наборов данных это должно быть достаточно эффективным.

Наконец, вот еще один подход rcpp, который кажется медленнее, чем @Waldi:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerMatrix loop_function_cpp(List x) {
    int x_size = x.size();
    IntegerMatrix ans(x_size, x_size);
    
    for (int m = 0; m < x_size - 1; m++) {
        DataFrame xm = x[m];
        CharacterVector mShop = xm[0];
        IntegerVector mOrder = xm[1];
        int nrows = mShop.size();
        for (int n = m + 1; n < x_size; n++) {
            DataFrame xn = x[n];
            CharacterVector nShop = xn[0];
            IntegerVector nOrder = xn[1];
            for (int i = 0; i < nrows; i++) {
                for (int j = 0; j < nrows; j++) {
                    if (mShop[i] == nShop[j]) {
                        if (mOrder[i] > nOrder[j])
                           ans(m, n)++;
                        else
                            ans(n, m)++;
                        break;
                    }
                }
            }
        }
    }
    return(ans);
}

loop_wrapper = function(df) {
  loop_function_cpp(split(df[, c(1L, 3L)], df[[2L]]))
}
loop_wrapper(df)
``

Хорошо, вот решение:

library(tidyverse)

# a dataframe with all fruit combinations
df_compare <-  expand.grid(row_fruit = unique(df$Fruit)
                           , column_fruit = unique(df$Fruit)
                           , stringsAsFactors = FALSE)

df_compare %>%
    left_join(df, by = c("row_fruit" = "Fruit")) %>%
    left_join(df, by = c("column_fruit" = "Fruit")) %>%
    filter(Shop.x == Shop.y &
               Order.x < Order.y) %>%
    group_by(row_fruit, column_fruit) %>%
    summarise(obs = n()) %>%
    pivot_wider(names_from = row_fruit, values_from = obs) %>%
    arrange(column_fruit) %>%
    mutate_if(is.numeric, function(x) replace_na(x, 0)) %>%
    column_to_rownames("column_fruit") %>%
    as.matrix()

       apple orange pear
apple      0      0    2
orange     2      0    1
pear       1      2    0

Если вы не знаете, что происходит во второй части кода (df_compare %>% ...), прочитайте вертикальную черту (%>%) как «затем». Запустите код от df_compare до любого канала, чтобы увидеть промежуточные результаты.