KSEcology: Functions of Ecology from KeepSimpler

#' @title compute [x^*](variables in equilibrium) from parameters
#' @description the equilibrium function is:
#'     r - sx + \frac{k_m m x}{1 + h k_m m x} = 0
#' @param r - intrinsic growth rate
#' @param m - inter-species mutualistic strength
#' @param km - mutualistic degree of species,
#' @param h - handling time
#' @param s - intra-species competitive strength
#' @return x^* - variables in equilibrium or NULL - if
solve_nstars <- function(r, m, km, h, s) {
  A = s * h * km * m
  B = s - km * m - r * h * km * m
  C = - r
  delta = B^2 - 4 * A * C
  if (abs(delta) < 1e-12) { # precision error, should be 0
    X1 = - B / (2 * A)
    X2 = - B / (2 * A)
  } else if (delta < 0) { # negative delta
    return(NULL)
  } else { # positive delta
    X1 = (-B + sqrt(B^2 - 4 * A * C)) / (2 * A)
    X2 = (-B - sqrt(B^2 - 4 * A * C)) / (2 * A)
  }
  #delta2 <- (h * gamma * degree * r + s + gamma * degree)^2 - 4 * s * gamma * degree
  lev1 <- X1 * (km * m / (1 + h * km * m * X1)^2 - s)
  lev2 <- X2 * (km * m / (1 + h * km * m * X2)^2 - s)
  drdx1 <- s - km * m / (1 + h * km * m * X1)^2
  dr2dx21 <- 2 * h * (km * m)^2 / (1 + h * km * m * X1)^3
  return(c(r = r, m = m, km = km, h = h, s = s, X1 = X1, X2 = X2, delta = delta, drdx1 = drdx1, dr2dx21 = dr2dx21, lev1 = lev1, lev2 = lev2))  #delta2 = delta2,
}

#' @title the cusp curve of cusp catastrophe
#' -(\sqrt(k_mm) - \sqrt(s))^2/(hk_mm)
cusp_curve <- function(s, km, m, h) {
  return(-(sqrt(km * m) - sqrt(s))^2 / (h * km * m))
}


#' @title calculate Jocobian matrix in equilibrium of model [model_lv2_cm] in mean case
#' @param m - inter-species mutualistic strength
#' @param c - inter-species competitive strength
#' @param km - mutualistic degree of species,
#' @param h - handling time
#' @param s - intra-species competitive strength
#' @param graph - the adjacency matrix of network structure, generated by function [gen_regular_competition_cooperation]
#' @param nstar - species densities in equilibrium
#' @return Jacobian matrix
get_jacobian_mean <- function(m, c, km, h, s, graph, nstar) {
  J = graph
  # competitive part of Jacobian
  J[J < 0] <- - c * nstar
  # mutualistic part of Jacobian
  m_tilde <- m * nstar / (1 + h * km * m * nstar)^2 # element of mutualistic part
  J[J > 0] <- m_tilde
  # diagonal elements of Jacobian
  diag(J) <- - s * nstar
  return(J)
  }

#' @title calculate the largest eigenvalue of Jacobian
#' @param J the Jacobian matrix
#' @return estimated and real largest eigenvalue
calc_lev <- function(J) {
  lambda1 <- eigen(J)$values[1]
  # competitive part of Jacobian
  Jc <- J
  Jc[Jc > 0] = 0
  lambda2.competition = eigen(Jc)$values[1] # (- s + c) * nstar
  # mutualistic part of Jacobian
  Jm <- J
  Jm[Jm < 0] = 0
  m_tilde <- unique(Jm[Jm > 0]) # element of mutualistic part
  Jm.bin <- Jm
  Jm.bin[Jm.bin > 0] = 1
  lambda2.mutual <- estimate_second_eigenvalue(Jm.bin) * m_tilde
  lambda_ellipse <- lambda2.mutual + lambda2.competition
  lambda_dot <- unique(rowSums(J))
  c(lambda1 = lambda1, lambda_dot = lambda_dot, lambda_ellipse = lambda_ellipse, lambda2.competition = lambda2.competition)
}

