The OP has requested to find the quickest way to recode ranges of day of year to month and count and has mentioned that the production dataset includes 10 M rows. The OP has conducted a benchmark test with a problem size of 2000 rows and test data which cover only 3 months instead of 12 months.
Although the question has an accepted answer I wondered
- how a data.table approach using
melt(), foverlaps(), and dcast() would compare to the other answers
- and how a more realistic benchmark with varying problem sizes would look like.
foverlaps()
library(data.table)
library(magrittr)
cols <- c("deb", "fin")
# reshape sub ranges from wide to long format
long <- melt(setDT(d)[, rn := .I], id.vars = "rn", measure.vars = patterns(cols),
value.name = cols)[deb > 0]
# create data.table with monthly breaks and set keys
b <- seq(as.IDate("2018-01-01"), as.IDate("2019-01-01"), "month")
mDT <- data.table(abb = forcats::fct_inorder(month.abb),
deb = yday(head(b, 12L)),
fin = yday(tail(b, 12L) - 1L),
key = c("deb", "fin"))
# find overlaps between sub ranges and monthly breaks
foverlaps(long, mDT)[
# compute days in overlaps
, days := pmin(fin, i.fin) - pmax(deb, i.deb) + 1L] %>%
# reshape to wide format for final result
dcast(rn ~ abb, sum, value.var = "days", fill = 0L, drop = FALSE)
rn Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1: 1 31 10 1 0 0 0 0 0 0 0 0 0
2: 2 0 28 0 0 0 0 0 0 0 0 0 0
3: 3 0 0 1 1 0 0 0 0 0 0 0 0
4: 4 0 0 0 1 0 0 0 0 0 0 0 0
Benchmark
The benchmark code below compares
AntoniosK's tidyverse answer (now deleted) has not been considered because the OP had commented that this approach is 30 times slower than his own g2 approach.
For the benchmark, the problem size is varied (n_rows = 100, 1000, 10000) as well as the share of rows with two ranges (p_sr_2 = 10%, 50%, 100%). The maximum problem size has been choosen with respect to the long run times of some of the approaches.
The different approaches differ also in the class of the result (i.e., data.table, matrix, tibble, and table) and the number of columns. Therefore, emphasis was put on writing an appropriate check function.
library(bench)
library(data.table)
library(magrittr)
library(ggplot2)
bm2 <- press(
p_sr_2 = c(0.1, 0.5, 1), # share or rows with 2nd sub range ]0, 1]
n_rows = 10^seq(2, 4, 1),
{ # create test data
set.seed(1L)
d0 <- t(replicate(n_rows, sample(365L, 4L) %>% sort())) %>% data.table()
setnames(d0, c("deb", "fin", "deb2", "fin2"))
idx <- sample(nrow(d0), (1 - p_sr_2) * nrow(d0))
d0[idx, c("deb2", "fin2") := 0L]
str(d0)
mark(
foverlaps = {
d <- copy(d0)
cols <- c("deb", "fin")
long <- melt(setDT(d)[, rn := .I], id.vars = "rn", measure.vars = patterns(cols),
value.name = cols)[deb > 0]
b <- seq(as.IDate("2018-01-01"), as.IDate("2019-01-01"), "month")
mDT <- data.table(abb = forcats::fct_inorder(month.abb),
deb = yday(head(b, 12L)),
fin = yday(tail(b, 12L) - 1L),
key = c("deb", "fin"))
foverlaps(long, mDT)[, days := pmin(fin, i.fin) - pmax(deb, i.deb) + 1L] %>%
dcast(rn ~ abb, sum, value.var = "days", fill = 0L, drop = FALSE)
# returns a data.table with 13 columns and 0 for missing values
},
g1 = {
f1 <- function(deb, fin, deb2, fin2, ...) {
f <- factor(c(deb:fin, deb2:fin2))
levels(f) <- list(jan = 1:31,
feb = 32:59,
mar = 60:90,
apr = 91:120,
may = 121:151,
jun = 152:181,
jul = 182:212,
aug = 213:243,
sep = 244:273,
oct = 274:304,
nov = 305:334,
dec = 335:365)
table(f)
}
do.call(rbind, d %>% purrr::pmap(f1))
# returns a matrix with 12 named columns and 0 for missing values
},
g2 = {
K <- vector(10, mode = "character")
K[1:31] <- "jan"
K[32:59] <- "feb"
K[60:90] <- "mar"
K[91:120] <- "apr"
K[121:151] <- "may"
K[152:181] <- "jun"
K[182:212] <- "jul"
K[213:243] <- "aug"
K[244:273] <- "sep"
K[274:304] <- "oct"
K[305:334] <- "nov"
K[335:365] <- "dec"
f2 <-
Vectorize(function(deb, fin, deb2, fin2)
table(c(K[deb:fin], K[deb2:fin2])))
template <- matrix(
integer(0), ncol = 12L,
dimnames = list(NULL, c("jan", "feb", "mar", "apr", "may", "jun",
"jul", "aug", "sep", "oct", "nov", "dec"))) %>%
tibble::as.tibble()
do.call(dplyr::bind_rows, c(list(template), f2(d$deb, d$fin, d$deb2, d$fin2)))
# returns a tibble with 12 columns and NA for missing values
},
g3 = {
K <- vector(10, mode = "character")
K[1:31] <- "jan"
K[32:59] <- "feb"
K[60:90] <- "mar"
K[91:120] <- "apr"
K[121:151] <- "may"
K[152:181] <- "jun"
K[182:212] <- "jul"
K[213:243] <- "aug"
K[244:273] <- "sep"
K[274:304] <- "oct"
K[305:334] <- "nov"
K[335:365] <- "dec"
names(K) <- 1:365
f3 <-
Vectorize(function(deb, fin, deb2, fin2)
c(K[deb:fin], K[deb2:fin2]))
m <- as.matrix(do.call(dplyr::bind_rows, f3(d$deb, d$fin, d$deb2, d$fin2)))
z <- unlist(purrr::map(list("jan", "feb", "mar", "apr", "may", "jun",
"jul", "aug", "sep", "oct", "nov", "dec"),
function(y)
apply(m, 1, function(x)
sum(x == y, na.rm = TRUE))))
dim(z) <- c(nrow(d), 12)
z
# returns a matrix with 12 columns and 0 for missing values
},
henrik = {
d <- copy(d0)
b <- as.numeric(format(seq(as.Date("2018-01-01"), as.Date("2018-12-31"),
by = "month"), "%j"))
l <- Map(
function(from, to, from2, to2) month.abb[findInterval(c(from:to, from2:to2), b)],
d$deb, d$fin, d$deb2, d$fin2)
i <- rep(1:nrow(d), lengths(l))
table(i, factor(unlist(l), levels = month.abb))
# returns an object of class table with 12 columns and 0 for missing values
},
chinsoon12 = {
d <- copy(d0)
firstOfMths <- seq(as.Date("2018-01-01"), as.Date("2019-01-01"), by="month")
daysPerMth <- c(1L, cumsum(as.integer(diff(firstOfMths))))
g <- ceiling(seq(1, ncol(d)) / 2)
t(apply(d, 1, function(x) {
x <- unlist(by(x, g, function(k) seq(k[1L], k[2L])), use.names=FALSE)
table(cut(x, daysPerMth, labels=month.abb, include.lowest=TRUE, right=TRUE))
}))
# returns a matrix with 12 named columns and 0 for missing values
},
check = function(x, y) {
cat("Run check: ")
xdt <- as.data.table(x) %>% .[, .SD, .SDcols = tail(seq_len(ncol(.)), 12L)]
if (tibble::is_tibble(y)) {
y <- dplyr::mutate_all(y, function(x) dplyr::if_else(is.na(x), 0L, x))
}
if (is.table(y)) y <- matrix(y, ncol = 12L)
ydt <- as.data.table(y) %>% .[, .SD, .SDcols = tail(seq_len(ncol(.)), 12L)]
result <- all.equal(xdt, ydt, check.attributes = FALSE)
if (!isTRUE(result)) {
print(result)
} else cat("OK\n")
return(result)
}
)
}
)
The timings of the benchmarks can be plotted :
library(ggplot2)
autoplot(bm)
Please, note the logarithmic time scale.
Apparently, the share of rows with two ranges doesn't have a visible effect on performance in contrast to the number of rows. For smaller problem sizes, henrik's approach is the fastest approach confirming OP's observations. However, for problem sizes of 1000 rows and more the foverlaps approach is significantly faster. For 10 k rows, foverlaps is about one to two magnitudes faster than the other approaches.
Also, the memory requirements differ strongly:
bm %>%
tidyr::unnest() %>%
ggplot(aes(expression, mem_alloc, color = gc)) +
ggbeeswarm::geom_quasirandom() +
coord_flip() +
facet_grid(p_sr_2 ~ n_rows, labeller = label_both)
Again, please note the log scale.
The foverlaps approach allocates about one to two magnitudes less memory than the other approaches.
Benchmark part 2
Due to the long run times (and my impatience), above benchmark was only conducted for up to 10 k rows. Also Henrik tested only 40 k rows. So, I wondered if the foverlaps approach would be capable to process 10 M rows (the size of OP's production dataset) in reasonable time.
Unfortunately, my code to create the input data turned out to be too slow for problem sizes of 1 M rows or greater. So, I had to optimize (and benchmark) this part separately.
Only the foverlaps approach is benchmarked with the share of a second range fixed to 25% as stated by the OP.
library(bench)
library(data.table)
library(magrittr)
library(ggplot2)
bm5 <- press(
n_rows = 10^seq(2, 7, 1),
{ # create test data
cat("Start to create test data", format(Sys.time()), "\n")
p_sr_2 <- 0.25 # share or rows with 2nd sub range ]0, 1]
set.seed(1L)
long <- data.table(rn = rep(seq_len(n_rows), each = 4L),
day = sample(365L, 4L * n_rows, replace = TRUE))
setorder(long, rn, day)
dups <- long[, which(anyDuplicated(day) > 0), by = rn]$rn
if (length(dups) > 0) long[
rn %in% dups,
day := replicate(length(dups), sample(365L, 4L) %>% sort(), simplify = FALSE) %>% unlist()]
d0 <- dcast(long, rn ~ rowid(rn), value.var = "day")[, rn := NULL]
setnames(d0, c("deb", "fin", "deb2", "fin2"))
idx <- sample(nrow(d0), (1 - p_sr_2) * nrow(d0))
d0[idx, c("deb2", "fin2") := 0L]
str(d0)
rm(long) # free the memory
tables()
cat("Test data created", format(Sys.time()), "\n")
mark(
foverlaps = {
d <- copy(d0)
cols <- c("deb", "fin")
long <- melt(setDT(d)[, rn := .I], id.vars = "rn", measure.vars = patterns(cols),
value.name = cols)[deb > 0]
b <- seq(as.IDate("2018-01-01"), as.IDate("2019-01-01"), "month")
mDT <- data.table(abb = forcats::fct_inorder(month.abb),
deb = yday(head(b, 12L)),
fin = yday(tail(b, 12L) - 1L),
key = c("deb", "fin"))
foverlaps(long, mDT)[, days := pmin(fin, i.fin) - pmax(deb, i.deb) + 1L] %>%
dcast(rn ~ abb, sum, value.var = "days", fill = 0L, drop = FALSE)
# returns a data.table with 13 columns and 0 for missing values
},
min_time = 2
)
}
)
The run time of the 10 M rows case was about 23.7 seconds on my system. The run time is increasing nearly linear for problem sizes above 1000 rows. The slight upward bend for 10 M rows might be due to memory limitations on my system.
bm4 %>%
tidyr::unnest() %>%
ggplot(aes(n_rows, time, colour = expression)) +
geom_smooth(data = . %>% dplyr::filter(n_rows > 10^3),
method = "lm", se = FALSE, colour = "blue", size = 1) +
geom_point(size = 2) +
scale_x_log10() +
stat_summary(fun.y = median, geom = "line", size = 1) +
ggtitle("time vs n_rows")
Please, note the double logarithmic scale.
