«Неравный брак», В. Пукирев, 1862 г.

Задача объединения табличных представлений очень часто встречается как в аналитике, так и в разработке (БД). Существует несколько различных типов слияний, фактически, это операции над множествами. Не будем погружаться в детали, на эту тему написано множество книг, семинаров, публикаций. Посмотрим на эти механизмы в преломлении практических задач. Будем смотреть по нарастающей сложности и пытаться решить их на «офисном» ноутбуке, не привлекая бесконечные мощности больших данных или реляционные БД.

Является продолжением серии предыдущих публикаций.

Задача 1. Корабли и рыбы

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

Есть два датафрейма:

1. точки телеметрии — отметки положения морских животных (меняются во времени и пространстве), их около 10 тысяч.
2. точки — отметки положения движущихся судов в том же районе (неск. млн) за те же месяцы.

Интересует вопрос: избегают ли животные судов или им все равно?

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

1. для каждой точки-отметки животного находим расстояния до N ближайших судов (в тот же самый момент времени) — получаем распределение расстояний;
2. в рассматриваемом районе генерируем случайные точки и находим расстояния до судов от этих случ. точек, получаем второе распределение;
3. сравниваем два распределения, делаем выводы касательно нулевой гипотезы.

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

И вот тут точно начинаются проблемы.
Алгоритм может и правильный для бумаги, но не совсем удачный для кодификации «в лоб». Классический подход «цикл в цикле» на интерпретируемых языках R/Python начинает давать сбой. В ожидании ответа можно провести не один час или день. Небольшая правка в методике расчета или данных — и опять ждем долго-долго. Надо что-то делать.

Для эффективного переложения на ЯП необходимо придерживаться ряда базовых правил:

1. Все что можно считать вне циклов и групп — надо считать вне.
2. Если оперативной памяти хватает, то лучше создать все сразу.

Начнем со второго пункта. Для начала проведем грубую предварительную оценку.

Несколько миллионов отметок по кораблям переведем в удельный суточный показатель.
Допустим, что наблюдения проводились ежедневно в течение полугода (~ 200 суток), а измерений у нас 10^6. Тогда в наблюдаемой зоне находится ~ 10^6 / 200 = 5*10^3. Для простоты предположим, что телеметрию собирают чуть чаще чем раз в сутки, что даст в зоне наблюдения ~ 10^3 кораблей.

10^4 животных * 10^3 кораблей = 10^7. Вполне нормальный показатель для современных ноутбуков. 12 штук int64/float64 колонок (координаты животного, корабля, дата и еще что-нибудь) дадут на строку ~ 100 байт, а на всю таблицу ~ 1 Гб

Итого, датасет должен получиться всего 1 Гб на один день. Можно практически поиграться, получим ~ 0.9 Гб с учетом степеней двойки:

data.table::CJ(a = as.numeric(1:10^3), b = as.numeric(1:10^4)) %>%
      pryr::object_size() * 6 %>%
      fs::fs_bytes()

Итого, штатный алгоритм такой:

1. Генерим реперный каркас («гребенка») для проведения расчетов в виде
   data.frame с перебором всех комбинаций {дата, id_животного, id_корабля}. Размах гребенки по датам будет определяться размером доступной RAM.
   Сама по себе задача распараллеливается на 100% — нет никакой взаимосвязности между разными датами на этих этапах.
2. Для этого реперного каркаса притягиваются координаты всех животных и кораблей (одним джойном на весь фрейм)
3. Векторно считаются расстояния между кораблями и животными.
4. Дальше в группе оставляются только Top N. Используем трюк с пакетной сортировкой, примерно такая конструкция

setorder(dt, -r) %>%
    .[, head(.SD, 10), by = .(date, an_id, sh_id)]

где an_id — id животного, sh_id — id корабля, r — расстояние. Это не самый быстрый (быстрее через однократную сборку через .I), но все еще понятный любому пользователю.
Все, получили набор измерений, можно строить распределения по датам.

Для генерации случайных точек тоже применяем специфику векторных вычислений и представлений данных внутри. Из соображений производительности тоже не стоит зацикливаться на циклах. Ключевая подсказка в следующих мантрах:

library(dqrng)
library(tictoc)
library(magrittr)

nn <- 1e4
tic()
ll <- 1:nn %>%
  purrr::map(~dqsample(1:12, 12, replace = FALSE))

