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

В большом наборе данных (~1M случаев) каждый случай имеет "созданный" и "подвергнутый цензуре" dateTime, Я хочу подсчитать количество других дел, открытых на момент создания каждого дела. Открыты дела между их "созданным" и "цензурированным" dataTimes,

Несколько решений хорошо работают с небольшими наборами данных (<100 000 случаев), но время вычислений растет в геометрической прогрессии. Моя оценка состоит в том, что время вычисления увеличивается как функция 3n^2. При n=100 000 случаев время вычислений>20 минут на моем сервере с ядрами 6 * 4 ГГц и 64 ГБ ОЗУ. Даже с многоядерными библиотеками, в лучшем случае, я мог бы сократить время в 8 или 10 раз. Этого недостаточно для обработки ~ 1 млн. Дел.

Я ищу более эффективный метод для этого расчета. Ниже я предоставил функцию, которая позволяет легко создавать большое количество "созданных" и "цензурированных" dateTime пары вместе с двумя решениями пробовали до сих пор, используя dplyr а также data.table библиотеки. Сроки сообщаются пользователю для простоты. Вы можете просто изменить переменную "CASE_COUNT" в верхней части, чтобы повторно выполнить и снова просмотреть время и легко сравнить время других решений, которые вам, возможно, придется предложить.

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

# Load libraries used in this example
library(dplyr);
library(data.table);
# Not on CRAN. See: http://bioconductor.org/packages/release/bioc/html/IRanges.html
library(IRanges);

# Set seed for reproducibility 
set.seed(123)

# Set number of cases & date range variables
CASE_COUNT  <<- 1000;
RANGE_START <- as.POSIXct("2000-01-01 00:00:00", 
                          format="%Y-%m-%d %H:%M:%S", 
                          tz="UTC", origin="1970-01-01");
RANGE_END   <- as.POSIXct("2012-01-01 00:00:00", 
                          format="%Y-%m-%d %H:%M:%S", 
                          tz="UTC", origin="1970-01-01");

# Select which solutions you want to run in this test           
RUN_SOLUTION_1 <- TRUE;     # dplyr::summarize() + comparisons
RUN_SOLUTION_2 <- TRUE;     # data.table:foverlaps()
RUN_SOLUTION_3 <- TRUE;     # data.table aggregation + comparisons
RUN_SOLUTION_4 <- TRUE;     # IRanges::IRanges + countOverlaps()
RUN_SOLUTION_5 <- TRUE;     # data.table::frank()

# Function to generate random creation & censor dateTime pairs
# The censor time always has to be after the creation time
# Credit to @DirkEddelbuettel for this smart function
# (https://stackoverflow.com/users/143305/dirk-eddelbuettel)

generate_cases_table <- function(n = CASE_COUNT, start_val=RANGE_START, end_val=RANGE_END) {
    # Measure duration between start_val & end_val
    duration <- as.numeric(difftime(end_val, start_val, unit="secs"));

    # Select random values in duration to create start_offset
    start_offset   <- runif(n, 0, duration);

    # Calculate the creation time list
    created_list  <- start_offset + start_val;

    # Calculate acceptable time range for censored values
    # since they must always be after their respective creation value
    censored_range <- as.numeric(difftime(RANGE_END, created_list, unit="secs"));

    # Select random values in duration to create end_offset
    creation_to_censored_times <- runif(n, 0, censored_range);

    censored_list <- created_list + creation_to_censored_times;

    # Create and return a data.table with creation & censor values
    # calculated from start or end with random offsets
    return_table  <- data.table(id       = 1:n,
                                created  = created_list,
                                censored = censored_list);

    return(return_table);
}

# Create the data table with the desired number of cases specified by CASE_COUNT above
cases_table <- generate_cases_table();

solution_1_function <- function (cases_table) { 
    # SOLUTION 1: Using dplyr::summarize:

    # Group by id to set parameters for summarize() function 
    cases_table_grouped <- group_by(cases_table, id);

    # Count the instances where other cases were created before
    # and censored after each case using vectorized sum() within summarize()

    cases_table_summary <- summarize(cases_table_grouped, 
                           open_cases_at_creation = sum((cases_table$created  < created & 
                                                         cases_table$censored > created)));
    solution_1_table <<- as.data.table(cases_table_summary, key="id");        
} # End solution_1_function

