#' Run the simulation
#' @param max_geno integer : maximum genotype
#' @param range array : range
#' @param No integer : initial pop size
#' @param log10Pextlim integer (negative) : Log10 of the probability of extinction under which we consider that pop is rescued
#' @param rmax double : growth rate in exponential phase at the optimum (ro ~ Bmax - d)
#' @param factsecuext double :
#' @param nbt double : number of birth before testing for a rescue
#' @param ne integer :
#' @param n integer :
#' @param Es double :
#' @param Uo double :
#' @param na integer :
#' @param natba integer : Min nb of genotype copies to consider this genotype as advantageous
#' @param nbreps integer :
#' @param have_rescuers logical : if you want rescuers file
#' @param have_rate logical : if you want rate file (you need have_rescuers=TRUE)
#' @param output_path string : the output path
#' @param verbose logical : if you want a verbose program
#' @return array :
simulation <- function(max_geno = 500, range = c(-0.05, -0.1, -0.15),
  No = 100, log10Pextlim = -10, rmax = 1, factsecuext = 10, nbt = 100,
  ne = 10, n = 5, Es = 0.01, Uo = 0.01, na = 2000, natba = 10, nbreps = 50000,
  have_rate = FALSE, have_rescuers = FALSE, output_path = "", verbose = FALSE)
{
  library(MASS);
  library(DistributionUtils);

  print("Simulation Start")

  rescuersL <- list()

  d <- rmax                          # DOUBLE.
  Bmax <- rmax + d                   # DOUBLE. bmax in order to get ro=Bmax - d
  dt <- 50 / (No * (Bmax + d))       # DOUBLE.
  rovals <- range * rmax
  sovals <- rmax - rovals

  nsovals <- length(sovals)
  simname <- paste("No", No, "_n", n, "_ne", ne, "_U", Uo, "_Es", Es, "_so", sep="")
  #Rprint(verbose, simname)
  Ur <- numeric()
  q <- numeric()

  # loop on the values of so
  for(iso in 1:nsovals) {
    # set input
    so <- sovals[iso]
    tmax <- factsecuext * log(No+1) / abs(Bmax - so - d);
    nbtimesteps = floor(tmax / dt) + 2; # nb time step printed (+ 2 because we include time 0 and t_max)
    param <- c(n, ne, na, No, Bmax, d, nbreps, Uo, tmax, log10Pextlim, nbt, so, nbtimesteps, Es, dt, natba)

    # build the landscape
    Ini <- Rbuiltlandscape(n, ne, Es, so, na)

    #  Ur <- c(Ur, Uo / 2 * (1 - exp(- (B(Bmax, Ini$A, Ini$xo) - d) / d)))
    L <- Ini$Sb %*% Ini$M
    q <- c(q, (Ini$xo %*% L %*% L %*% Ini$xo) / (Ini$xo %*% L %*% Ini$xo) * Rtrace(L) / Rtrace(L %*% L))

    # birth and death process replicated nbreps times
    print(paste("running so=", so, "... Start", sep=""))

    # ============================
    # Start C++
    # ============================

    rescuersL[[iso]] <- Rrun(max_geno, n, ne, na, No, Bmax, d, nbreps, Uo, tmax, log10Pextlim, nbt, so, nbtimesteps, Es, dt, natba, Ini$Sb, Ini$A, Ini$xo, have_rate, have_rescuers, output_path, verbose)

    # ============================
    # End C++
    # ============================

    print(paste("running so=", so, "... End", sep=""))
  }

  pdf(paste(output_path,"distT1stepNo",No,"Uo",Uo,"n",n,"ne",ne,".pdf",sep=""), height=5, width=15)
  layout(matrix(1:3, 1, 3))

  for(iso in 1:nsovals) {

      Prescuers <- rescuersL[[iso]][["P"]]
      Trescuers <- rescuersL[[iso]][["T"]]
      NMrescuers <- rescuersL[[iso]][["NM"]]

      Trescuershist <- numeric()
      for (i in 1:length(Trescuers))

        if (Prescuers[i]==1) {
          Trescuershist <- c(Trescuershist, Trescuers[i])
        }
      hist(Trescuershist, freq=FALSE, breaks=50, col="orange", border=FALSE, main=paste("so =", sovals[iso]))
      xx <- density(Trescuershist)$x;

      points(dexp(0:500, abs(rovals[iso])), col="black", type="l", lwd=2)
  }

  dev.off()

  Prescuers <- rescuersL[[nsovals]][["P"]]
  NMrescuers <- rescuersL[[nsovals]][["NM"]]
  rescued <- NMrescuers[Prescuers == 1]
  res <- c()
  for (i in 0:max(rescued)) {
    res <- c(res, length(rescued[rescued==i]))
  }

  print(res)

  print("Simulation Stop")

  return(res)
}

Rprint <- function (verbose, message)
{
  if(verbose) {
    print(message)
  }
}

Rtrace <- function(matrice)
{
  return(sum(diag(matrice)))
}

# building M & S as random matrices
# the distance to optimum xo
# the set of nba mutation phenotypic effects drawn into MVN(0,M)
# it is an K-allele-model with very large K = 2000
Rbuiltlandscape <- function(n, ne, Es, so, na)
{
  if(ne >= n) m <- 1000 else m <- round(2 * n * ne / (n-ne) )

  Xs <- matrix(rnorm(n*m),ncol=m)
  S1 <- (1/m) * Xs %*% t(Xs)
  Xm <- matrix(rnorm(n*m), ncol=m)
  M1 <- (1/m) * Xm %*% t(Xm)

  #scaling such that E(s)=bmax*tr(SM)/2-d
  xx <- sqrt(2 * Es / Rtrace(S1 %*% M1) )
  M <- xx * M1
  S <- xx * S1

  # matrix of allelic effects
  A <- mvrnorm(na, numeric(n), M);

  xo <- numeric(n)
  if (so != 0) { #if so = 0, all the genotypic values stay at 0
    x1 <- rnorm(n)
    xo <- x1 * as.vector(sqrt(2 * so / (t(x1) %*% S %*% x1) ))
  }

  return(list(M=M, Sb=S, A=A, xo=xo))
}