如何基于R tidyverse计算离散时序曲线的山谷区域面积

计算逻辑说明 #

和需求匹配的计算规则如下：

• 首行数据、当日y值小于等于前一日y值（下降或持平）时，当日a值直接为0
• 当日y值大于前一日时：
  1. 找到当前行左侧最近的、y值大于等于当日y值的行位置j
  2. 取j和当前行中间的所有日期区间，计算当日y值乘以区间天数，减去区间内所有y值的总和，得到的结果就是当日的a值
     该逻辑完全匹配给出的测试数据预期结果，计算结果的单位也符合回淹规模m²·天的业务定义。

基于tidyverse的实现代码 #

# 加载依赖包
library(tidyverse)
library(lubridate)

# 构造测试数据
dat <- tibble(
  d = seq.Date(ymd("2021-01-01"), ymd("2021-01-10"), by = "1 day"),
  y = c(100,95,90,70,50,55,60,75,85,90),
  expected_a = c(0,0,0,0,0,5,15,65,115,145)
)

# 线性时间复杂度的计算函数，适配18万行级别的大数据量
calc_reflood_area <- function(y) {
  n <- length(y)
  a <- numeric(n)
  # 单调栈存储左侧递减的y值对应行号和数值，保证查找最近符合要求的j的效率
  stack <- list()
  
  for (i in 1:n) {
    # 下降/持平/首行直接返回0
    if (i == 1 || y[i] <= y[i-1]) {
      a[i] <- 0
    } else {
      # 弹出栈内所有小于当前y值的元素
      while (length(stack) > 0 && stack[[length(stack)]][2] < y[i]) {
        stack <- stack[-length(stack)]
      }
      j <- stack[[length(stack)]][1]
      interval_sum <- sum(y[(j+1):(i-1)])
      a[i] <- y[i] * (i - j - 1) - interval_sum
    }
    # 维护单调栈的递减特性
    while (length(stack) > 0 && stack[[length(stack)]][2] <= y[i]) {
      stack <- stack[-length(stack)]
    }
    stack[[length(stack) + 1]] <- c(i, y[i])
  }
  return(a)
}

# 执行计算
dat <- dat %>%
  mutate(a = calc_reflood_area(y))

运行后得到的a列和测试数据给出的预期值完全一致。

内容的提问来源于stack exchange，提问作者fantanaman