test_that("airGR::Calibration should work", {
  ## loading catchment data
  data(L0123001)

  ## preparation of InputsModel object
  InputsModel <- CreateInputsModel(
    RunModel_GR4J,
    DatesR = BasinObs$DatesR,
    Precip = BasinObs$P,
    PotEvap = BasinObs$E
  )

  ## calibration period selection
  Ind_Run <- seq(
    which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
    which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31")
  )
  Ind_WarmUp <- seq(
    which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1989-01-01"),
    which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1989-12-31")
  )

  ## preparation of RunOptions object
  RunOptions <- CreateRunOptions(
    RunModel_GR4J,
    InputsModel = InputsModel,
    IndPeriod_Run = Ind_Run,
    IndPeriod_WarmUp = Ind_WarmUp
  )

  ## calibration criterion: preparation of the InputsCrit object
  InputsCrit <- CreateInputsCrit(
    ErrorCrit_NSE,
    InputsModel = InputsModel,
    RunOptions = RunOptions,
    Obs = BasinObs$Qmm[Ind_Run]
  )

  ## preparation of CalibOptions object
  CalibOptions <- CreateCalibOptions(
    RunModel_GR4J,
    FUN_CALIB = Calibration_Michel
  )

  ## calibration
  OutputsCalib <- Calibration(
    InputsModel = InputsModel,
    RunOptions = RunOptions,
    InputsCrit = InputsCrit,
    CalibOptions = CalibOptions,
    FUN_MOD = RunModel_GR4J,
    FUN_CALIB = Calibration_Michel
  )

  expect_length(OutputsCalib$ParamFinalR, 4)
})

# data set up
e <- runCalibration(runRunModel = TRUE, FUN_CRIT = ErrorCrit_NSE)
for (x in ls(e)) {
  assign(x, get(x, e))
}

test_that("Calibrated parameters equals reference parameters without regularization", {
  InputsCrit <- CreateInputsCrit(
    InputsModel,
    FUN_CRIT = ErrorCrit_NSE,
    RunOptions = RunOptions,
    Obs = Qobs[
      IndPeriod_Run,
      np$id[np$calibration == "Gauged"],
      drop = FALSE
    ],
    AprioriIds = NULL
  )
  e <- runCalibration(runRunModel = TRUE, InputsCrit = InputsCrit)
  for (x in ls(e)) {
    assign(x, get(x, e))
  }
  lapply(names(Param), function(id) {
    expect_equal(Param[[!!id]], ParamMichel[[id]])
  })
})

skip_on_cran()

test_that("Calibration with regularization is OK", {
  InputsCrit <- CreateInputsCrit(
    InputsModel = InputsModel,
    RunOptions = RunOptions,
    Obs = Qobs[IndPeriod_Run, ],
    AprioriIds = c(
      "54057" = "54032",
      "54032" = "54001",
      "54001" = "54095"
    ),
    transfo = "sqrt"
  )

  OC <- Calibration(
    InputsModel = InputsModel,
    RunOptions = RunOptions,
    InputsCrit = InputsCrit,
    CalibOptions = CalibOptions
  )

  ParamLavenne <- extractParam(OC)
  expect_equal(
    OC[["54095"]]$CritFinal,
    ErrorCrit(
      InputsCrit[["54095"]],
      RunModel(InputsModel, RunOptions, ParamLavenne)[["54095"]]
    )$CritValue
  )
  OM <- RunModel(InputsModel, RunOptions, ParamLavenne)
  lapply(names(OC), function(id) {
    expect_gt(
      ErrorCrit(
        InputsCrit[[id]],
        OM[[id]]
      )$CritValue,
      0.89
    )
  })
})

