#eq 1
eq.Ms <- function(Ms, Gs, Ls){
  Ms + Gs - Ls
}

#eq 2
eq.Mr <- function(Mr, Gr, Lr){
  Mr + Gr - Lr
}

#eq 3
eq.Gs <- function(g, Cs, Ns, Ms){
  g * Cs * Ns / Ms
}

#eq 4
eq.Gr <- function(g, Cr, Nr, Mr){
  g * Cr * Nr / Mr
}

#eq 5
eq.Ls <- function(m, Ms, Km){
  m * Ms / (1 + Km / Ms)
}

#eq 6
eq.Lr <- function(m, Mr, Km){
  m * Mr / (1 + Km / Mr)
}

#eq 7
eq.Uc <- function(a, Ms, Ka, Cs, Jc){
  a * Ms / ((1 + Ms/Ka)*(1 + Cs/(Ms * Jc)))
}

#eq 8
eq.Un <- function(b, Mr, Ka, Nr, Jn){
  b * Mr / ((1 + Mr/Ka)*(1 + Nr/(Mr*Jn)))
}

#eq 9
eq.Tauc <- function(Ms, Mr, Cs, Cr){
  (Ms * Mr / (Ms + Mr)) * (Cs / Ms - Cr / Mr)
}

#eq 10
eq.Taun <- function(Ms, Mr, Ns, Nr){
  (Ms * Mr / (Ms + Mr)) * (Nr / Mr - Ns / Ms)
}

#eq 11
eq.Cs <- function(Cs, Uc, fc, Gs, Tauc){
  Cs + Uc - fc * Gs - Tauc
}

#eq 12
eq.Cr <- function(Cr, Tauc, fc, Gr){
  Cr + Tauc - fc * Gr
}

#eq 13
eq.Ns <- function(Ns, Taun, fn, Gs){
  Ns + Taun - fn * Gs
}

#eq 14
eq.Nr <- function(Nr, Un, fn, Gr, Taun){
  Nr + Un - fn * Gr - Taun
}

#step function
S <-function(X, b1, b2){
  max(
    min(
      (X - b1) / (b2 - b1),
      1
    ),
    0
  )
}
#trapezoid function
Z <- function(X, b1, b2, b3, b4){
  max(
    min(
      (X - b1)/(b2 - b1),
      1,
      (b4 - X)/(b4 - b3)
    ),
    0
  )
}

#eq 15
a.std.eq <- function(amax, Z_T, S_Q, S_M, S_NsMs){
  amax * min(Z_T, S_Q, S_M, S_NsMs)
}

#eq 16
a.fqr.eq <- function(afqr, S_M, S_NsMs){
  afqr * min(S_M, S_NsMs)
}

#eq 17
a.red.eq <- function(afqr, S_M, S_NsMs){
  afqr * S_M * S_NsMs
}

#eq 18
a.oak.eq <- a.red.eq

#eq 19
b.std.eq <- function(bmax, S_T, Z_M, S_Nsoil){
  bmax * min(S_T, Z_M, S_Nsoil)
}
b.fqr.eq <- b.std.eq

#eq 20
b.red.eq <- function(bmax, S_T, Z_M, S_Nsoil){
  bmax * S_T * Z_M  * S_Nsoil
}
b.oak.eq <- b.red.eq

#eq 21
g.std.eq <- function(gmax, Z_T){
  gmax * Z_T
}
g.fqr.eq <- g.std.eq

#eq 22
g.red.eq <- function(gmax, Z_T, S_M){
  gmax * Z_T * S_M
}
g.oak.eq <- g.red.eq

#eq 23
#m.std.eq <- function(mmin, mmax, S_T){
#  max(mmin, mmax * S_T)
#}
#m.fqr.eq <- m.std.eq
#m.red.eq <- m.std.eq
#changed in commit 73640cf: mortality with L_mix
m.std.eq <- function(p_1, p_2, mmax, S_T){
  p_1 * mmax + p_2 * mmax * S_T
}
m.fqr.eq <- m.std.eq
m.red.eq <- m.std.eq

#eq 24
# m.oak.eq <- function(mmin, mmax, S_M, S_T){
#   max(mmin, mmax * (1 - S_M), mmax * (1 - S_T))
# }
#changed in commit 73640cf: mortality with L_mix
m.oak.eq <- function(p_1, p_2, p_3, mmax, S_M, S_T){
  p_1 * mmax + p_2 * mmax * S_T + p_3 * mmax * S_M
}


#addendum
positive <- function(v){
  v[which(v < 0)] <- 0.0
  v
}