df <- matrix(unlist(ll), byrow = TRUE, ncol = 12) %>%
  tibble::as_tibble()
toc()

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

• «R, Монте-Карло и enterprise задачи, часть 2»
• «Ещё немного о параллельных вычислениях в R»

Итого, вся задача сжимается на секунды расчетов. Для больших объемов эффективно применение параллелизации. За день можно провести не одну сотню формирования гипотез, прогона расчетов и проверки результатов. Повышение эффективности работы налицо. Простые джойны и векторизация нам помогли.

Задача 2. Дома и магазины.

Масшатабируем условия. Пришла задачка из ГИС тематики и звучит она примерно следующим образом.

• Есть два датафрейма, в одном («А») объекты с координатами. В другом, («Б»), другие объекты с координатами.
• Один датасет содержит 10^5 точек, другой — 5*10^5 точек. Как сделать функцию, чтобы в для каждого объекта из А узнать количество объектов из датафрейма Б в радиусе N метров?

Решение

Можно было бы воспользоваться уже проверенным выше решением cross join + фильтрация.
Однако, для начала Используем информацию о размерах — один датасет 10^5 точек, другой — 5*10^5 точек.

Из таких показателей сразу проистекает ответ о сложности создания кросс-джойна простыми средствами — это будет 5*10^10 строк — даже превышает емкость Int32 (если используется как индекс), да и две колонки с пересеченными id даже для Int32 в id составит ~400 Gb. И это мы еще даже считать не приступили, только стол очищаем для работы.

fs::fs_bytes(5e10 * 4 * 2)

Казалось бы, вот оно, пространство для применения BigData, хадупов, спарков и прочих обогревателей атмосферы. Засчитаем на 100 нодах и все будет круто!

Если же немного подумать, то задачу можно утрамбовать до размера калькулятора. Было бы желание.

Предположим, что точки у нас евклидова геометрия (точки относительно близко). Тогда задача легко и быстро решается в 2 итерации. На первый взгляд, это тупое усложнение задачи, но именно такое разбиение задачи позволяет решать ее за секунды на «офисном ноутбуке» без привлечения какой-либо BigData.

1. Сначала грубо посечём множество по описанным квадратам. Для каждой точки из Б делаем слияние по условию (conditional join) abs(x - x0) < N & abs(y - y0) < N. Тем самым мы резко сократили множество возможных точек. В силу того, что условное слияние работает для простых операций сравнения, мы и делаем первичное приближение через квадратное множество.
2. Для оставшихся точек делаем тонкую фильтрацию через честное расстояние (x-x0)^2 + (y-y0)^2 <= N^2

Все, задача решена.

library(data.table)
library(magrittr)

# сэмулируем данные
a_dt <- 5e5 %>% 
  {data.table(x_a = runif(., -100, 100), 
              y_a = runif(., -100, 100),
              id_a = 1:.)}

b_dt <- 1e5 %>% 
  {data.table(x_b = runif(., -100, 100), 
              y_b = runif(., -100, 100),
              id_b = 1:.)}

# решение
nn = 5 # устанавливаем расстояние
# 1. грубое слияние
b_dt[, `:=`(x_l = x_b - nn, x_h = x_b + nn, 
            y_l = y_b - nn, y_h = y_b + nn)]
res_dt <- a_dt[b_dt, .(id_a, x.x_a, x.y_a, id_b, x_b, y_b), 
               on = c("x_a >= x_l", "x_a <= x_h", "y_a >= y_l", "y_a <= y_h")]

# получили ~120 млн.
# 2. тонкая фильтрация
res_dt[, r := sqrt((x.x_a - x_b)^2 + (x.y_a - y_b)^2)]
# имеем на выходе ~96 млн
final_dt <- res_dt[r <= nn]

> final_dt[, .N, by = id_b]

          id_b    N
     1:      1  977
     2:      2 1019
     3:      3  942
     4:      4  989
     5:      5  966
    ---            
 99996:  99996  995
 99997:  99997  930
 99998:  99998 1014
 99999:  99999  979
100000: 100000  937

По бенчмаркам:

expression                               median mem_alloc 
<bch:expr>                             <bch:tm> <bch:byt> 
res_dt <- a_dt[b_dt, on = .(x >=  ...]    41.9s    8.47GB 
res_dt[, `:=`(r, sqrt((x_a - x_b) ...]    1.36s    1.82GB

Немного ссылок по возможностям data.table в контексте join-ов:

• «EFFICIENT IN-MEMORY NON-EQUI JOINS using the #rdatatable package, Arun Srinivasan», UseR! 2016
• «The unequalled joy of non-equi joins»
• «non-equi merge in data.table and epidemiology
• «R – Data.Table Rolling Joins»
• «Understanding data.table Rolling Joins»

Тут сложно удержаться от небольшого комментария в сторону «популярного» в DS среде языка. В дополнение к печальным моментам, упомянутым в «R vs Python в продуктивном контуре», приходится отметить что даже для таких тривиальных задачек питон оказывается неудобен. В пандасе нет и, похоже, не предвидится conditional join. Использовать БД или ссылаться на Polars, который имплементирует идеи Apache Arrow. Тут пока сложно сказать, что лучше. Arrow развивается семимильными шагами и является кроссплатформенным и кроссязычным. А если использовать полностью сторонние средства, то вопрос питона плавно исчезает — эти средства можно использовать откуда угодно, хоть из командной строки.

Задача 3. Ищем e-mail ботов

Тоже вполне бытовая задачка. Есть какой-то портал. На нем регистрируются действия пользователей. У пользователей есть электронная почта. Иногда бывают не пользователи, а боты-скрипты. Один из признаков — их почтовые адреса генерируются по шаблону. Хотелось бы кластеризовать и локализовать пачки адресов, которые потенциально генерируются скриптами. Теоретически все просто и понятно. Есть адреса, строим матрицу расстояний и кластеризуем. В качестве метрики расстояния оптимально по смыслу использовать расстояние Левеншнейтна.

Теперь от теории опускаемся к практике. Имеем ~10^7адресов. Матрица 10^7 x 10^7 — вот и BigData подъехала. Финиш?

Подключим голову. Какие могут быть практические соображения?

1. E-mail адреса должны соответствовать определенным стандартам. Возьмем RFC 5321. «Simple Mail Transfer Protocol» и RFC 5322. «Internet Message Format» в качестве отправной точки.
2. Потенциально кластеризуемые адреса злостных скриптов должны принадлежать одному домену. Сложно, глядя только на почтовый адрес, предположить, что адреса admin@yandex.ru и mimimi@aol.com попадут в один кластер.
3. Поскольку мы говорим о расстоянии Левенштейна, то, очевидно, есть какой-то разумный порог похожести. Сложно заподозрить в принадлежности к одному кластеру два таких email 11@ya.ru и ivan.ivanovich.ivanov-1990@Megatool.ipipe.ru.

Из этих соображений получается вполне простой и компактный алгоритм. Называется «разреженные матрицы».

1. Сгрупируем почтовые адреса по доменам и одинаковому количеству символов.
2. Для каждой группы определим потенциальное подмножество для поиска. Например, тот же домен +- 3 символа в длине.

Вот и все. Вселенских размеров матрица развалилась на совокупность крошечных космических частиц. Задачу можно делать в многопоточном режиме, каждый из элементов выступает в качестве самодостаточной сущности. А в конце мы склеим (reduce) все результаты в единое целое. «Офисный» ноутбук вполне успешно пережевывает эту задачку за несколько минут.

Полезные ссылки:

• «The stringdist Package for Approximate String Matching» by Mark P.J. van der Loo
• «Approximate text matching with the stringdist package» by Mark van der Loo. UseR! 2014

Кластеризация e-mail является всего-лишь частным примером, подобные задачи возникают регулярно в том или ином виде. В частности, в биоинформатике есть задача Sequence Alignment, реализованная, например, в пакете bioconductor:msa (Multiple Sequence Alignment). Аналогичная задача существует и в предметной области process mining.

Задача 4. Поиск похожих документов.

Тут как с луковицей, можно снимать слой за слоем… Но не будем здесь погружаться во вселенную поисковых машин. Этой тематике посвящено множество материалов. Просто отмечу, что и эту невообразимую по объемам задачу, утрамбовали до осязаемых размеров путем ряда размышлений и допущений. Речь идет об алгоритме «Locality Sensitive Hashing». Его можно использовать и для своих задачек по матчингу.

Предыдущая публикация — «Лущим веб с помощью R».