test_that("Calibration with Diversion works", {
  n_div <- rbind(
    nodes,
    data.frame(
      id = "54029",
      down = "54002",
      length = 50,
      area = NA,
      model = "Diversion"
    )
  )
  g_div <- CreateGRiwrm(n_div)
  Qmin = matrix(1E5, nrow = length(DatesR), ncol = 1)
  colnames(Qmin) = "54029"
  Qdiv <- -Qmin
  IM_div <- CreateInputsModel(
    g_div,
    DatesR,
    Precip,
    PotEvap,
    Qinf = Qdiv,
    Qmin = Qmin
  )
  RO_div <- setupRunOptions(IM_div)$RunOptions
  P_div <- ParamMichel
  P_div$`54002` <- c(1, ParamMichel$`54002`)
  IC_div <- CreateInputsCrit(
    InputsModel = IM_div,
    RunOptions = RO_div,
    Obs = Qobs[IndPeriod_Run, ],
  )
  CO_div <- CreateCalibOptions(IM_div)
  OC <- Calibration(
    InputsModel = IM_div,
    RunOptions = RO_div,
    InputsCrit = IC_div,
    CalibOptions = CO_div
  )
  expect_length(OC$`54002`$ParamFinalR, 5)
})

test_that("Derivation and normal connection should return same calibration", {
  n_2ol <- nodes[nodes$id %in% c("54095", "54001"), ]
  n_2ol[n_2ol$id %in% c("54095", "54001"), c("down", "length")] <- c(NA, NA)
  meteoIds <- n_2ol$id
  n_2ol$area[n_2ol$id == "54001"] <-
    n_2ol$area[n_2ol$id == "54001"] - n_2ol$area[n_2ol$id == "54095"]
  n_2ol <- rbind(
    n_2ol,
    data.frame(
      id = "54095",
      down = "54001",
      length = 42,
      area = NA,
      model = "Diversion"
    ),
    data.frame(
      id = "upstream",
      down = "54095",
      length = 0,
      area = NA,
      model = NA
    )
  )
  g_2ol <- CreateGRiwrm(n_2ol)

  # Add upstream flow on 54095 that is removed by the Diversion
  # and derive previously simulated flow in order to get the same Qsim as before
  Qinf = matrix(0, nrow = length(DatesR), ncol = 2)
  Qinf[IndPeriod_Run, 1] <- OM_GriwrmInputs[["54095"]]$Qsim_m3
  Qinf[IndPeriod_Run, 2] <- -OM_GriwrmInputs[["54095"]]$Qsim_m3
  Qinf[IndPeriod_WarmUp, 1] <- OM_GriwrmInputs[[
    "54095"
  ]]$RunOptions$WarmUpQsim_m3
  Qinf[IndPeriod_WarmUp, 2] <- -OM_GriwrmInputs[[
    "54095"
  ]]$RunOptions$WarmUpQsim_m3

  colnames(Qinf) <- c("upstream", "54095")

  Qmin = matrix(0, nrow = length(DatesR), ncol = 1)
  colnames(Qmin) <- "54095"

  IM_2ol <- CreateInputsModel(
    g_2ol,
    DatesR,
    Precip[, meteoIds],
    PotEvap[, meteoIds],
    Qinf = Qinf,
    Qmin = Qmin
  )

  # Copy area of upstream node to downstream node in order to get
  # correct conversion of Qsim in mm
  IM_2ol[["54001"]]$BasinAreas[1] <- tail(IM_2ol[["54095"]]$BasinAreas, 1)

  RO_2ol <- setupRunOptions(IM_2ol)$RunOptions
  IC_2ol <- CreateInputsCrit(
    InputsModel = IM_2ol,
    FUN_CRIT = ErrorCrit_NSE,
    RunOptions = RO_2ol,
    Obs = Qobs[IndPeriod_Run, ]
  )
  CO_2ol <- CreateCalibOptions(IM_2ol)
  CO_2ol[["54095"]]$FixedParam[1] <- 1
  OC_2ol <- Calibration(
    InputsModel = IM_2ol,
    RunOptions = RO_2ol,
    InputsCrit = IC_2ol,
    CalibOptions = CO_2ol
  )
  ParamRef <- Param[names(IM_2ol)]
  ParamRef[["54095"]] <- c(1, ParamRef[["54095"]])
  ParamFinalR <- extractParam(OC_2ol)
  lapply(names(ParamFinalR), function(id) {
    expect_equal(
      OC_2ol[[id]]$CritFinal,
      OutputsCalib[[id]]$CritFinal,
      tolerance = 1E-5
    )
  })
  #Excepted parameter #2 of GR4J all others are equal (precision 3/1000)
  lapply(names(ParamFinalR), function(id) {
    expect_equal(
      ParamFinalR[[!!id]][-3] / ParamRef[[!!id]][-3],
      rep(1, 4),
      tolerance = 3E-2
    )
  })
})