a.eq <- function(version, state, parameters, constants, forcing_data){
  amax <- constants["CUmax"]
  Ns <- state["Ns"]
  Ms <- state["Ms"]
  switch(
    version,
    std = a.std.eq(
      amax,
      Z(forcing_data["tcur"],
        parameters["CU_tcur_1"],
        parameters["CU_tcur_2"],
        parameters["CU_tcur_3"],
        parameters["CU_tcur_4"]),
      S(forcing_data["ppfd"],
        parameters["CU_ppfd_1"],
        parameters["CU_ppfd_2"]),
      S(forcing_data["swc"],
        parameters["CU_swc_1"],
        parameters["CU_swc_2"]),
      S(Ns / Ms,
        parameters["CU_Ns_1"],
        parameters["CU_Ns_2"])
    ),
    fqr = a.fqr.eq(
      amax / 30 * forcing_data["afqr"],
      S(forcing_data["swc"],
        parameters["CU_swc_1"],
        parameters["CU_swc_2"]),
      S(Ns/Ms,
        parameters["CU_Ns_1"],
        parameters["CU_Ns_2"])
    ),
    red = a.red.eq(
      amax / 30 * forcing_data["afqr"],
      S(forcing_data["swc"],
        parameters["CU_swc_1"],
        parameters["CU_swc_2"]),
      S(Ns/Ms,
        parameters["CU_Ns_1"],
        parameters["CU_Ns_2"])
    ),
    oak = a.oak.eq(
      amax / 30 * forcing_data["afqr"],
      S(forcing_data["swc"],
        parameters["CU_swc_1"],
        parameters["CU_swc_2"]),
      S(Ns/Ms,
        parameters["CU_Ns_1"],
        parameters["CU_Ns_2"])
    )
  )
}

b.eq <- function(version, state, parameters, constants, forcing_data){
  bmax <- constants["NUmax"]
  switch(
    version,
    std = b.std.eq(
      bmax,
      S(forcing_data["tnur"],
        parameters["NU_tnur_1"],
        parameters["NU_tnur_2"]),
      Z(forcing_data["swc"],
        parameters["NU_swc_1"],
        parameters["NU_swc_2"],
        parameters["NU_swc_3"],
        parameters["NU_swc_4"]),
      S(forcing_data["nsoil"],
        parameters["NU_Nsoil_1"],
        parameters["NU_Nsoil_2"])
    ),
    fqr = b.fqr.eq(
      bmax,
      S(forcing_data["tnur"],
        parameters["NU_tnur_1"],
        parameters["NU_tnur_2"]),
      Z(forcing_data["swc"],
        parameters["NU_swc_1"],
        parameters["NU_swc_2"],
        parameters["NU_swc_3"],
        parameters["NU_swc_4"]),
      S(forcing_data["nsoil"],
        parameters["NU_Nsoil_1"],
        parameters["NU_Nsoil_2"])
    ),
    red = b.red.eq(
      bmax,
      S(forcing_data["tnur"],
        parameters["NU_tnur_1"],
        parameters["NU_tnur_2"]),
      Z(forcing_data["swc"],
        parameters["NU_swc_1"],
        parameters["NU_swc_2"],
        parameters["NU_swc_3"],
        parameters["NU_swc_4"]),
      S(forcing_data["nsoil"],
        parameters["NU_Nsoil_1"],
        parameters["NU_Nsoil_2"])
    ),
    oak = b.oak.eq(
      bmax,
      S(forcing_data["tnur"],
        parameters["NU_tnur_1"],
        parameters["NU_tnur_2"]),
      Z(forcing_data["swc"],
        parameters["NU_swc_1"],
        parameters["NU_swc_2"],
        parameters["NU_swc_3"],
        parameters["NU_swc_4"]),
      S(forcing_data["nsoil"],
        parameters["NU_Nsoil_1"],
        parameters["NU_Nsoil_2"])
    )
  )
}

g.eq <- function(version, state, parameters, constants, forcing_data){
  gmax <- constants["gmax"]
  switch(
    version,
    std = g.std.eq(
      gmax,
      Z(forcing_data["tgrowth"],
        parameters["g_tgrowth_1"],
        parameters["g_tgrowth_2"],
        parameters["g_tgrowth_3"],
        parameters["g_tgrowth_4"])
    ),
    fqr = g.fqr.eq(
      gmax,
      Z(forcing_data["tgrowth"],
        parameters["g_tgrowth_1"],
        parameters["g_tgrowth_2"],
        parameters["g_tgrowth_3"],
        parameters["g_tgrowth_4"])
    ),
    red = g.red.eq(
      gmax,
      Z(forcing_data["tgrowth"],
        parameters["g_tgrowth_1"],
        parameters["g_tgrowth_2"],
        parameters["g_tgrowth_3"],
        parameters["g_tgrowth_4"]),
      S(forcing_data["swc"],
        parameters["g_swc_1"],
        parameters["g_swc_2"])
    ),
    oak = g.oak.eq(
      gmax,
      Z(forcing_data["tgrowth"],
        parameters["g_tgrowth_1"],
        parameters["g_tgrowth_2"],
        parameters["g_tgrowth_3"],
        parameters["g_tgrowth_4"]),
      S(forcing_data["swc"],
        parameters["g_swc_1"],
        parameters["g_swc_2"])
    )
  )
}