solution_2_function <- function (cases_table) {
    # SOLUTION 2: Using data.table::foverlaps:

    # Adapted from solution provided by @Davidarenburg
    # (https://stackoverflow.com/users/3001626/david-arenburg)

    # The foverlaps() solution tends to crash R with large case counts
    # I suspect it has to do with memory assignment of the very large objects
    # It maxes RAM on my system (64GB) before crashing, possibly attempting
    # to write beyond its assigned memory limits.
    # I'll submit a reproduceable bug to the data.table team since
    # foverlaps() is pretty new and known to be occasionally unstable

    if (CASE_COUNT > 50000) {
        stop("The foverlaps() solution tends to crash R with large case counts. Not running.");
    }

    setDT(cases_table)[, created_dupe := created];
    setkey(cases_table, created, censored);

    foverlaps_table  <- foverlaps(cases_table[,c("id","created","created_dupe"), with=FALSE],
                                  cases_table[,c("id","created","censored"),    with=FALSE], 
                                  by.x=c("created","created_dupe"))[order(i.id),.N-1,by=i.id];

    foverlaps_table  <- dplyr::rename(foverlaps_table, id=i.id, open_cases_at_creation=V1);

    solution_2_table <<- as.data.table(foverlaps_table, key="id");
} # End solution_2_function

solution_3_function <- function (cases_table) {    
    # SOLUTION 3: Using data.table aggregation instead of dplyr::summarize

    # Idea suggested by @jangorecki
    # (https://stackoverflow.com/users/2490497/jangorecki)

    # Count the instances where other cases were created before
    # and censored after each case using vectorized sum() with data.table aggregation

    cases_table_aggregated <- cases_table[order(id), sum((cases_table$created  < created & 
                                                     cases_table$censored > created)),by=id];   

    solution_3_table <<- as.data.table(dplyr::rename(cases_table_aggregated, open_cases_at_creation=V1), key="id");

} # End solution_3_function

solution_4_function <- function (cases_table) { 
    # SOLUTION 4: Using IRanges package

    # Adapted from solution suggested by @alexis_laz
    # (https://stackoverflow.com/users/2414948/alexis-laz)

    # The IRanges package generates ranges efficiently, intended for genome sequencing
    # but working perfectly well on this data, since POSIXct values are numeric-representable
    solution_4_table <<- data.table(id      = cases_table$id,
                     open_cases_at_creation = countOverlaps(IRanges(cases_table$created, 
                                                                    cases_table$created), 
                                                            IRanges(cases_table$created, 
                                                                    cases_table$censored))-1, key="id");

} # End solution_4_function

solution_5_function <- function (cases_table) {
    # SOLUTION 5: Using data.table::frank()

    # Adapted from solution suggested by @danas.zuokas
    # (https://stackoverflow.com/users/1249481/danas-zuokas)

    n <- CASE_COUNT;

    # For every case compute the number of other cases
    # with `created` less than `created` of other cases
    r1 <- data.table::frank(c(cases_table[, created], cases_table[, created]), ties.method = 'first')[1:n];

    # For every case compute the number of other cases
    # with `censored` less than `created`
    r2 <- data.table::frank(c(cases_table[, created], cases_table[, censored]), ties.method = 'first')[1:n];

    solution_5_table <<- data.table(id      = cases_table$id,
                     open_cases_at_creation = r1 - r2, key="id");

} # End solution_5_function;

# Execute user specified functions;
if (RUN_SOLUTION_1)
    solution_1_timing <- system.time(solution_1_function(cases_table)); 
if (RUN_SOLUTION_2) {
    solution_2_timing <- try(system.time(solution_2_function(cases_table)));
    cases_table <- select(cases_table, -created_dupe);
}
if (RUN_SOLUTION_3)
    solution_3_timing <- system.time(solution_3_function(cases_table)); 
if (RUN_SOLUTION_4)
    solution_4_timing <- system.time(solution_4_function(cases_table));
if (RUN_SOLUTION_5)
    solution_5_timing <- system.time(solution_5_function(cases_table));         

