maxEigen <- function(n,p,sd) {
  A <- matrix(0,nrow=n,ncol=n)
  for (i in seq(n))
    for (j in seq(n)) {
      if (i != j && runif(1) < p) A[i,j] <- rnorm(1,0,sd)
    }
  diag(A) <- -1
  #diag(A) <- -abs(rnorm(n,mean=1,sd=0.5))
  return(max(Re(eigen(A)$values)))
}

# n is the number of "species", p the connectivity and s the standard deviation:
n <- 100
p <- 1
s <- 0.25
maxEigen(n,p,s)   # maxEigen returns the largest eigenvalue
s*sqrt(n*p)       # this returns the value of Eq. (10.13)

# Run this for nPs different values of p until a maximum of maxP
# make nTrials different samples for every value of p:

maxP <- 0.25; nPs <- 10; nTrials <- 20
nstable <- rep(0,nPs)  # this will store the number of stable systems for every value of p
for (i in seq(nPs)) {
  p <- maxP*i/nPs
  print(c(P=p,Eq_May=s*sqrt(n*p)))
  for (j in seq(nTrials))
    nstable[i] <- nstable[i] + ifelse(maxEigen(n,p,s)<0,1,0)
}

plot(maxP*seq(nPs)/nPs,nstable/nTrials,xlab="P",ylab="Fraction stable",main=paste(n,s))
curve(s*sqrt(n*p),from=0,to=maxP,xname="p",add=TRUE)  # Plot Eq. (10.13) 


# Cui and Metha from MIT in Boston performed a similar analysis, but they kept the total
# interaction strength per species the same when they increased the connectivity of the
# matrix.  Highly connected systems would then have smaller off-diagonal elements than
# lowly connected systems.  This can be achieved by scaling the standard deviation of
# the normal distribution. Since the average of drawing n elements out of a normal
# distribution |N(0,1)|=0.7979, where 0 the mean, 1 the standard deviation, and |x| 
# means the absolute value, one can scale the standard deviation s=x/((n-1)*p*0.7979),
# where x is the desired total of a row of off-diagonal elements (called interTot in
# the loop below). Try different values of this interspecific competitions.

# interTot <- 4
# maxP <- 0.4; nPs <- 10; nTrials <- 20
# nstable <- rep(0,nPs)
# for (i in seq(nPs)) {
#   p <- maxP*i/nPs
#   s <- interTot/((n-1)*p*0.7979)
#   print(c(P=p,StDev=s,EqMay=s*sqrt(n*p)))
#   for (j in seq(nTrials))
#     nstable[i] <- nstable[i] + ifelse(maxEigen(n,p,s)<0,1,0)
# }
# plot(maxP*seq(nPs)/nPs,nstable/nTrials,xlab="P",ylab="Fraction stable",main=interTot)
# curve(interTot/((n-1)*p*0.7979)*sqrt(n*p),from=0,to=maxP,xname="p",add=TRUE)