R: эффективное вычисление резюме подмножеств значений, содержание которых определяется соотношением между двумя переменными

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

# Load the package
install.packages("sqldf")
library(sqldf)

# Create tables
A <- data.frame("loc"=c(150,250,400))
B <- data.frame("loc"=c(25,77,170,320),"value"=c(7,19,10,15))


# Join tables
df0 <- sqldf('select a.loc
                    ,count(b.value) as n_value
                    ,sum(b.value) as sum_value
                    ,avg(b.value) as avg_value
              from A as a
              left join B as b
              on abs(a.loc - b.loc) <= 100
              group by a.loc')

# Print data frame
df0

Я не уверен, насколько хорошо это решение будет масштабироваться - это зависит от того, поместится ли матрица фильтра в память.

A <- within(A,{
 B.filter <- outer(B$loc, A$loc, function(x, y) abs(x - y) <= 100) 

 n <- colSums(B.filter)
 sum <- colSums(B$value * B.filter)
 avg <- sum / n
 rm(B.filter)
})

Если местоположения в A и / или B повторяются, вы можете уменьшить размер матрицы фильтра, используя только уникальные значения:

A <- within(A,{
 B.filter <- outer(unique(B$loc), unique(A$loc), function(x, y) abs(x - y) <= 100) 
 colnames(B.filter) <- unique(A$loc)
 rownames(B.filter) <- unique(B$loc)

 n <- colSums(B.filter[,as.character(A$loc)])
 sum <- colSums(B$value * B.filter[as.character(B$loc),])
 avg <- sum / n
 rm(B.filter)
})

Вот решение tidyverse

library(tidyverse)

A = read.table(text = "
loc
150
250
400
", header=T)

B = read.table(text = "
loc    value
25     7
77     19
170    10
320    15
", header=T)

A %>%
  mutate(B = list(B)) %>%              # create all combinations of rows of A and B
  unnest() %>%
  filter(abs(loc - loc1) <= 100) %>%   # keep rows that satisfy your condition
  group_by(loc) %>%                    # for each loc values
  summarise(sum = sum(value),          # calculate sum
            avg = mean(value))         # calculate mean

# # A tibble: 3 x 3
#     loc   sum   avg
#    <int> <int> <dbl>
# 1   150    29  14.5
# 2   250    25  12.5
# 3   400    15  15  

Возможно, это не лучшее решение, если у вас большие таблицы A и B, поскольку вам нужно создать все комбинации строк, а затем отфильтровать.

Мое чистое решение R (ниже) все еще значительно медленное, в моей системе требовалось 32 секунды, чтобы закончить большой пример Мэтта Саммерсгилла, но по сравнению с другими решениями, это все еще разумно.

Логика моего решения такова: поскольку входы отсортированы, когда вы рассматриваете новые элементы A_loc, диапазон значений в B_loc останется прежним, если новый элемент A_loc идентичен предыдущему, или он сместится вправо в B_loc, возможно сокращение или расширение. Обратите внимание, что если вы работали с входами double, вам нужно быть немного осторожнее со сравнениями.

Эта версия C++, естественно, быстрее. Если вы можете Rcpp::sourceCpp этот код:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
DataFrame foo(IntegerVector A_loc, IntegerVector B_loc, IntegerVector B_val) {
    IntegerVector n(A_loc.length());
    IntegerVector sum(A_loc.length());
    NumericVector avg(A_loc.length());

    int lower = 0;
    int upper = 0;
    int count = 0;
    int current_sum = 0;
    for (int i = 0; i < A_loc.length(); i++) {
        checkUserInterrupt();

        while (lower < B_loc.length()) {
            if (B_loc[lower] >= A_loc[i] - 100) {
                break;
            }

            if (count > 0) {
                count--;
                current_sum -= B_val[lower];
            }

            lower++;
        }

        if (upper < lower) {
            upper = lower;
        }

        while (upper < B_loc.length()) {
            if (B_loc[upper] > A_loc[i] + 100) {
                break;
            }

            count++;
            current_sum += B_val[upper++];
        }

        n[i] = count;
        sum[i] = current_sum;
        avg[i] = static_cast<double>(current_sum) / count;
    }

    DataFrame df = DataFrame::create(
        Named("loc") = A_loc,
        Named("n") = n,
        Named("sum") = sum,
        Named("avg") = avg
    );

    return df;
}

тогда это:

A <- data.frame(loc = sample.int(1000, size = 1e4, replace = TRUE))


B <- data.frame(loc = sample.int(1000, size = 1e6, replace = TRUE),
                value = sample.int(100, size = 1e6, replace = TRUE))

test <- function() {
    # remove unique if you want to consider duplicated values
    A_loc <- sort(unique(A$loc), decreasing = FALSE)
    B <- B[order(B$loc),]
    out <- foo(A_loc, B$loc, B$value)
}

microbenchmark::microbenchmark(test())

показывает эти тайминги:

Unit: milliseconds
   expr      min      lq     mean   median       uq      max neval
 test() 44.74469 45.8118 51.35361 47.34657 48.99376 95.00938   100

Если вы не можете использовать Rcpp, затем рассмотрите версию R ниже, или решение Фрэнка с data.table, Я думаю, что сортировка входных данных также может помочь в этом случае?

Циклы for обычно избегают в R, но я не думаю, что они всегда медленные, вам просто нужно быть осторожным, чтобы не копировать слишком много данных. Кроме того, начиная с R v3.5.0, запись чего-то вроде for i in 1:10 больше не выделяет сначала весь вектор, он поддерживает компактное представление.

A_loc <- sort(unique(A$loc), decreasing = FALSE)
B <- B[order(B$loc),]

out <- data.frame(loc = A_loc,
                  n = 0L,
                  sum = 0L,
                  avg = 0)

lower <- 1L
upper <- 1L
count <- 0L
sum <- 0L
upper_limit <- nrow(B)
for (i in seq_along(A_loc)) {
  current_loc <- A_loc[i]

  while (lower <= upper_limit) {
    if (B$loc[lower] >= current_loc - 100L) {
      break
    }

    if (count > 0L) {
      count <- count - 1L
      sum <- sum - B$value[lower]
    }

    lower <- lower + 1L
  }

  if (upper < lower) {
    upper <- lower
  }

  while (upper <= upper_limit) {
    if (B$loc[upper] > current_loc + 100L) {
      break
    }

    count <- count + 1L
    sum <- sum + B$value[upper]
    upper <- upper + 1L
  }

  out$n[i] <- count
  out$sum[i] <- sum
  out$avg[i] <- sum / count
}

Это возможно с помощью функции foverlaps в data.table, и следующий метод действительно требует завершения вашего фактического варианта использования - A - тысячи строк и таблица B - почти миллионы строк. - в разумные сроки.

С вашим примером игрушки:

library(data.table)

A <- fread("
           loc
           150
           250
           400")

B <- fread("
           loc    value
           25     7
           77     19
           170    10
           320    15")

## Create a 'dummy' value to create an interval w/same start and end in A
A[,loc_Dummy := loc]

## Create values bounding the match range for loc in B
B[,loc_Plus100 := loc + 100]
B[,loc_Minus100 := loc - 100]

## Set up for the overlap join
setkey(A,loc,loc_Dummy)
setkey(B,loc_Minus100, loc_Plus100)

## Create a table of only matches instead of doing a full cartesian join of all cases
Matches <- foverlaps(A[,.(loc, loc_Dummy)],
                     B[,.(loc_Minus100,loc_Plus100,value)])

## Create a summary table
Matches[,.(n = .N, sum = sum(value), avg = mean(value)), by = .(loc)]

#    loc n sum  avg
# 1: 150 2  29 14.5
# 2: 250 2  25 12.5
# 3: 400 1  15 15.0

Увеличение масштаба - ура!

Однако на самом деле это проблема с интенсивными вычислениями очень сильно. Масштабирование до ваших фактических размеров случая показывает здесь проблему - используя 10 000 строк для таблицы A и 1 000 000 строк для таблицы B, этот метод завершается в 91 секунда на сервере, на котором я работаю, но использует более 112 ГБ памяти!

A <- data.table(loc = sample.int(1000, size = 1e4, replace = TRUE))


B <- data.table(loc = sample.int(1000, size = 1e6, replace = TRUE),
                value = sample.int(100, size = 1e6, replace = TRUE))

system.time({
  A[,loc_Dummy := loc]
  B[,loc_Plus100 := loc + 100]
  B[,loc_Minus100 := loc - 100]

  setkey(A,loc,loc_Dummy)
  setkey(B,loc_Minus100, loc_Plus100)

  Matches <- foverlaps(A[,.(loc, loc_Dummy)],
                       B[,.(loc_Minus100,loc_Plus100,value)])

  Summary  <- Matches[,.(n = .N, sum = sum(value), avg = mean(value)), by = .(loc)]

})

## Warning: Memory usage peaks at ~112 GB!

# user  system elapsed 
# 56.407  46.258  91.135

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

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

Насколько я могу судить, ваша проблема фактически похожа на проблему, поставленную (и решенную) Лоренцо Бузетто и подробно описанную в сообщении блога: Ускорение пространственного анализа за счет интеграции sf и data.table: тестовый пример.

Кусочки на помощь

Требование более ~ 100 гигабайт памяти не является возможным решением В самом деле, особенно если вы хотите в какой-то момент масштабировать A или B на порядок выше.

Однако метод фрагментации (вдохновленный сообщением Лоренцо, ссылка на который приведена выше), который разбивает проблему на 100 фрагментов, фактически только увеличивает во время выполнения тривиальную величину до 116 секунд, но снижает пиковое использование памяти до менее 3 ГБ! Если бы я планировал сделать это в продакшене, я бы выбрал что-то вроде следующего.

Одно замечание: я действительно не проводил углубленного аудита точности результатов (я мог указать одну из границ диапазона неправильно открыт или закрыт), поэтому я бы внимательно изучил вывод с данными, которые вам знакомы. с перед запуском в производство.

A <- data.table(loc = sample.int(1000, size = 1e4, replace = TRUE))

B <- data.table(loc = sample.int(1000, size = 1e6, replace = TRUE),
                value = sample.int(100, size = 1e6, replace = TRUE))

system.time({

  A[,loc_Dummy := loc]
  B[,loc_Plus100 := loc + 100]
  B[,loc_Minus100 := loc - 100]

  setkey(A,loc)
  setkey(B,loc)

  ChunkCount <- 100
  ChunkSize <- A[,.N/ChunkCount]

  ResultList <- vector("list", ChunkCount) 

  for (j in seq_len(ChunkCount)){

    A_loc_Min <- A[((j-1)*ChunkSize + 1):(min(nrow(A),(j)*ChunkSize)), min(loc)]
    A_loc_Max <- A[((j-1)*ChunkSize + 1):(min(nrow(A),(j)*ChunkSize)), max(loc)]

    A_Sub <- A[loc >= A_loc_Min & loc < A_loc_Max]
    B_Sub <- B[loc_Plus100 >= A_loc_Min & loc_Minus100 < A_loc_Max]

    setkey(A_Sub,loc,loc_Dummy)
    setkey(B_Sub,loc_Minus100, loc_Plus100)

    Matches <- foverlaps(A_Sub[,.(loc, loc_Dummy)],
                         B_Sub[,.(loc_Minus100,loc_Plus100,value)])

    ResultList[[j]]  <- Matches[,.(n = .N, sum = sum(value), avg = mean(value)), by = .(loc)]

  }

  Summary  <- rbindlist(ResultList)

})

#    user  system elapsed 
# 109.125  16.864 116.129 

Проверка

Update: @Alexis and @Frank's suggestion in the comments generate the same result set, mine comes out slightly different, but only on the count. If someone else can verify that the correct answer is actually that provided by @Alexis/@Frank, then I'd be happy to retract my answer as both methods execute faster than the one I proposed.

set.seed(1234)

A <- data.table(loc = sample.int(1000, size = 1e3, replace = TRUE))

B <- data.table(loc = sample.int(1000, size = 1e4, replace = TRUE),
                value = sample.int(10, size = 1e4, replace = TRUE))



## Matt 
Matt_A <- copy(A)
Matt_B <- copy(B)

Matt_A[,loc_Dummy := loc]
Matt_B[,loc_Plus100 := loc + 100]
Matt_B[,loc_Minus100 := loc - 100]

setkey(Matt_A,loc,loc_Dummy)
setkey(Matt_B,loc_Minus100, loc_Plus100)

Matches <- foverlaps(Matt_A[,.(loc, loc_Dummy)],
                     Matt_B[,.(loc_Minus100,loc_Plus100,value)])

Summary_Matt  <- Matches[,.(n = .N, sum = sum(value), avg = mean(value)), keyby = .(loc)]


## Alexis

Rcpp::sourceCpp("RowRanges.cpp")

A_loc <- sort(A$loc, decreasing = FALSE)
B <- B[order(B$loc),]
Alexis <- foo(unique(A_loc), B$loc, B$value)

Summary_Alexis <- as.data.table(Alexis)
colnames(Summary_Alexis) <- c("n","sum","avg")

Summary_Alexis[,loc := unique(A_loc)]
setcolorder(Summary_Alexis, c("loc","n","sum","avg"))

## Frank

Frank <- A[, up := loc + 100][
  , dn := loc - 100][
    , c("n", "s", "m") := B[copy(.SD), on=.(loc >= dn, loc <= up), .(.N, sum(value), mean(value)), by=.EACHI][
      , .(N, V2, V3)]][]

Summary_Frank <- unique(Frank[,.(loc,n, sum = s, avg = m)][order(loc)])

## Comparing

all.equal(Summary_Frank,Summary_Alexis)
# [1] TRUE

all.equal(Summary_Frank,Summary_Matt)
# [1] "Column 'n': Mean relative difference: 1.425292"

Подобно ответу Мэтта Саммерсгилла, вы можете выполнить неэквивалентное соединение для обновления A:

A[, up := loc + 100]
A[, dn := loc - 100]
A[, c("n", "s", "m") := 
  B[copy(.SD), on=.(loc >= dn, loc <= up), .(.N, sum(value), mean(value)), by=.EACHI][, .(N, V2, V3)]
]

Или в одной связанной команде:

A[, up := loc + 100][, dn := loc - 100][, c("n", "s", "m") := 
  B[copy(.SD), on=.(loc >= dn, loc <= up), 
    .(.N, sum(value), mean(value)), by=.EACHI][, 
    .(N, V2, V3)]
]

Думаю, это должно быть достаточно эффективно.

Как это работает

Внутри jx[i, j], .SD относится к подмножеству данных из x (в данном случае это все A).

x[i, on=, j, by=.EACHI] - это соединение, в котором каждая строка i (в данном случае copy(.SD) == A) используется для поиска совпадающих строк x (в данном случае B) с использованием условий в on=. Для каждой строки i вычисляется j (что означает by=.EACHI).

Когда у j нет имен, они присваиваются автоматически. V1, V2 и так далее. .N по умолчанию получает имя N.