hydro-topmodel 0.1.0

#!/usr/bin/env Rscript
# 生成 TOPMODEL bit-exact oracle:跑 R `topmodel` 包 on huagrahuma 数据集 →
# dump 6 过程序列(Q/qo/qs/S/fex/Ea)+ 输入(parameters/topidx/delay/rain/ETp)→ JSON。
#
# 用法(从 hydro-topmodel/ crate 根运行):
#   Rscript validation/topmodel/gen_reference.R
# 依赖:R + topmodel 包(从 CRAN archive 装:
#   install.packages("https://cran.r-project.org/src/contrib/Archive/topmodel/topmodel_0.7.5.tar.gz", repos=NULL, type="source")
#   )+ jsonlite。
# 输出:tests/topmodel_data/oracle.json(供 tests/topmodel_crosscheck.rs 比对)。

library(topmodel)
library(jsonlite)

data(huagrahuma)

cat("=== huagrahuma 结构(确认变量名)===\n")
print(ls.str(huagrahuma))

# 标准变量名(若 ls.str 显示不同,按实际调整)
parameters <- huagrahuma$parameters   # 11 参 c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,K0,CD,dt)
topidx     <- huagrahuma$topidx       # 2 列矩阵 (atb, Aatb_r)
delay      <- huagrahuma$delay        # 2 列矩阵 (d, Ad_r)
rain       <- huagrahuma$rain
ETp        <- huagrahuma$ETp

# 跑模型 verbose → 6 序列(Q/qo/qs/S/fex/Ea)
res <- topmodel(parameters, topidx, delay, rain, ETp, verbose = TRUE)

oracle <- list(
  parameters     = as.numeric(parameters),
  topidx_atb     = as.numeric(topidx[, 1]),
  topidx_Aatb_r  = as.numeric(topidx[, 2]),
  delay_d        = as.numeric(delay[, 1]),
  delay_Ad_r     = as.numeric(delay[, 2]),
  rain           = as.numeric(rain),
  ETp            = as.numeric(ETp),
  Q              = as.numeric(res$Q),
  qo             = as.numeric(res$qo),
  qs             = as.numeric(res$qs),
  S              = as.numeric(res$S),
  fex            = as.numeric(res$fex),
  Ea             = as.numeric(res$Ea)
)

out <- "tests/topmodel_data/oracle.json"
dir.create(dirname(out), recursive = TRUE, showWarnings = FALSE)
writeLines(toJSON(oracle, auto_unbox = TRUE, digits = 17), out)
cat("\n写入", out, "\n")
cat("序列长度:", length(rain), " 参数数:", length(parameters), "\n")
cat("NSE(R 自评, Qsim vs Qobs):", topmodel(parameters, topidx, delay, rain, ETp, Qobs = if (is.null(huagrahuma$Qobs)) NA else huagrahuma$Qobs), "\n")