#' @title generate simulation coefficient
gen_sim_coeffs <- function(n1, n2, r.max, s, c, kc, m, km, h) {
  s1 = n1  # plants number
  s2 = n2 # animal number
  k = km  # node degree
  alpha.row.mu = r.max # initial intrinsic growth rate
  alpha.row.sd = 0.
  alpha.col.mu = r.max
  alpha.col.sd = 0.
  beta0.mu = s # intra-species competition
  beta0.sd = 0
  beta1.mu = c  # inter-species competition
  beta1.sd = 0
  gamma.mu = m # mutualism
  gamma.sd = 0
  h.mu = h  # handling time
  h.sd = 0.
  delta = -0. # trade-off between mutualistic interaction strength and number
  list(s1 = s1, s2 =s2, k = k, alpha.row.mu = alpha.row.mu, alpha.row.sd = alpha.row.sd, alpha.col.mu = alpha.col.mu, alpha.col.sd = alpha.col.sd, beta0.mu = beta0.mu, beta0.sd = beta0.sd, beta1.mu = beta1.mu, beta1.sd = beta1.sd, gamma.mu = gamma.mu, gamma.sd = gamma.sd, h.mu = h.mu, h.sd = h.sd, delta = delta)
}

#' @title simulate the pressure experiment
sim_bipartite <- function(r, s, c, kc, m, km, h, n1, n2, graph) {
  r.max <- r[1]
  r.stepwise <- r[1] - r[2]
  r.steps <- length(r)
  coeffs <- gen_sim_coeffs(n1, n2, r.max, s, c, kc, m, km, h)

  #graph[graph < 0] = 0 # keep mutualistic part, remove competitive part
  params <- params_lv2_bipartite_new(graph = graph, coeffs = coeffs)

  # initial nstar
  #nstars <- solve_nstars(r.max, m, km, h, (s + kc * c))
  init <- runif2(n1 + n2, 1, 0)
  times <- 1:2000  # 1:100  1:500

  out <- sim_ode_press(model = model_lv2_cm, init = init, params = params, times = times, perturb = perturb_growthrate, perturbNum = r.steps + 10, isout = T, extinct_threshold = extinct_threshold, r.delta.mu = r.stepwise, r.delta.sd = 0) #r.stepwise/2

}

#' @title solve r2 in the case of B^2 - 4AC = 0
fold2_delta_root <- function(m12, m21, s1, s2, h, r1) {
  alpha1 <- s1 * m12 * m21 + s1 * s2 * m21
  alpha2 <- s2 * m12 * m21 + s1 * s2 * m12
  A0 <- r1 * h * s2 * m12 * m21 + s2 * m12 * m21 + s1 * s2 * m12
  A = (s1 * h * s1 * m12 * m12 + m12 * (1 + h * r1) * h * s1 * m12 * m21)^2
  B = 2 * (s1 * h * s1 * m12 * m12 + m12 * (1 + h * r1) * h * s1 * m12 * m21)       * (s1 * alpha2 + m12 * (1 + h * r1) * alpha1) -
      4 * s1 * m12 * A0 * h * s1 * m12 * m21
  C = (s1 * alpha2 + m12 * (1 + h * r1) * alpha1)^2 -
      4 * s1 * m12 * A0 * alpha1
  r2min1 = (-B + sqrt(B^2 - 4 * A * C)) / (2 * A)
  r2min2 = (-B - sqrt(B^2 - 4 * A * C)) / (2 * A)
  c(m12 = m12, m21 = m21, s1 = s1, s2 = s2, h = h, r1 = r1, r2min1 = r2min1, r2min2 = r2min2)
}

# fold2_delta_root <- function(m12, m21, s1, s2, h, r1) {
#   alpha1 <- s1 * m12 * m21 + s1 * s2 * m21
#   alpha2 <- s2 * m12 * m21 + s1 * s2 * m12
#   A = m12^2 * (1 + h * r1)^2 + s1^2 + 2 * s1 * m12 * (1 + h * r1)
#   B = m12^2 * (1 + h * r1)^2 * 2*alpha1 + s1^2 * 2*alpha2 +
#     2 * s1 * m12 * (1 + h * r1) * (alpha1 + alpha2) -
#     4 * s1 * m12 * (r1 * h * s2 * m12 * m21 + alpha2)
#   C = m12^2 * (1 + h * r1)^2 * alpha1^2 + s1^2 * alpha2^2 +
#     2 * s1 * m12 * (1 + h * r1) * alpha1 * alpha2 -
#     4 * s1 * m12 * (r1 * h * s2 * m12 * m21 + alpha2) * alpha1
#   r2min1 = (-B + sqrt(B^2 - 4 * A * C)) / (2 * A)
#   r2min1 = r2min1 / (h * s1 * m12 * m21)
#   r2min2 = (-B - sqrt(B^2 - 4 * A * C)) / (2 * A)
#   r2min2 = r2min2 / (h * s1 * m12 * m21)
#   c(m12 = m12, m21 = m21, s1 = s1, s2 = s2, h = h, r1 = r1, r2min1 = r2min1, r2min2 = r2min2)
# }

