# 各輸入檔案乃採用原發問者上傳之 csv 和 nc 檔

library(magrittr)
library(data.table)
library(ncdf4)

# 各種讀資料
analog.indices <-
  fread("HIRAM_WRF_analog.indices(1979-2005).csv",
        col.names = c("day", "i1", "i2")) %>%
  .[1:9862, -1] %>% # 發問者要求刪除 9862 列之後；第 1 欄刪除
  as.matrix # 轉 matrix 之後較好懂
analog.weights <-
  fread("HIRAM_WRF_weights(1979-2005).csv",
        col.names = c("day", "w1", "w2")) %>%
  .[1:9862, -1] %>% # 發問者要求刪除 9862 列之後；第 1 欄刪除
  as.matrix # 轉 matrix 之後較好懂
HIRAM_WRF_data <- nc_open("197901-200512_pr_axis_time_domain.nc")
hiram_wrf_lon <- ncvar_get(HIRAM_WRF_data, "lon")
hiram_wrf_lat <- ncvar_get(HIRAM_WRF_data, "lat")
hiram_wrf_time <- ncvar_get(HIRAM_WRF_data, "time")
hiram_wrf_pr <- ncvar_get(HIRAM_WRF_data, "pr")
dim(hiram_wrf_pr)

pr.new <-
  array(
    NA_real_,
    dim = c(dim(hiram_wrf_pr)[1:2], nrow(analog.indices)),
    # 填上合適的 dimnames（尚無實際作用）
    dimnames = list(
      hiram_wrf_lon %>% sprintf("%1.2f", .),
      hiram_wrf_lat %>% sprintf("%1.2f", .),
      1:nrow(analog.indices)
    )
  )

##########################################
# 方法 A：直接列式子
# analog.indices 逐列處理（1:10227）
# 每一圈會為 pr.new 填上一層資料
pr.new.A <- pr.new
for (i in 1:nrow(analog.indices)) {
  pr.new.A[, , i] <-
    hiram_wrf_pr[, , analog.indices[i, 1]] * analog.weights[i, 1] +
    hiram_wrf_pr[, , analog.indices[i, 2]] * analog.weights[i, 2]
}

##########################################
# 方法 B：二層巢式迴圈
# 外層 i 照舊走直的，內層 j 走橫的
# 每一小圈會計算部份累積的總和
# 每一大圈會為 pr.new 填上一層資料
pr.new.B <- pr.new
for (i in 1:nrow(analog.indices)) {
  partialSum <- 0.0
  for (j in 1:ncol(analog.indices)) {
    partialSum <-
      partialSum +
      hiram_wrf_pr[, , analog.indices[i, j]] * analog.weights[i, j]
  }
  pr.new.B[, , i] <- partialSum
  # 以上用了一個以迴圈求累積加總的常見小技巧
  # 沒看過的話，看看這個求 4+6+7 的例子：
  # a <- c(4,6,7); b <- 0; for(i in 1:3){b <- b + a[i]}; print(b)
}

###################################################
# 檢查第11日第 (1,1) 格的結果
analog.indices[11,]
# i1 i2 
# 10 11 
analog.weights[11,]
#        w1        w2 
# 0.2276517 0.7622840 
print(hiram_wrf_pr[1, 1, 10] * 0.2276517 + hiram_wrf_pr[1, 1, 11] * 0.762284) # 0.007988894
print(pr.new.A[1, 1, 11]) # 0.007988894
print(pr.new.B[1, 1, 11]) # 0.007988894
identical(pr.new.A, pr.new.B) # 二者全等