m.eq <- function(version, state, parameters, constants, forcing_data){
  mmax <- constants["mmax"]
  if(version == "oak") {
    #m3 <- max(0.0, 100.0 - (parameters["L_mix_1"]+parameters["L_mix_2"]))
    sum_m <- parameters["L_mix_1"] + parameters["L_mix_2"] + parameters["L_mix_3"]
  }
  switch(
    version,
    std = m.std.eq(
      parameters["L_mix_1"]/100.0,
      1.0 - parameters["L_mix_1"]/100.0,
      mmax,
      S(forcing_data["tloss"],
        parameters["r_tloss_1"],
        parameters["r_tloss_2"])
    ),
    fqr = m.fqr.eq(
      parameters["L_mix_1"]/100.0,
      1.0 - parameters["L_mix_1"]/100.0,
      mmax,
      S(forcing_data["tloss"],
        parameters["r_tloss_1"],
        parameters["r_tloss_2"])
    ),
    red = m.red.eq(
      parameters["L_mix_1"]/100.0,
      1.0 - parameters["L_mix_1"]/100.0,
      mmax,
      S(forcing_data["tloss"],
        parameters["r_tloss_1"],
        parameters["r_tloss_2"])
    ),
    oak = m.oak.eq(
      parameters["L_mix_1"] / sum_m,
      parameters["L_mix_2"] / sum_m,
      parameters["L_mix_3"] / sum_m,
      mmax,
      1.0 - S(forcing_data["swc"],
        parameters["L_swc_1"],
        parameters["L_swc_2"]),
      1.0 - S(forcing_data["tloss"],
        parameters["L_tloss_1"],
        parameters["L_tloss_2"])
    )
  )
}

induce <- function(version, state, parameters, constants, forcing_data){
  Ms <- state["Ms"]
  Mr <- state["Mr"]
  Cs <- state["Cs"]
  Cr <- state["Cr"]
  Ns <- state["Ns"]
  Nr <- state["Nr"]

  Km <- constants["KM"]
  Ka <- constants["KA"]
  Jc <- constants["Jc"]
  Jn <- constants["Jn"]
  fc <- constants["Fc"]
  fn <- constants["Fn"]

  #Divergence from technical specification
  a <- a.eq(version, state, parameters, constants, forcing_data)
  b <- b.eq(version, state, parameters, constants, forcing_data)

  Uc <- eq.Uc(a, Ms, Ka, Cs, Jc)
  Un <- eq.Un(b, Mr, Ka, Nr, Jn)

  g <- g.eq(version, state, parameters, constants, forcing_data)

  Gs <- eq.Gs(g, Cs, Ns, Ms)
  Gr <- eq.Gr(g, Cr, Nr, Mr)

  Tauc <- eq.Tauc(Ms, Mr, Cs, Cr)
  Taun <- eq.Taun(Ms, Mr, Ns, Nr)

  Cs <- eq.Cs(Cs = Cs, Uc = Uc, fc = fc, Gs = Gs, Tauc = Tauc)
  Cr <- eq.Cr(Cr, Tauc, fc, Gr)

  Ns <- eq.Ns(Ns, Taun, fn, Gs)
  Nr <- eq.Nr(Nr, Un, fn, Gr, Taun)

  m <- m.eq(version, state, parameters, constants, forcing_data)

  Ls <- eq.Ls(m, Ms, Km)
  Lr <- eq.Lr(m, Mr, Km)

  Ms <- eq.Ms(Ms, Gs, Ls)
  Mr <- eq.Mr(Mr, Gr, Lr)

  list(
    next.state = positive(c(
      Ms = Ms,
      Mr = Mr,
      Cs = Cs,
      Cr = Cr,
      Ns = Ns,
      Nr = Nr
    )),
    out = c(
      Ms = state["Ms"],
      Mr = state["Mr"],
      Cs = state["Cs"],
      Cr = state["Cr"],
      Ns = state["Ns"],
      Nr = state["Nr"],
      TAUc = Tauc,
      TAUn = Taun,
      Uc = Uc,
      Un = Un,
      Gs = Gs,
      Gr = Gr,
      Ms_loss = Ls,
      Mr_loss = Lr,
      CUR = a,
      NUR = b,
      g = g,
      loss = m
    )
  )
}