# Check generated tables for comparison
if (RUN_SOLUTION_1 && RUN_SOLUTION_2 && class(solution_2_timing)!="try-error") {
    same_check1_2 <- all(solution_1_table$open_cases_at_creation == solution_2_table$open_cases_at_creation);
} else {same_check1_2 <- TRUE;}
if (RUN_SOLUTION_1 && RUN_SOLUTION_3) {
    same_check1_3 <- all(solution_1_table$open_cases_at_creation == solution_3_table$open_cases_at_creation);
} else {same_check1_3 <- TRUE;}
if (RUN_SOLUTION_1 && RUN_SOLUTION_4) {
    same_check1_4 <- all(solution_1_table$open_cases_at_creation == solution_4_table$open_cases_at_creation);
} else {same_check1_4 <- TRUE;}
if (RUN_SOLUTION_1 && RUN_SOLUTION_5) {
    same_check1_5 <- all(solution_1_table$open_cases_at_creation == solution_5_table$open_cases_at_creation);
} else {same_check1_5 <- TRUE;}
if (RUN_SOLUTION_2 && RUN_SOLUTION_3 && class(solution_2_timing)!="try-error") {
    same_check2_3 <- all(solution_2_table$open_cases_at_creation == solution_3_table$open_cases_at_creation);
} else {same_check2_3 <- TRUE;}
if (RUN_SOLUTION_2 && RUN_SOLUTION_4 && class(solution_2_timing)!="try-error") {
    same_check2_4 <- all(solution_2_table$open_cases_at_creation == solution_4_table$open_cases_at_creation);
} else {same_check2_4 <- TRUE;}
if (RUN_SOLUTION_2 && RUN_SOLUTION_5 && class(solution_2_timing)!="try-error") {
    same_check2_5 <- all(solution_2_table$open_cases_at_creation == solution_5_table$open_cases_at_creation);
} else {same_check2_5 <- TRUE;}
if (RUN_SOLUTION_3 && RUN_SOLUTION_4) {
    same_check3_4 <- all(solution_3_table$open_cases_at_creation == solution_4_table$open_cases_at_creation);
} else {same_check3_4 <- TRUE;}
if (RUN_SOLUTION_3 && RUN_SOLUTION_5) {
    same_check3_5 <- all(solution_3_table$open_cases_at_creation == solution_5_table$open_cases_at_creation);
} else {same_check3_5 <- TRUE;}
if (RUN_SOLUTION_4 && RUN_SOLUTION_5) {
    same_check4_5 <- all(solution_4_table$open_cases_at_creation == solution_5_table$open_cases_at_creation);
} else {same_check4_5 <- TRUE;}


same_check    <- all(same_check1_2, same_check1_3, same_check1_4, same_check1_5,
                     same_check2_3, same_check2_4, same_check2_5, same_check3_4,
                     same_check3_5, same_check4_5);

# Report summary of results to user
cat("This execution was for", CASE_COUNT, "cases.\n",
    "It is", same_check, "that all solutions match.\n");
if (RUN_SOLUTION_1)
    cat("The dplyr::summarize() solution took", solution_1_timing[3], "seconds.\n");
if (RUN_SOLUTION_2 && class(solution_2_timing)!="try-error")
    cat("The data.table::foverlaps() solution took", solution_2_timing[3], "seconds.\n");
if (RUN_SOLUTION_3)
    cat("The data.table aggregation solution took", solution_3_timing[3], "seconds.\n");
if (RUN_SOLUTION_4)
    cat("The IRanges solution solution took", solution_4_timing[3], "seconds.\n");
if (RUN_SOLUTION_5)
    cat("The data.table:frank() solution solution took", solution_5_timing[3], "seconds.\n\n");

data.table::foverlaps() решение быстрее для меньшего числа случаев (<5000 или около того; зависит от случайности в дополнение к n, поскольку оно использует бинарный поиск для оптимизации). dplyr::summarize() решение быстрее для большего количества случаев (>5000 или около того). За пределами 100 000 ни одно из решений не является жизнеспособным, поскольку оба они слишком медленные.

РЕДАКТИРОВАТЬ: Добавлено третье решение, основанное на идее, предложенной @jangorecki, которая использует data.table агрегация вместо dplyr::summarize()и в остальном аналогичен dplyr решение. Для примерно 50 000 случаев это самое быстрое решение. За 50 000 случаев dplyr::summarize() Решение немного быстрее, но не намного. К сожалению, для 1М случаев это все еще не практично.

EDIT2: добавлено четвертое решение, адаптированное из решения, предложенного @alexis_laz, которое использует IRanges пакет и его countOverlaps функция. Это значительно быстрее, чем 3 других решения. С 50 000 случаев это было почти на 400% быстрее, чем решения 1 и 3.

РЕДАКТИРОВАТЬ3: Модифицированная функция создания регистра для правильного выполнения условия "цензуры". Спасибо @jangorecki за ограничение предыдущей версии.

EDIT4: переписано, чтобы позволить пользователю выбрать, какие решения для выполнения и использовать system.time() для сравнения времени с сборкой мусора перед каждым выполнением для более точной синхронизации (согласно проницательному наблюдению @ jangorecki) - также добавлены некоторые проверки условий для случаев сбоя.

EDIT5: Добавлено пятое решение, адаптированное из решения, предложенного @ danas.zuokas с использованием rank(), Мои эксперименты показывают, что это всегда, по крайней мере, на порядок медленнее, чем другие решения. В 10 000 случаев это занимает 44 секунды против 3,5 секунд для dplyr::summarize и 0,36 секунды для IRanges решения.