#' @title plot a fold bifurcation, for the nonstructural model II,
#'  control param is r2, state variable is x
#' @param m12, m21 total mutualistic strengths between group 1 and 2
#' @param s1, s2 total competitive strength in group 1 and 2
#' @param h, handling time
#' @param r1, intrinsic growth rate of group 1
#' @example plot_fold2_rx(m12 = , m21 = , s1 = , s2 = , h = 0.3, r1 = )
plot_fold2_rx <- function(m12, m21, s1, s2, h, r1) {
  r2min <- fold2_delta_root(m12, m21, s1, s2, h, r1)['r2min1'] + 1e-10
  #xmin = sqrt(s) * (sqrt(m) - sqrt(s)) / (s * h * m)
  #xr0 <- -(s - m) / (s * h * m)
  r2s = seq(from = r2min, to = 0, length.out = 200)
  require(plyr)
  tmp <- ldply(r2s, function(r2) {
    A = r1 * h * s2 * m12 * m21 + s2 * m12 * m21 + s1 * s2 * m12
    B = r1 * s2 * m21 - r2 * s1 * m12
    C = r2 * h * s1 * m12 * m21 + s1 * m12 * m21 + s1 * s2 * m21
    A2 = s1 * h * m12 * A / C # Avoid C == 0 !
    B2 = s1 * h * m12 * B / C - h * r1 * m12 + s1 * A / C - m12
    C2 = s1 * B / C - r1
    delta = B2^2 - 4 * A2 * C2
    X21 = (-B2 + sqrt(B2^2 - 4 * A2 * C2)) / (2 * A2)
    X22 = (-B2 - sqrt(B2^2 - 4 * A2 * C2)) / (2 * A2)
    X11 = (A * X21 + B) / C
    X12 = (A * X22 + B) / C
    c(r2 = r2, X21 = X21, X22 = X22, X11 = X11, X12 = X12)
  })
  matplot(tmp$r2, tmp[,c('X11', 'X12')], type = 'l', lwd = 2, xlab = 'r2', ylab = '(X1,X2)')
  matlines(tmp$r2, tmp[,c('X21', 'X22')], type = 'l', lwd = 2)
  lines(c(r2min*1.05, 0), c(0, 0), lwd = 2)
#   r1 <- tmp[nrow(tmp)/3,]$r
#   x1r1 <- tmp[nrow(tmp)/3,]$X1
#   x2r1 <- tmp[nrow(tmp)/3,]$X2
#   r2 <- tmp[nrow(tmp)*2/3,]$r
#   x1r2 <- tmp[nrow(tmp)*2/3,]$X1
#   x2r2 <- tmp[nrow(tmp)*2/3,]$X2
#   arrows(-0.0, xr0*0.2, -0.0, xr0*0.8, length = 0.1, lty = 1)
#   arrows(r1, x2r1+(x1r1-x2r1)*0.1, r1, x2r1+(x1r1-x2r1)*0.9, length = 0.1, lty = 1)
#   arrows(r2, x2r2+(x1r2-x2r2)*0.1, r2, x2r2+(x1r2-x2r2)*0.9, length = 0.1, lty = 1)
#   arrows(rmin, xmin*0.8, rmin, xmin*0.2, length = 0.1, lty = 1)
#   arrows(r1, x2r1*0.9, r1, x2r1*0.1, length = 0.1, lty = 1)
#   arrows(r2, x2r2*0.9, r2, x2r2*0.1, length = 0.1, lty = 1)
#   points(0, 0, pch = 16)
#   text(0, 0, '(0,0)', pos = 3)
#   points(0, xr0, pch = 16)
#   text(0, xr0, expression(paste("(0,",x[0],")")), pos = 1)
#   points(rmin, xmin, pch = 16)
#   text(rmin, xmin, expression(paste("(", r[min], ',', x[min],")")), pos = 4)
}