### Function to generate, fit and plot Fisher's Logseries ###

## Discrete birth-death process with imigration
## Births, deaths and migrant arrival are Poisson processes
## Note that changes are not restricted to one at each time-step
## m= colonization ratio: parameter of a Poisson distributed number of arriving immigrants. Each immigrant belongs to a different species
## b = per capita birth rate
## d = per capita death rate
## tmax = number of time intervals to run
## N0 = starting abundances of each species
## S = starting number of species in the model
## When b<=d and m is small this simulation approaches Caswell's model 1 and thus produces a logseries distribution
disc.bde <- function(b=0.11,d=0.1,m=0.01,N0=1,S=20,tmax=1e4){
  j <- rep(N0,S)
    for(i in 1:tmax){
      n.migr <- rpois(1,m)
      j <- c(j,rep(1,n.migr))
      S <- length(j)
      bir <- rpois(S,b*j)
      dea <- rbinom(S,size=j,p=d)
      j <- j+bir-dea
      j <- j[j>0]
    }
  j
}

## Same function, but birth and death rates are uniform random variables
## death rate has expected value of 0.5
## b = expected birth rate
disc.bder <- function(b=0.11,m=0.01,N0=1,S=20,tmax=1e4){
  j <- rep(N0,S)
    for(i in 1:tmax){
      n.migr <- rpois(1,m)
      j <- c(j,rep(1,n.migr))
      S <- length(j)
      br <- runif(S,b-b/2,b+b/2)
      dr <- runif(S,0,1)
      bir <- rpois(S,br*j)
      dea <- rbinom(S,size=j,p=dr)
      j <- j+bir-dea
      j <- j[j>0]
    }
  j
}



## Predicted number of species for a given abundance (r)
pred.logser=function(r,alfa,N){
	X=N/(alfa+ N)
	alfa*X^r/r
}

## Logseries PDF
## Suppl Mat. of Alonso et al. 2008 Ecol.Lett. 11:93-105
dlr <- function(n,N,alfa){
		x<-N/(N+alfa)
		gama<- function(x) {1/log(1/(1-x))}
		gama(x)*(x^n)/n
              }

## Fit logseries given an abundance vector or
## just the number of individuals and species
## if abundance vector is given, returns also the log-likelihood and AIC 
fit.logser=function(x=NULL,size=NULL,rich=NULL){
  if(!is.null(x)){
    S <- length(x)
    N <- sum(x)
    if(!is.null(size)|!is.null(rich)){
      warning(paste("Model fitted with size = ",N," and rich = ",S," \n calculated from supplied abundances"))
    }
  }
  if(is.null(x)&!is.null(size)&!is.null(rich)){
    S <- rich
    N <- size
  }

  if(is.null(x)&is.null(size)&is.null(rich)){
    stop("Please provide size and species number or a vector of species abundances")
  }
  f1 <- function(a) {
        S + a * log((a/(a + N))) #thanks to function fishers.alpha of untb package
    }
  sol <- uniroot(f1, interval = c(1/N, N))
  alfa=sol$root
  X <- N/(N+alfa)
  results <- list(summary=c(N=as.integer(N),S=as.integer(S),alfa=alfa,X=X,loglik=NULL,AIC=NULL),sol=sol)
  if(!is.null(x)){
    LL <- sum(log(dlr(x,N,alfa)))
    aic <- -2*LL+2
    results$summary <- c(results$summary,loglik=LL,AIC=aic)
    results
  }
  results
}

## PDF for the continuous version of logseries (Tab.1 de Green & Plotkin 2006 Ecol Letters) ##
## Acessory functions: Incomplete Gamma upper limit
## Integral of x from infinity of the function  e^-t * t^(alfa-1)
## See also the function Igamma in the zipfR package,
## and the recipe in the help page of the function dgamma
## But none of these works for alfa=0, which is necessary to the equation of Green & Plotkin
igama.up <- function(alfa,x){
  f1 <- function(t1,a) {exp(-t1)*t1^(alfa-1)}
  integrate(f1,upper=Inf,lower=x)$value}
## Continuous PDF as function of parameter X
clr <- function(n,X){
  c1 <- igama.up(0,-log(X))
  (X^n) / (n*c1) 
}
##Continuous PDF as function of parameter alfa and total number of individual (N)
clr2 <- function(n,alfa,N){
  X<-N/(N+alfa)
  c1 <- igama.up(0,-log(X))
  (X^n) / (n*c1) 
}


## Species x abundance functions:
## Expected number of species for a given number of individuals, given alpha and N
S.logser <- function(N,alfa){
  alfa*log(1+N/alfa)
}
## Expected number of individuals, given alpha and number of species
N.logser <- function(S,alfa){
alfa*exp(S/alfa)-alfa
}

## Fit the expected curve of the logseries in a rank-abundance plot
## This function just return the predicted points form the cumulative number of species predicted by the logseries
rad.logser <- function(x,alfa=NULL){
  if(is.null(alfa)){
    a <- fit.logser(x)$summary["alfa"]
  }
else a <- alfa
  N <- sum(x)
  data.frame(x=cumsum(pred.logser(N:1,a,N)),y=N:1)
}

##This function creates a spline function that can be used with the function curve to add the curve to the rank-abundance plot
rad.logserf <- function(x,alfa=NULL){
  if(is.null(alfa)){
    a <- fit.logser(x)$summary["alfa"]
  }
else a <- alfa
  N <- sum(x)
  splinefun(x=cumsum(pred.logser(N:1,a,N)),y=N:1,"natural")
}