ЗАКЛЮЧИТЕЛЬНОЕ РЕДАКТИРОВАНИЕ: Я внес небольшие изменения в решение 5, предложенное @ danas.zuokas, и сопоставил наблюдения @Khashaa с типами. Я установил тип as.numeric в dataTime функция генерации, которая резко ускоряет rank как это действует на integers или же doubles вместо dateTime объекты (увеличивает скорость других функций тоже, но не так резко). После некоторого тестирования, настройки ties.method='first' дает результаты в соответствии с намерением. data.table::frank быстрее, чем оба base::rank а также IRanges::rank, bit64::rank является самым быстрым, но, кажется, обрабатывать связи иначе, чем data.table::frank и я не могу заставить его обращаться с ними так, как хотелось бы. однажды bit64 загружается, маскирует большое количество типов и функций, изменяя результаты data.table::frank по пути. Конкретные причины почему выходят за рамки этого вопроса.

POST END NOTE: Оказывается, data.table::frank ручки POSIXctdateTimes эффективно, тогда как ни base::rank ни IRanges::rank кажется. Таким образом, даже as.numeric (или же as.integer) настройка типа не требуется с data.table::frank и нет потери точности при преобразовании, поэтому их меньше ties.method расхождения. Спасибо всем, кто внес свой вклад! Я многому научился! Очень признателен!:) Кредит будет включен в мой исходный код.

ENDNOTE: Этот вопрос представляет собой уточненную и уточненную версию с более простым в использовании и более удобочитаемым примером кода более эффективного метода подсчета открытых случаев по времени создания каждого случая - я выделил его здесь, чтобы не перегружать исходную статью слишком много правок и для упрощения создания большого количества dataTime пары в примере кода. Таким образом, вам не нужно усердно работать, чтобы ответить. Еще раз спасибо!

2 ответа

Решение

Ответ обновляется в связи с комментарием автора вопроса.

Я бы предложил решение с использованием рангов. Таблицы создаются как в ответ на этот вопрос, или с использованием dateTime функция генерации пар в настоящем вопросе. Оба должны работать.

n <- cases_table[, .N]

# For every case compute the number of other cases
# with `created` less than `creation` of other cases
r1 <- data.table::frank(c(cases_table[, created], cases_table[, created]),
           ties.method = 'first')[1:n]

# For every case compute the number of other cases
# with `censored` less than `created`
r2 <- data.table::frank(c(cases_table[, created], cases_table[, censored]),
           ties.method = 'first')[1:n]

Принимая разницу r1 - r2 (-1 не требуется с связями. Метод ='первый') дает результат (исключая ряды created). С точки зрения эффективности требуется только нахождение рангов вектора длины этого числа строк в cases_table, data.table::frank ручки POSIXctdateTime объекты так быстро, как numeric объекты (в отличие от base::rank), поэтому преобразование типов не требуется.

Это, вероятно, не будет точно отвечать на ваш вопрос, так как воспроизводимый пример не подвергается cases_table$censored > created состояние, см min а также max ниже. Создание меньшего примера поможет вам обнаружить такие проблемы. Также вы должны использовать set.seed в вашем примере.

set.seed(123)
library(data.table)
CASE_COUNT  <- 1000;
RANGE_START <- as.POSIXct("2000-01-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC", origin="1970-01-01");
RANGE_END   <- as.POSIXct("2012-01-01 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="UTC", origin="1970-01-01");
generate_cases_table <- function(n = CASE_COUNT, start=RANGE_START, end=RANGE_END) {
    half_duration <- as.numeric(difftime(end, start, unit="sec")) / 2;
    start_offset  <- runif(n, 0, half_duration);
    end_offset    <- runif(n, 0, half_duration);
    data.table(id       = 1:n,created  = start + start_offset,censored = end   - end_offset)
}
cases_table = generate_cases_table()

cases_table[, .(min_censored = min(censored), max_created = max(created))]
#          min_censored         max_created
#1: 2006-01-01 13:02:12 2005-12-30 04:40:49

setorder(cases_table, created)[, created_so_far := .I - 1L]
cases_table[, censored_after := cases_table[cases_table, on = c("created" = "censored"), roll = Inf, which = TRUE]]

roll Возможно, необходимо изменить соединение, но я не смог проверить его из-за упомянутой проблемы с примерами данных.
which Аргумент - это просто извлечение номера строки из скользящего соединения, для отсортированных данных это также означает количество наблюдений после соединения. Упомянутая проблема приводит к тому, что значение всегда 1000 вызывает все created меньше чем censored,
Подробное описание переходов в data.table см. В этом посте: http://gormanalysis.com/r-data-table-rolling-joins/
Как только вам удастся применить это решение, пожалуйста, поделитесь разницей во времени в комментариях

Другие вопросы по тегам