state_variables <- c("Ms", "Mr", "Cs", "Cr", "Ns", "Nr")

produce.expectation <- function(parameters, Data){
  output <- rep(0, prod(Data$out.dim))
  dim(output) <- Data$out.dim
  dimnames(output) <- Data$out.dimnames
  nspecies <- dim(output)["species"]
  for(site in 1:dim(output)["sites"]){
    initial_mass <- matrix(
      data = ifelse(
        Data$options$initial_mass_par,
        parameters$beta["iM",],
        rep(Data$options$initial_mass, nspecies)
      ),
      nrow = 6,
      ncol = nspecies
    )
    state <- initial_mass * matrix(
      data = c(1,1,.05,.05,.01,.01),
      nrow = length(state_variables),
      ncol = nspecies)
    dimnames(state) <- list(state_variables,Data$options$species)
    for(step in 1:dim(output)["steps"]){
      for(species in 1:nspecies){
        ts <- Data$timeseries[,step,site]
        ti <- Data$time.invariant[,site]
        afqr <- unname(
          switch(
            as.character(Data$options$photo[species]),
            c3 = ts["C3p"],
            c4 = ts["C4p"]
          )
        )
        forcing_data <- c(ts, ti, afqr = afqr)
        constants <- Data$globals
        out <- induce(
          as.character(Data$options$ve),
          state[,species],
          c(parameters$alpha, parameters$beta[,species]),
          constants,
          forcing_data)
        output[,species,step,site] <- out$out
        state[,species] <- out$next.state
      }
    }
  }
  output[,,Data$options$steps,]
}

test_that("Compiled version has no significant errors", {
  #create simulated data
  nspecies <- 3
  nsites <- 6
  ntime <- 99

  observations <- 1:(nspecies*nsites*ntime)
  obsdim <- c(ntime,nspecies,nsites)
  names(obsdim) <- c("time","species","site")
  dim(observations) <- obsdim

  wp <- rnorm(nsites, 1, 0.1)
  fc <- rnorm(nsites, 10,0.5)

  gauss_noise <- function(sd){
    rnorm(ntime, 0, sd)
  }

  maxtemp<- rnorm(nsites, 15, 5)
  mintemp<- rnorm(nsites, 8, 5)
  moist <- pmin(100, pmax(rnorm(nsites, 50,40), 0))

  per <- function(min, max, freq, sd){
    min + (max-min)/2 * (sin((1:ntime)*2*pi/(ntime/freq))+1) + gauss_noise(sd)
  }
  swc <- sapply(moist, function(p){pmin(100,pmax(per(p/10,p,2,10),0))})
  rsds <- sapply(1:nsites, function(a){pmax(per(0,100,3,10),0)})
  minmax <- rbind(mintemp,maxtemp)
  cols <- split(minmax, col(minmax))
  tcur <- sapply(
    cols,
    function(col){per(col[[1]],col[[2]],5,5)+5}
  )
  tnur <- sapply(
    cols,
    function(col){per(col[[1]],col[[2]],5,5)}
  )
  tgrowth <- sapply(
    cols,
    function(col){per(col[[1]],col[[2]],5,5)+2}
  )
  tloss <- sapply(
    cols,
    function(col){per(col[[1]],col[[2]],5,5)-3}
  )
  inp <- get_input(
    observations = observations,
    seconds = 1,
    p3 = TTR.PGM::p3,
    p4 = TTR.PGM::p4,
    lon = c(1:nsites),
    lat = c(1:nsites),
    swc=swc,
    catm=matrix(200,nrow=ntime, ncol=nsites),
    rsds=rsds,
    tcur=tcur,
    tnur=tnur,
    tgrowth = tgrowth,
    tloss = tloss
  )
  for(ver in c("std", "fqr", "red", "oak")){
    opt <- standard_options
    opt$steps <- 1:(ntime-1)
    opt$species <- c("1", "2","3")
    opt$photo <- c("c3","c4","c3")
    opt$ve <- ver

    optmin <- opt
    optmin$steps <- nsites

    Data <- make_data(inp, options = opt,globals = standard_globals)
    some.parms <- rnorm(n=Data$n.parm, mean=(Data$bounds[,"lower"] + Data$bounds[,"upper"]) / 2, sd=1.2)
    par <- get_parms(some.parms, Data)
    x <- run_ttr(par, Data)
    xexp <- produce.expectation(par, Data)
    expect(all(abs(x["Ms",,,]-xexp["Ms",,,])/mean(x["Ms",,,]) < 1e-10) , paste0(ver, ": Ms error is OK"))
    expect(all(abs(x["Mr",,,]-xexp["Mr",,,])/mean(x["Mr",,,]) < 1e-10) , paste0(ver, ": Mr error is OK"))
  }
})