flan: Fluctuation Analysis of Mutant Cell Counts

    ## --------------------------------------------------------------------------------------------------------------------------
    ############################################## Estimations and tests functions ##############################################
    ## --------------------------------------------------------------------------------------------------------------------------

  # Estimation function of the mean number of mutations. Compute also the probability of mutation if with.prob is TRUE and the fitness parameter if fitness is empty.
    # Arguments:
    # mc:  a (non-empty) numeric vector of mutants counts.
    # fn:  an optional (non-empty) numeric vector with same length as mc of final numbers of cells.
    # mfn:  mean final number of cells which. Ignored if fn is non-missing.
    # cvfn:  coefficient of variation of final number of cells. Ignored if fn is non-missing.
    # fitness:  fitness parameter: ratio of growth rates of normal and mutant cells. If fitness is NULL (default) then the fitness will be estimated. Otherwise, the given value will be used to estimate the mean mutation number mutations
    # death:  death probability. Must be smaller than 0.5.
    # method:  estimation method as a character string: one of ML (default), P0, or GF.
    # winsor:  winsorization parameter: positive integer. Only used when method is ML or when method is P0 and fitness is NULL.
    # model:  statistical lifetime model as a character string: one of LD (default) for Luria-Delbruck model with exponential lifetimes, or H for Haldane model with constant lifetimes.
    # Assumptions for the estimation
    #   - fitness (if given)
    #   - death probability
    #   - lifetimes distribution model : exponentialy distributed lifetime ("LD") or constant lifetime "H"
    #   - sample of mutants counts
    #   - sample of finals counts (or a mean and coefficient of variation for finals counts)
    # Three possible methods :
    #   - p0-method ("P0")   (winsorize the mc sample if rho is estimated with the parameter winsor)
    #   - Generating function method ("GF")
    #   - Maximum of Likelihood method ("ML")  (winsorize the mc sample with the parameter winsor)
    #
    # Returns :
    #   - Estimate of alpha (or pi) and rho if non given
    #   - Standard deviation of alpha, (or pi) and rho if non given

mutestim <- function(mc,fn=NULL,mfn=NULL,cvfn=NULL,                 # user's data
                  fitness=NULL,death=0.,                # user's parameters
                  method=c("ML","GF","P0"),winsor=512, # estimation method
                  model=c("LD","H"))                   # clone growth model
                  {

    if(missing(method)){method <- "ML"}
    if(missing(model)){model <- "LD"}
    method <- match.arg(method)
    model <- match.arg(model)

    if(is.null(mc)){
      stop("'mc' is empty...")
    }
    if(sum(floor(mc)==mc) != length(mc)) stop("'mc' must be a vector of integer.")
    if(!is.null(fn)){
      if(length(fn) != length(mc)){
	stop("'fn must have the same length as 'mc'.")
      }
      mfn <- mean(fn)
      cvfn <- sd(fn)/mfn
      if(sd(fn) == 0) fn <- NULL
    }
    if(!is.null(mfn)){
      if(mfn < 0 | length(mfn) > 1){
	stop("'mfn' must be a single positive number.")
      }
      if(is.null(cvfn)) cvfn <- 0
    }
    if(!is.null(cvfn)){
      if(cvfn < 0 | length(cvfn) > 1){
	stop("'cvfn' must be a single positive number.")
      }
      if(is.null(mfn)) mfn <- 1e9
    }
    if(!is.null(fitness)){
      if(fitness < 0 | length(fitness) > 1){
	stop("'fitness' must be empty or a single positive number.")
      }
    }
    if(death < 0 | death >= 0.5 | length(death) > 1){
      stop("'death must be a single positive and < 0.5 number.")
    }
    if(winsor < 0 | trunc(winsor) != winsor | length(winsor) > 1){
      stop("'winsor' must be a single positive integer.")
    }
    if(method=="P0" & sum(mc==0) == 0 & death == 0){
      stop(paste("P0 method can not be used if 'mc' does not contain any null counts and if 'death' is null.",sep=" "))
    }
    if(is.null(fitness)){   # If fitness is empty : compiute estimates of mean number of mutations (mutation probaility), fitness

      # Maximum of Likelihood estimators
      if(method == "ML"){
	
	if(!is.null(fn)) output <- MutationProbabilityFitnessMLOptimization(mc=mc,fn=fn,model=model,death=death,winsor=winsor)
	else output <- MutationFitnessMLOptimization(mc=mc,mfn=mfn,cvfn=cvfn,model=model,death=death,winsor=winsor)
	
# 	if(!output$succeeds) warning("Initialization of 'fitness' with 'GF'-method may have failed.")
# 	output$succeeds <- NULL
      }
      # P0 method
      if(method == "P0") output <- MutationFitnessP0Estimation(mc=mc,fn=fn,mfn=mfn,cvfn=cvfn,model=model,death=death,winsor=winsor)
	
      # GF method
      if(method == "GF") {
	output <- MutationFitnessGFEstimation(mc=mc,mfn=mfn,cvfn=cvfn,death=death,model=model)
	if(!output$succeeds) warning(paste("Impossible to estimate 'fitness' with 'GF'-method : 'fitness' is set to default value 1 and only",if(!is.null(mfn)){"mutation probability"}else{"mutation number"},"is estimated.",sep=" "))
	output$succeeds <- NULL
      }
    } else {    # Else : compute estimate(s) of mean number of mutations, or mutation probability

      # Maximum of Likelihood estimator of mutations or mutprob
      if(method == "ML"){
	if(!is.null(fn)) output <- MutationProbabilityMLOptimization(mc=mc,fn=fn,model=model,fitness=fitness,death=death,winsor=winsor)
	else output <- MutationMLOptimization(mc=mc,mfn=mfn,cvfn=cvfn,model=model,fitness=fitness,death=death,winsor=winsor)


      }
      # P0 method
      if(method == "P0") output <- MutationP0Estimation(mc,mfn=mfn,cvfn=cvfn,death=death)     # P0 estimator of mutations
      

      # GF method
      if(method == "GF")  output <- MutationGFEstimation(mc=mc,mfn=mfn,cvfn=cvfn,fitness=fitness,death=death,model=model)

    }

    output
}


    # One-sample or two-sample tests for mean number of mutations, fitness and mutation probability using mutestim estimates
    # The possible assumptions are the same as for the mutestim function.
    # Arguments :
    # Returns a "flantest" object, which contains:
    # Tstat:  the value of the computed statistic.
    # parameter:  the values of the parameter of the model: fitness(if not tested), death, mfn (if needed) and cvfn
    # p.value:  the p-value(s) of the test.
    # conf.int:  a confidence interval for the parameter(s) of interest appropriate to the specified alternative hypothesis.
    # estimates:  the estimate(s) of interest.
    # null.value:  the specified hypothesized value(s).
    # alternative:  a character string describing the alternative hypothesis.
    # model:  the statistical lifetime model.
    # method:  method used to compute the estimation(s) of the parameter(s) of interest.
    # data.name:  a character string giving the name of the complete data.


flan.test <- function(mc,fn=NULL,mfn=NULL,cvfn=NULL,                      # user's data
               fitness=NULL,death=0.,                         # user's parameters
               mutations0=1,mutprob0=NULL,fitness0=1,       # null hypotheses
               conf.level=0.95,                              # confidence level
               alternative=c("two.sided","less","greater"),  # alternative
               method=c("ML","GF","P0"),winsor=512,          # estimation method
               model=c("LD","H"))                            # clone growth model
               {

  with.prob <- FALSE               # Boolean: if TRUE (if fn, mfn, or cvfn are given), mutprob is tested instead of mutations

  if(is.null(mc)){
    stop("'mc' is empty...")
  } else {
    if(is.list(mc)){
      if(length(mc) == 1){
	nsamples <- 1
	dname <- list(deparse(substitute(mc)))
	if(is.null(mc[[1]])) stop("'mc[[1]]' is empty...")
	if(is.list(fn)){
	  if(length(fn) != nsamples) stop("'mc' and 'fn' must have the same length.")
	  if(!is.null(fn[[1]])){
	    dname <- c(dname,deparse(substitute(fn)))
	    with.prob <- TRUE
	    mfn <- mean(fn)
	    cvfn <- sd(fn)/mfn
	  }
	} else {
	  if(!is.null(fn)) stop("'fn' must have the same type as 'mc'.")
	  fn <- list(fn)
	}
      }
      if(length(mc) == 2){
	nsamples <- 2
	dname <- list(c(paste(deparse(substitute(mc)),"1",sep=""),paste(deparse(substitute(mc)),"2",sep="")))
	if(is.null(mc[[1]])) stop("'mc[[1]]' is empty...")
	if(is.null(mc[[2]])) stop("'mc[[2]]' is empty...")

	if(is.list(fn)){
	  if(length(fn) != nsamples) stop("'fn' must have the same length as 'mc'.")

	  if(!is.null(fn[[1]])) {
	    if(length(fn[[1]]) != length(mc[[1]])) stop("'fn[[1]]' must have the same length as 'mc[[1]]'.")

	    if(is.null(fn[[2]])){
	      warning("'fn[[2]]' is empty : 'fn[[1]]' is ignored.")
	      fn <- list(NULL,NULL)
	    } else {
	      if(length(fn[[2]]) != length(mc[[2]])) stop("'fn[[2]]' must have the same length as 'mc[[2]]'.")

	      dname <- c(dname,list(c(paste(deparse(substitute(fn)),"1",sep=""),paste(deparse(substitute(fn)),"2",sep=""))))
	      with.prob <- TRUE
	      mfn <- unlist(lapply(fn,mean))
	      cvfn <- unlist(lapply(fn,sd))/mfn
	    }
	  } else {
	    if(!is.null(fn[[2]])) {
	      warning("'fn[[1]]' is empty : 'fn[[2]]' is ignored.")
	      fn <- list(NULL,NULL)
	    }
	  }
	} else {
	  if(!is.null(fn)) stop("'fn' must have the same type as 'mc'.")
	  fn <- list(fn,fn)
	}
      }
    } else {
      nsamples <- 1
      dname <- list(deparse(substitute(mc)))
      if(!is.null(fn)){
	if(is.list(fn)) stop("'fn' must have the same type as 'mc'.")
	if(length(fn) != length(mc)) stop("'fn' must have the same length as 'mc'.")
	dname <- c(dname,deparse(substitute(fn)))
	with.prob <- TRUE
	mfn <- mean(fn)
	cvfn <- sd(fn)/mfn
      }
      mc <- list(mc)
      fn <- list(fn)
    }
  }

  if(!is.null(mfn)){
    if(sum(mfn < 0) != 0 | length(mfn) > 2) stop("if given, 'mfn' must be a vector with size <= 2 of positive numbers.")

    if(is.null(cvfn)) cvfn <- rep(0,length(mfn))
    with.prob <- TRUE
  }
  if(!is.null(cvfn)){
    if(sum(cvfn < 0) != 0 | length(cvfn) > 2){
      stop("if given, 'cvfn' must be a vector with size <= 2 of positive numbers.")
    }
    if(is.null(mfn)) mfn <- rep(1e9,length(cvfn))
    with.prob <- TRUE
  }
  if(!is.null(fitness)){
    if(sum(fitness < 0) != 0 | length(fitness) > 2){
      stop("if given, 'fitness' must be a vector with size <= 2 of positive numbers.")
    }
    fitness0 <- NULL
  }
  if(length(mutations0) > 1 | mutations0 < 0){
    stop("'mutations0' must be a single positive number.")
  }
  if(!is.null(fitness0)){
    if(length(fitness0) > 1 | fitness0 < 0){
      stop("'fitness0' must be a single positive number.")
    }
  }
  if(!is.null(mutprob0)){
    if(length(mutprob0) > 1 | mutprob0 < 0 | mutprob0 >= 1){
      stop("if given, 'mutprob0' must be a positive and <= 1 number.")
    }
    with.prob <- TRUE
    if(is.null(mfn)) mfn <- 1e9
    if(is.null(cvfn)) cvfn <- 0
  } else {
    if(with.prob){
      if(is.null(mfn)) mfn <- 1e9
      if(is.null(cvfn)) cvfn <- 0
      if(nsamples == 1) mutprob0 <- mutations0/mfn
    }
  }


  # Default values if two-samples test
  if(nsamples == 2){
    if(missing(mutations0)) mutations0 <- 0
    if(is.null(fitness)){
      if(missing(fitness0)) fitness0 <- 0
    }
    if(with.prob) {
      if(missing(mutprob0)) mutprob0 <- 0
    }
  }

  if((length(conf.level) > 1) | conf.level > 1 | conf.level < 0){
    stop("'conf.level' must be a single positive and <= 1 numbers.")
  }
  if(sum(death < 0 | death >= 0.5) != 0 | length(death) > 2){
    stop("'death' must be a vector with size <= 2 of positive and < 0.5 numbers.")
  }
  if(winsor < 0 | trunc(winsor) != winsor | length(winsor) > 1){
    stop("'winsor' must be a single positive integer.")
  }

  H0 <- c(if(with.prob) mutprob0 else mutations0,fitness0)     # Vector of null hypothesises

  np <- length(H0)                                             # Number of tested values

  if(missing(alternative)) {alternative <- rep("two.sided",np)}
  if(missing(method)) {method <- "ML"}
  if(missing(model)) {model <- "LD"}

  alternative <- match.arg(alternative,several.ok=TRUE)
  method <- match.arg(method)
  model <- match.arg(model)


  if(nsamples == 1){

    parameter <- NULL
    names <- character()

    if(np == 1){
      names(H0) <- if(with.prob) "mutation probability" else "mutation number"
      parameter <- c(fitness,death)
      names <- c("fitness","death")
      if(with.prob){
	parameter <- c(parameter,mfn,cvfn)
	names <- c(names,"mfn","cvfn")
      }
    }
    if(np == 2) {
      names(H0) <- c(if(with.prob) "mutation probability" else "mutation number", "fitness")
      parameter <- death
      names <- "death"
      if(with.prob){
	parameter <- c(parameter,mfn,cvfn)
	names <- c(names,"mfn","cvfn")
      }
    }
    names(parameter) <- names
  }
  if(nsamples == 2){
    if(np == 1 & length(fitness) == 1){fitness <- c(fitness,fitness)}
    if(length(death) == 1){death <- c(death,death)}
    if(length(mfn) == 1){mfn <- c(mfn,mfn)}
    if(length(cvfn) == 1){cvfn <- c(cvfn,cvfn)}

    parameter <- NULL
    names <- character()


    if(np == 1){
      names(H0) <- if(with.prob) "mutprob difference" else "mutations difference"
      parameter <- cbind(fitness,death)
      names <- c("fitness","death")
      if(with.prob) {
	parameter <- cbind(parameter,mfn,cvfn)
	names <- c(names,"mfn","cvfn")
      }
    }
    if(np == 2){
      names(H0) <- c(if(with.prob) "mutprob difference" else "mutations difference","fitness difference")
      parameter <- cbind(death)
      names <- "death"
      if(with.prob){
	parameter <- cbind(parameter,mfn,cvfn)
	names <- c(names,"mfn","cvfn")
      }
    }
    colnames(parameter) <- names
  }

  if(is.null(mfn)) mfn <- list(mfn)
  if(is.null(cvfn)) cvfn <- list(cvfn)
  if(is.null(fitness)) fitness <- list(fitness)

   # Estimates mean number of mutations alpha, given fitness parameter rho
  ests <- mapply(function(x,y,m,c,f,d){
    mutestim(mc=x,fn=y,mfn=m,cvfn=c,method=method,model=model,fitness=f,death=d,winsor=winsor)
  },mc,fn,mfn,cvfn,fitness,death)


  if(np == 1){   # If only mutations (or mutprob) is tested
    if(length(alternative) > 1) alternative <- alternative[1]

    Tstat <- unlist(ests[1,])    # Extract the estimate to compute statistic of the test
    sds <- unlist(ests[2,])      # Standard deviation of the estimate
    names(Tstat) <- NULL
    ests <- Tstat                # Keep the estimate

    if(nsamples == 2){    # If Two-sample test
      ests <- cbind(ests)
      Tstat <- -diff(Tstat)      # Statistic of the test is build with difference of estimates
      sds <- sqrt(sum(sds^2))    # Standard deviation of the difference between estimates
      colnames(ests) <- if(with.prob) "mutprob" else "mutations"
    }
    names(ests) <- if(with.prob) "mutprob" else "mutations"


  }
  if(np == 2){    # If fitness is also tested
    if(length(alternative) == 1) alternative <- rep(alternative,2)
    if(length(alternative) > 2) alternative <- alternative[c(1,2)]

    Tstat <- rbind(unlist(ests[1,]),unlist(ests[3,]))   # Extract estimates to compute statistic of the test
    sds <- rbind(unlist(ests[2,]),unlist(ests[4,]))     # Standard deviation of the estimates
    colnames(Tstat) <- NULL
    ests <- cbind(Tstat[1,],Tstat[2,])                  # Keep the estimates

    if(nsamples == 2){               # If Two-sample test
      Tstat <-- apply(Tstat,1,diff)                     # Statistic of the test is build with difference of estimates
      sds <- apply(sds,1,function(s){sqrt(sum(s^2))})   # Standard deviation of the difference between estimates
    }
    colnames(ests) <- c(if(with.prob) "mutprob" else "mutations","fitness")
  }

  cint <- mapply(function(e,s,alt){
    if(alt == "less"){                                        # Confidence interval(s) of the test
      c(-Inf,e+s*qnorm(conf.level))                           # with respect to the confidence level
    } else if(alt == "greater"){                              # and the alternative
      c(e-s*qnorm(conf.level),Inf)
    } else if(alt == "two.sided"){
      e+s*c(-1,1)*qnorm((1+conf.level)/2)
    }
  },Tstat,sds,alternative)



  Tstat <- (Tstat-H0)/sds                                      # Statistic(s) of the test

  pval <- mapply(function(alt,tstat){
      if(alt == "less"){                                        # p-value(s) of the test
	pnorm(tstat)                                            # with respect to the alternative
      } else if(alt == "greater"){
	pnorm(tstat,lower.tail=FALSE)
      } else if(alt == "two.sided"){
	2*pnorm(-abs(tstat))
      }
    },alternative,Tstat)

    if(nsamples == 2){
      names(Tstat) <- sapply(colnames(ests),function(noun){paste(noun,"difference")})
    } else {
      names(Tstat) <- names(ests)
    }

    colnames(cint) <- names(Tstat)
    names(pval) <- names(Tstat)
    names(alternative) <- names(Tstat)

  rownames(cint) <- c("bInf","bSup")


  attr(cint,"conf.level") <- conf.level

  # Build object with 'flantest' class

  flantest(Tstat=Tstat,estimate=ests,parameter=parameter,conf.int=cint,p.value=pval,
  null.value=H0,alternative=alternative,data.name=dname,model=model,method=method,nsamples=nsamples)

}

    ## --------------------------------------------------------------------------------------------------------------
    ############################ Density, distribution function, quantile function and random #######################
    ################################## generation for the mutants counts ############################################
    ## --------------------------------------------------------------------------------------------------------------

    # There are two cases :
    #  - If alpha is given, then returns a sample mcs of mutants counts with parameters :
    #   * mutations : mean number of mutation of mutation
    #   * fitness : fitness parameter
    #   * death : probability of death
    #   * dist : lifetimes distribution
    #
    #  - If alpha is not given, then compute first a sample fn of final counts with a Log-Normal distribution with mean mfn and
    #    coefficient of variation cvfn. Then compute a sample mcs of mutants counts, where the ith count has the following parameters :
    #   * mutprob*fni : mean number of mutation
    #                      (where fni is the corresponding of final count, and pi the probability of mutation)
    #   * fitness : fitness parameter
    #   * death : probability of death
    #   *dist : lifetimes distribution
    #
    #   Returns the two samples mcs and fn
    #
rflan <- function(n,mutations=1,mutprob=NULL,fitness=1,death=0.,
        dist=list(name="lnorm",meanlog=-0.3795851,sdlog=0.3016223),
        mfn=1e9,cvfn=0) {

    if(n <= 0) stop("'n' must be a positive integer.")

    if(length(n) > 1) n <- length(n)

    if(mutations < 0 | length(mutations) > 1) stop("'mutations' must be a single positive number")

    if(fitness < 0 | length(fitness) > 1) stop("'fitness' must be a single positive number")

    if(mfn < 0 | length(mfn) > 1) stop("'mfn' must be a single positive number")

    if(cvfn < 0 | length(cvfn) > 1) stop("'cvfn' must be a single positive number")

    if(!is.list(dist)) stop("'dist' must be a list of a character chain followed by its arguments")


    names(dist)[1] <- "name"
    if(dist$name == "exp"){
      names(dist)[2] <- "rate"
    } else if(dist$name == "dirac"){
      names(dist)[2] <- "location"
    } else if(dist$name == "lnorm"){
      names(dist)[2] <- "meanlog"
      names(dist)[3] <- "sdlog"
    } else if(dist$name == "gamma"){
      names(dist)[2] <- "shape"
      names(dist)[3] <- "scale"
    } else stop("'dist[[1]]' must be a character chain among 'exp', 'dirac', 'lnorm', 'gamma'")


    if(death < 0 | death >= 0.5 | length(death) > 1) stop("'death' must be a single positive and < 0.5 number ")


    if(!is.null(mutprob)){
      if(mutprob < 0 | mutprob > 1 | length(mutprob) > 1) stop("'mutprob' must be a single positive and <= 1 number")
      mutations <- mutprob*(if(cvfn == 0) mfn else 1)              # Sample with mutprob instead of mutations if cvfn > 0
    } else {
      if(cvfn > 0) mutations <- mutations/mfn    # Default value of mutprob if missing and cvfn > 0
    }



      # output <- .Call("rflan",list(
			#   n=as.integer(n),
			#   mutations=as.double(mutations),
			#   fitness=as.double(fitness),
			#   death=as.double(death),
			#   distribution=dist,
			#   mfn=as.double(mfn),
			#   cvfn=as.double(cvfn)
			#   )
		  # )

    if(dist$name == "lnorm" | dist$name == "gamma")  dist <- adjust.rate(dist,1,death)
    flan.sim <- new(FlanSim,list(
      mutations=mutations,
      fitness=fitness,
      death=death,
      dist=dist,
      mfn=mfn,
      cvfn=cvfn
      ))
      
    output <-flan.sim$rflan(n)


    mc <- output$mc
    indNA <- which(mc < 0)
    if(length(indNA) > 0){
      warning("Production of NaNs in geometric sampling.")
      output$mc[indNA] <- NaN
    }
    
    if(cvfn == 0) output$fn <- rep(mfn,n)

    output
}

  # Quantile function of the mutants count with parameter
  #   mutations : mean number of mutations
  #   fitness : fitness parameter
  #   death : death probability
  #   model : lifetimes distribution model : exponentialy distributed lifetimes ("LD") or constant lifetimes ("H")

qflan <- function(p,mutations=1,fitness=1,death=0.,model=c("LD","H"),lower.tail=TRUE){

  if((sum(p < 0)+sum(p > 1)) != 0){
    stop("'p' must be a vector of positive and <= 1 numbers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  if(missing(model)){model="LD"}
  model <- match.arg(model)
#   if(!lower.tail) p <- 1-p


  m <- 100
  P <- pflan(0:m,mutations,fitness,death,model,lower.tail)
  sapply(p,function(pp){
    if(pp == 1) Inf
    if (lower.tail & pp <= P[1]) 0
    if (!lower.tail & pp >= P[1]) 0
    else {
      k <- if(lower.tail) max(which(P<pp)) else max(which(P>pp))          # quantile if in the actual table
      while (k>=m){                   # if p not yet in the table
	m2 <- 2*m                       # double the table
	P2 <- pflan(m:m2,mutations,fitness,death,model,lower.tail)      # truncated distribution function
	if (lower.tail & pp <= P2[1]) k
	if (!lower.tail & pp >= P2[1]) k
	else {
	  k <- k-1+if(lower.tail) max(which(P2<pp)) else max(which(P2>pp))           # quantile if in the table
	  m <- m2
	}
      }                              # end while
    }                                  # end if
    k
  })

}

  # Distribution function of the mutants count with parameter
  #   mutations : mean number of mutations
  #   fitness : fitness parameter
  #   death : death probability
  #   model : lifetimes distribution model : exponentialy distributed lifetimes ("LD") or constant lifetimes ("H")


pflan <- function(m,mutations=1,fitness=1,death=0.,model=c("LD","H"),lower.tail=TRUE){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  if(missing(model)){model="LD"}
  model <- match.arg(model)

#   output=.Call("pflan",
#     list(m=as.integer(max(m)),mutations=as.double(mutations),fitness=as.double(fitness),death=as.double(death),model=as.character(model),lowerTail=as.logical(lower.tail),bool=as.logical(length(m)>1)))

  flan.mutmodel <- new(FlanMutMod,list(
    mutations=mutations,
    fitness=fitness,
    death=death,
    model=model,
    lt=lower.tail
  ))
  
#   print(flan.mutmodel$getfcts())
  
  M=max(m)
  output <- flan.mutmodel$pflan(M)

  output[m+1]
}


  # Density function of the mutants count with parameter
  #   mutations : mean number of mutations
  #   fitness : fitness parameter
  #   death : death probability
  #   model : lifetimes distribution model : exponentialy distributed lifetimes ("LD") or constant lifetimes ("H")


dflan <- function(m,mutations=1,fitness=1,death=0.,model=c("LD","H")){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  if(missing(model)){model="LD"}
  model <- match.arg(model)

  
  flan.mutmodel <- new(FlanMutMod,list(
    mutations=mutations,
    fitness=fitness,
    death=death,
    model=model,
    lt=TRUE
  ))
  
  M <- max(m)
  output <- flan.mutmodel$dflan(M)

  output[m+1]
}



    ## --------------------------------------------------------------------------------------------------------------
    ############################################## Hidden functions #################################################
    ## --------------------------------------------------------------------------------------------------------------


## Sampling
adjust.rate <- function(dist,fitness=1,death=0.){
#   rescales the parameter(s) of distribution dist to obtain growh rate equal to fitness.
#   The distribution is given and returned as a list of a character chain 
#   followed by the list of parameters. 

  adist <- dist				# adjusted distribution
  m <- 2*(1-death)			# mean offspring number
  switch(dist[[1]],
    dirac={                              # Dirac distribution
	adist$location <- log(m)/fitness      # new location parameter
	return(adist)
	  },                              # end Dirac distribution
    exp={                                # Exponential distribution
	adist$rate <- fitness/(m-1)
	return(adist)
	  },   
    gamma={				# Gamma distribution
      a <- dist$shape			# according to the shape parameter
      adist$scale <- (m^(1/a)-1)/fitness	# adjust scale parameter 
    },
    lnorm={					# Log-normal distribution
#       s <- growth.rate(dist,death)		# according to the rate
      a <- dist$sdlog
      l <- dist$meanlog
      lap <- function(s){               # Laplace transform
	f <- function(x){exp(-(s*exp(a*x*sqrt(2)+l)+x^2))}

	I <- {integrate(f,
	      lower=-Inf,upper=Inf,rel.tol=.flantol,subdivisions=.flansubd)}      # integrate from - to + infinity
	
	I <- I$value                      # get the value
  	I <- I/sqrt(pi)			# rescale
        I-1/m                     # return shifted Laplace transform
      }
      lwb <- 1                          # lower bound of search interval
      while(lap(lwb)<0){lwb<-lwb/2}     # adjust lower bound
      upb <- 1                          # upperbound of search interval
      while(lap(upb)>0){upb<-upb+1}     # adjust upperbound
      s <- uniroot(lap,c(lwb,upb))$root       # find root of lap 
      adist$meanlog <- log(s/fitness)+adist$meanlog	# adjust scale parameter 
     
    }
  )
  adist
}

## Estimation

	      #//////////////////////////////////ML method//////////////////////////////////#
	      
# Returns the ML estimate of mean number of mutation for a sample mc, given the fitness and death
# If mfn or cvfn are non-empty, returns the estimate of mutation probability instead of the mean number, after decreasing the bias
# induced if cvfn > 0
MutationMLOptimization <- function(mc,mfn=NULL,cvfn=NULL,model=c("LD","H"),fitness=1,death=0.,winsor=512){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)

  
  # Initialization
  a.est <- MutationGFEstimation(mc,fitness=fitness,death=death,model=model)$mutations
  
  # Winsorization
  if(max(mc) > winsor) mc[mc > winsor] = winsor
  
  
  dC <- dclone(m=0:max(mc),fitness=fitness,death=death,model=model)
  
  ll <- function(a){
    p <- log(deduce.dflan(m=mc,mutations=a,fitness=fitness,death=death,model=model,clone=dC))
    -sum(p)
  }
  # gradient of log-likelihood of the sample
  dll <- function(a){
    p <- deduce.dflanda(m=mc,mutations=a,fitness=fitness,death=death,model=model,clone=dC)
    res <- p$dQ_da/p$Q
    -sum(res)
  }
  
  
  lower = 0.1*a.est
  upper = 10*a.est
  
  a.est <- lbfgsb3(prm=a.est,fn=ll,gr=dll,lower = lower,upper=upper,control=list(trace=0,iprint=-1))$prm
#   a.est <- optimx(par=a.est,fn=ll,gr=dll,lower = lower,upper=upper)$par_1
  
  dldd <- deduce.dflanda(m=mc,mutations=a.est,fitness=fitness,death=death,model=model,clone=dC)
  
  I <- sum((dldd$dQ_da)^2/(dldd$Q)^2)
  
  sda <- sqrt(1/I)
  
  if(!is.null(mfn)) {
    if(cvfn > 0){
      z4=.tunings$z4
      Mutmodel <- new(FlanMutMod,list(
		mutations=a.est,
		fitness=fitness,
		death=death,
		model=model,
		lt=TRUE
	      ))
      pm.est <- Mutmodel$unbias.mutprob(sda,z4,mfn,cvfn)
    } else pm.est <- list(mutprob=a.est/mfn,sd.mutprob=sda/mfn)
    pm.est
  } else list(mutations=a.est,sd.mutations=sda)
  
  
}

# Returns the ML estimate of mutation probability for a sample of couple (mc,fn), given the fitness and death
MutationProbabilityMLOptimization <- function(mc,fn,model=c("LD","H"),fitness=1,death=0.,winsor=512){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)
  
  mfn <- mean(fn)
  cvfn <- sd(fn)/mfn
  
  # Initialization
  pm.est <- MutationGFEstimation(mc=mc,mfn=mfn,cvfn=cvfn,fitness=fitness,death=death,model=model)$mutprob*mfn

  # Winsorization
  if(max(mc) > winsor) mc[mc > winsor] = winsor
  
  dC <- dclone(m=0:max(mc),fitness=fitness,death=death,model=model)
  
  ll <- function(pm){
    p <- mapply(function(x,y){
      y <- y/mfn
      log(deduce.dflan(m=x,mutations=pm*y,fitness=fitness,death=death,model=model,clone=dC))
    },mc,fn)
    -sum(p)
  }
  
  # gradient of log-likelihood of the sample
  dll <- function(pm){
    res <- mapply(function(x,y){
      y <- y/mfn
      p<-deduce.dflanda(m=x,mutations=pm*y,fitness=fitness,death=death,model=model,clone=dC)
      p$dQ_da*y/p$Q
    },mc,fn)
    -sum(res)
  }
  
  lower = 0.1*pm.est
  upper = 10*pm.est
  
  pm.est <- lbfgsb3(prm=pm.est,fn=ll,gr=dll,lower = lower,upper=upper,control=list(trace=0,iprint=-1))$prm/mfn
#   pm.est <- optimx(par=pm.est,fn=ll,gr=dll,lower = lower,upper=upper)$ar_1/mfn
  
  dldd <- mapply(function(x,y){
	    deduce.dflanda(m=x,mutations=pm.est*y,fitness=fitness,death=death,model=model,clone=dC)
	  },mc,fn)
  
  I <- sum((unlist(dldd[2,])*fn)^2/unlist(dldd[1,])^2)            # Fisher information
  
  sdpm <- sqrt(1/I)
  
  list(mutprob=pm.est,sd.mutprob=sdpm)

}

# Returns the ML estimates of mean number of mutation and fitness for a sample mc, given the death
# If mfn or cvfn are non-empty, returns the estimate of mutation probability instead of the mean number, after decreasing the bias
# induced if cvfn > 0
MutationFitnessMLOptimization <- function(mc,mfn=NULL,cvfn=NULL,model=c("LD","H"),death=0.,winsor=512){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)
  
  # Initialization
  est <- MutationFitnessGFEstimation(mc=mc,death=death,model=model)
  if(!est$succeeds) warning("Initialization of 'fitness' with 'GF'-method has failed.")
  a.est <- est$mutations ; r.est <- est$fitness

  # Winsorization
  if(max(mc) > winsor) mc[mc > winsor] = winsor
  
  ll <- function(param){
    a <- param[1]
    r <- param[2]
    p <- log(dflan(m=mc,mutations=a,fitness=r,death=death,model=model))
    -sum(p)
  }
  
  # gradient of log-likelihood of the sample
  dll <- function(param){
    a = param[1]
    r = param[2]
#     cat("alpha =",a,"| rho =",r,"\n")
    p <- dflan.grad(m=mc,mutations=a,fitness=r,death=death,model=model,dalpha=TRUE,drho=TRUE)
    res <- rbind(p$dQ_da,p$dQ_dr)/p$Q
    -apply(res,1,sum)
  }
  
  lower = 0.1*c(a.est,r.est)
  upper = 10*c(a.est,r.est)
  
  if(!est$succeeds) {
    lower[2] <- 10*r.est ; 
    upper[2] <- 500*r.est
  }
  
  est <- lbfgsb3(prm=c(a.est,r.est),fn=ll,gr=dll,lower = lower,upper=upper,control=list(trace=0,iprint=-1))$prm
#   est <- optimx(par=c(a.est,r.est),fn=ll,gr=dll,lower = lower,upper=upper,method="L-BFGS-B")  
  a.est = est[1]					# Update alpha estimate
  r.est = est[2]					# Update rho estimate
  
  dldd <- dflan.grad(m=mc,mutations=a.est,fitness=r.est,death=death,model=model,dalpha=TRUE,drho=TRUE)
  
  p2 <- (dldd$Q)^2
  dpa <- dldd$dQ_da
  dpr <- dldd$dQ_dr
  
  Ia <- sum((dpa^2/p2))
  Ir <- sum(dpr^2/p2)
  Iar <- sum(dpa*dpr/p2)
# #   
  det <- Ia*Ir-Iar^2
  sda <- sqrt(Ir/det)
  sdr <- sqrt(Ia/det)
  
  if(!is.null(mfn)) {
    if(cvfn > 0){
      z4=.tunings$z4
      Mutmodel <- new(FlanMutMod,list(
		mutations=a.est,
		fitness=r.est,
		death=death,
		model=model,
		lt=TRUE
	      ))
      pm.est <- Mutmodel$unbias.mutprob(sda,z4,mfn,cvfn)
    } else pm.est <- list(mutprob=a.est/mfn,sd.mutprob=sda/mfn)
    
    c(pm.est,fitness=r.est,sd.fitness=sdr)
  } else list(mutations=a.est,sd.mutations=sda,fitness=r.est,sd.fitness=sdr)
  

}


# Returns the ML estimates of mutation probability and fitness for a sample of couple (mc,fn), given the death
MutationProbabilityFitnessMLOptimization <- function(mc,fn,model=c("LD","H"),death=0.,winsor=512){
  
  if(missing(model)) model <- "LD"
  model <- match.arg(model)
  
  mfn <- mean(fn)
  cvfn <- sd(fn)/mfn
  
  
  # Initialization
  est <- MutationFitnessGFEstimation(mc,mfn=mfn,cvfn=cvfn,death=death,model=model)
  if(!est$succeeds) warning("Initialization of 'fitness' with 'GF'-method has failed.") 
  pm.est <- est$mutprob*mfn ; r.est <- est$fitness
  
  # Winsorization
  if(max(mc) > winsor) mc[mc > winsor] = winsor
  
  ll <- function(param){
    pm <- param[1]
    r <- param[2]
    p <- mapply(function(x,y){
      y <- y/mfn
      log(dflan(m=x,mutations=pm*y,fitness=r,death=death,model=model))
    },mc,fn)
    -sum(p)
  }
  
  # gradient of log-likelihood of the sample
  dll <- function(param){
    pm = param[1]
    r = param[2]
    res <- mapply(function(x,y){
      y <- y/mfn
      p<-dflan.grad(m=x,mutations=pm*y,fitness=r,death=death,model=model,dalpha=TRUE,drho=TRUE)
      rbind(p$dQ_da*y,p$dQ_dr)/p$Q
    },mc,fn)
    -apply(res,1,sum)
  }
  
  
  lower = 0.1*c(pm.est,r.est)
  upper = 10*c(pm.est,r.est)
  
  if(!est$succeeds) {
    lower[2] <- 10*r.est ; 
    upper[2] <- 500*r.est
  }
  
  est <- lbfgsb3(prm=c(pm.est,r.est),fn=ll,gr=dll,lower = lower,upper=upper,control=list(trace=0,iprint=-1))$prm
  
  pm.est = est[1]/mfn					# Update alpha estimate
  r.est = est[2]					# Update rho estimate
  
  dldd <- mapply(function(x,y){
	    dflan.grad(m=x,mutations=pm.est*y,fitness=r.est,death=death,model=model,dalpha=TRUE,drho=TRUE)
	  },mc,fn)
  
  p2 <- (unlist(dldd[1,]))^2
  dppm <- unlist(dldd[2,])
  dpr <- unlist(dldd[3,])
  
  Ipm <- sum((dppm^2/p2))
  Ir <- sum(dpr^2/p2)
  Ipmr <- sum(dppm*dpr/p2)
# #   
  det <- Ipm*Ir-Ipmr^2
  sdpm <- sqrt(Ir/det)
  sdr <- sqrt(Ipm/det)
  
  list(mutprob=pm.est,sd.mutprob=sdpm,fitness=r.est,sd.fitness=sdr)

}


	      #//////////////////////////////////P0 method//////////////////////////////////#

# Returns the P0 estimate of mean number of mutations for a sample of couple mc, given the death
MutationsNumberP0Estimation <- function(mc,death=0.){

# Return the estimate of alpha for a sample mc
# by p0-method under cell deaths with probability delta
# # when delta is known.
  if(death > 0){
    dstar <- death/(1-death)			# extinction probability
    epgf <- function(z) mean(z^mc)		# empirical PGF

    a <- -log(epgf(dstar))/(1-dstar)	#  estimate alpha

    sda <- (1-dstar)^(-2)*(epgf(dstar^2)/(epgf(dstar)^2)-1)	# variance of alpha's estimate

    sda <- sqrt(sda/length(mc))

  } else {
    p0 <- mean(mc==0)
    a <- -log(p0)
    sda <- sqrt((1/p0-1)/length(mc))
  }
  
  list(mutations=a,sd.mutations=sda)
  
}


# Returns the P0 estimate of mean number of mutations for a sample of couple mc, given the death
# If mfn or cvfn are non-empty, returns the estimate of the mutation probability instaed of the mean number
# after decreasing the induced bias if cvfn > 0
MutationP0Estimation <- function(mc,mfn=NULL,cvfn=NULL,death=0.){

# Return the estimate of alpha for a sample mc
# by p0-method under cell deaths with probability delta
# # when delta is known.
 
 
  a <- MutationsNumberP0Estimation(mc,death)
  
  sda <- a$sd.mutations
  a <- a$mutations
  
  if(!is.null(mfn)){
    pm <- a/mfn
    sdpm <- sda/mfn
    if(cvfn != 0){
      umds <- 1-death/(1-death)			# extinction probability
      temp <- a*(cvfn^2)*umds         # relative bias on mutprob
      pm <- pm*(1+temp/2)		# unbiased estimator of mutrate
      sdpm <- sdpm*(1+temp)		# standard deviation of unbiased estimator
    } 
    list(mutprob=pm,sd.mutprob=sdpm)
  } else list(mutations=a,sd.mutations=sda)  
}

# Returns the P0 estimate of mean number of mutations and fitness for a sample of couple mc,given the death
# If mfn or cvfn are non-empty, returns the estimate of the mutation probability instaed of the mean number
# after decreasing the induced bias if cvfn > 0
MutationFitnessP0Estimation <- function(mc,fn=NULL,mfn=NULL,cvfn=NULL,model=c("LD","H"),death=0.,winsor=512){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)
  
  a <- MutationsNumberP0Estimation(mc,death)
  
  if(!is.null(fn)){
    mfn <- mean(fn)
    cvfn <- sd(fn)/mfn
  }
  
  sda <- a$sd.mutations
  a <- a$mutations
  
  if(!is.null(mfn)){
    pm <- a/mfn
    sdpm <- sda/mfn
    if(cvfn != 0){
      umds <- 1-death/(1-death)			# extinction probability
      temp <- a*(cvfn^2)*umds         # relative bias on mutprob
      pm <- pm*(1+temp/2)		# unbiased estimator of mutrate
      sdpm <- sdpm*(1+temp)		# standard deviation of unbiased estimator
    }
    output <- list(mutprob=pm,sd.mutprob=sdpm)
  } else output <- list(mutations=a,sd.mutations=sda)
  
  if(!is.null(fn)) r <- FitnessP0Optimization(mc=mc,fn=fn,model=model,mut=pm,death=death,winsor=winsor)
  else r <- FitnessP0Optimization(mc=mc,fn=fn,model=model,mut=a,death=death,winsor=winsor)
    
  c(output,fitness=r$fitness,sd.fitness=r$sd.fitness)  
  
}


# Returns the ML estimate of fitness for a sample of couple mc, given the mean number of mutations (or mutation probability) and death
# If fn is non-empty, 
FitnessP0Optimization <- function(mc,fn=NULL,model=c("LD","H"),mut,death=0.,winsor=512){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)

  # Initialization
  est <- FitnessGFEstimation(mc,death,model)
  if(!est$succeeds) warning("Initialization of 'fitness' with 'GF'-method has failed.")
  r.est <- est$fitness
  
  # Winsorization
  if(max(mc) > winsor){ mc[which(mc>winsor)] = winsor}
  
  
  if(is.null(fn)){
    if(missing(mut)) mut <- 1
    
    ll <- function(r){
      p <- log(dflan(m=mc,mutations=mut,fitness=r,death=death,model=model))
      -sum(p)
    }
    
    # gradient of log-likelihood of the sample
    dll <- function(r){
      p <- dflan.grad(m=mc,mutations=mut,fitness=r,death=death,model=model,dalpha=FALSE,drho=TRUE)
      res <- p$dQ_dr/p$Q
      -sum(res)
    }
  } else {
    if(missing(mut)) mut <- 1e-9
    
    ll <- function(r){
      p <- mapply(function(x,y){
	    log(dflan(m=x,mutations=mut*y,fitness=r,death=death,model=model))
	  },mc,fn)
      -sum(p)
    }
    
    # gradient of log-likelihood of the sample
    dll <- function(r){
      res <- mapply(function(x,y){
	      p<-dflan.grad(m=x,mutations=mut*y,fitness=r,death=death,model=model,dalpha=FALSE,drho=TRUE)      
	      p$dQ_dr/p$Q
	    },mc,fn)
      -sum(res)
    }
  
  }
    
  
  lower = 0.1*r.est
  upper = 10*r.est
  
  if(est$succeeds) {
    lower <- 10*r.est ; 
    upper <- 500*r.est
  }
  
  est <- lbfgsb3(prm=r.est,fn=ll,gr=dll,lower = lower,upper=upper,control=list(trace=0,iprint=-1))
  r.est <- est$prm
  
  if(is.null(fn)) dldd <- dflan.grad(m=mc,mutations=mut,fitness=r.est,death=death,model=model,dalpha=FALSE,drho=TRUE)
  
  else dldd <- mapply(function(x,y){
		  dflan.grad(m=x,mutations=mut*y,fitness=r.est,death=death,model=model,dalpha=FALSE,drho=TRUE)
	       },mc,fn)
  
  I <- sum((dldd$dQ_dr)^2/(dldd$Q)^2)
  
  sdr <- sqrt(1/I)
  list(fitness=r.est,sd.fitness=sdr)
}

	      #//////////////////////////////////GF method//////////////////////////////////#

MutationGFEstimation <- function(mc,mfn=NULL,cvfn=NULL,fitness=1,death=0.,model=c("LD","H")){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)

  q <- .tunings$q
  b <- quantile(mc,q,names=FALSE)+1           # scaling factor
  
  Mutmodel= new(FlanMutMod,list(
    mc=mc,
    mfn=mfn,cvfn=cvfn,
    fitness=fitness,death=death,
    model=model,
    scale=b
  ))
  
  Mutmodel$MutationGFEstimation()
  
}

FitnessGFEstimation <- function(mc,death=0.,model=c("LD","H")){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)

  tu <- .tunings
  q <- tu$q
  b <- quantile(mc,q,names=FALSE)+1           # scaling factor
  z1 <- tu$z1                            # lower value
  z2 <- tu$z2                            # higher value
  z1<-z1^(1/b); z2<-z2^(1/b) ;            # rescale variables 
  
  g1 <- mean(z1^mc)                   # empirical generating function at z1
  g2 <- mean(z2^mc)                   # empirical generating function at z2
  y <- log(g1)/log(g2)                   # get ratio of logs
  
  if(model == "LD") clone=new(FlanExpClone,list(death=death))
  if(model == "H") clone=new(FlanDirClone,list(death=death))
  
  f <- function(r){
    PGF <- clone$pgf2(r,c(z1,z2))
    ((1-PGF[1])/(1-PGF[2]))-y
  }
		  
  binf <- 0.01
  bsup <- 100
  
  GFrho.succeeds <- TRUE
  if(f(binf)*f(bsup)>0) {
    rho <- 1
    GFrho.succeeds <- FALSE
  } else rho <- uniroot(f,interval=c(binf,bsup),tol=1e-6,maxiter=1e3)$root
  
  list(fitness=rho,succeeds=GFrho.succeeds)

}


MutationFitnessGFEstimation <- function(mc,mfn=NULL,cvfn=NULL,death=0.,model=c("LD","H")){

  if(missing(model)) model <- "LD"
  model <- match.arg(model)

  tu <- .tunings
  q <- tu$q
  b <- quantile(mc,q,names=FALSE)+1           # scaling factor
  z1 <- tu$z1                            # lower value
  z2 <- tu$z2                            # higher value
  z3 <- tu$z3
  z1<-z1^(1/b); z2<-z2^(1/b) ; z3 <- z3^(1/b)             # rescale variables 
  
  rho=FitnessGFEstimation(mc,death=death,model=model)
  
  if(rho$succeeds){
    
    a.est=MutationGFEstimation(mc,fitness=rho$fitness,death=death,model=model)
    
    Mutmodel <- new(FlanMutMod,list(
		  mutations=a.est$mutations,
		  fitness=rho$fitness,
		  death=death,
		  model=model,
		  lt=TRUE
		))
    
    Cov <- Mutmodel$CovGFEstimation(z1,z2,z3)
    
    sd <- sqrt(diag(Cov)/length(mc))
    
    if(!is.null(mfn)) {
      if(cvfn > 0) pm.est <- Mutmodel$unbias.mutprob(sd[1],z3,mfn,cvfn)
      else pm.est <- list(mutprob=a.est$mutations/mfn,sd.mutprob=sd[1]/mfn)
      
      c(pm.est,fitness=rho$fitness,sd.fitness=sd[2],
		succeeds=rho$succeeds)
      
    } else list(mutations=a.est$mutations,sd.mutations=sd[1],
		fitness=rho$fitness,sd.fitness=sd[2],
		succeeds=rho$succeeds)
  } else {
    res <- MutationGFEstimation(mc,mfn,cvfn,rho$fitness,death,model)
    c(res,fitness=rho$fitness,sd.fitness=0,succeeds=rho$succeeds)
  }
   
}

dclone <- function(m,fitness=1,death=0.,model=c("LD","H")){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  
  if(missing(model)){model="LD"}
  model <- match.arg(model)

  
  if(model=="LD") clone <- new(FlanExpClone,list(fitness=fitness,death=death))
  if(model == "H") clone <- new(FlanDirClone,list(fitness=fitness,death=death))
  
  
  M <- max(m)
  
  
  clone$dclone(M)[m+1]

}



dflan.grad <- function(m,mutations=1,fitness=1,death=0.,model=c("LD","H"),dalpha=TRUE,drho=FALSE){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  
  if(!(dalpha | drho)) stop("Need to precise which derivative has to be computed.")
  
  if(missing(model)){model="LD"}
  model <- match.arg(model)

#   output=.Call("dflan",
#   list(m=as.integer(max(m)),mutations=as.double(mutations),fitness=as.double(fitness),death=as.double(death),model=as.character(model),bool=as.logical(length(m)>1)))
  
  flan.mutmodel <- new(FlanMutMod,list(
    mutations=mutations,
    fitness=fitness,
    death=death,
    model=model,
    lt=TRUE
  ))
  
  M <- max(m)
  
  
  if(dalpha){
    if(drho) output <- flan.mutmodel$dflangrad(M)
    else output <- flan.mutmodel$dflanda(M)
  } else output <- flan.mutmodel$dflandr(M)
  
  lapply(output,function(l) l[m+1])
}


deduce.dflan <- function(m,mutations=1,fitness=1,death=0.,model=c("LD","H"),clone){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  if(missing(model)){model="LD"}
  model <- match.arg(model)

  
  flan.mutmodel <- new(FlanMutMod,list(
    mutations=mutations,
    fitness=fitness,
    death=death,
    model=model,
    lt=TRUE
  ))
  
  M <- max(m)
  
  output <- flan.mutmodel$deduce.dflan(M,clone)

  output[m+1]
}

deduce.dflanda <- function(m,mutations=1,fitness=1,death=0.,model=c("LD","H"),clone){

  if(sum(m < 0) > 0 | sum(trunc(m) != m) > 0){
    stop("'m' must be a vector of positive integers")
  }
  if(mutations < 0 | length(mutations) > 1){
    stop("'mutations' must be a single positive number.")
  }
  if(fitness < 0 | length(fitness) > 1){
    stop("'fitness' must be a single positive number.")
  }
  if(death < 0 | death >= 0.5 | length(death) > 1){
    stop("'death' must be a single positive and < 0.5 number.")
  }
  
  if(missing(model)){model="LD"}
  model <- match.arg(model)

#   output=.Call("dflan",
#   list(m=as.integer(max(m)),mutations=as.double(mutations),fitness=as.double(fitness),death=as.double(death),model=as.character(model),bool=as.logical(length(m)>1)))
  
  flan.mutmodel <- new(FlanMutMod,list(
    mutations=mutations,
    fitness=fitness,
    death=death,
    model=model,
    lt=TRUE
  ))
  
  M <- max(m)
  
  
  output <- flan.mutmodel$deduce.dflanda(M,clone)
  
  lapply(output,function(l) l[m+1])
}


## --------------------------------------------------------------------------------------------------------------
  ############################################## flantest object class ##############################################
  ## --------------------------------------------------------------------------------------------------------------
    
# Constructor of the class flantest
flantest <- function(Tstat,parameter,estimate,p.value,conf.int,estimates,null.value,alternative,model,method,data.name,sample.name,nsamples){
  obj <- list(Tstat=Tstat,estimate=estimate,parameter=parameter,conf.int=conf.int,p.value=p.value,
  null.value=null.value,alternative=alternative,data.name=data.name,model=model,method=method,nsamples=nsamples)
  class(obj) <- "flantest"
  obj
}

# Print function for flantest objects, inspired from print function for htest objects
print.flantest <- function(x,...){
  digits <- getOption("digits")
  prefix <- "\t"
  cat("\n")
  if(x$nsamples == 2){
    cat(strwrap(paste("Two sample ",x$method,"-test"," (",x$model," model)",sep=""), prefix=prefix), sep="\n")
  } else {
    cat(strwrap(paste("One sample ",x$method,"-test"," (",x$model," model)",sep=""), prefix=prefix), sep="\n")
  }
  cat("\n")
  cat("------------------------------------- Data -----------------------------------\n")

  if(x$nsamples == 2){
    if(length(x$data.name) == 1){
      cat("data:  ",x$data.name[[1]][1]," and ",x$data.name[[1]][2],"\n", sep="")
    }
    if(length(x$data.name) == 2){
      if(length(x$data.name[[2]]) == 2){
	cat("data:  ",x$data.name[[1]][1]," with ",x$data.name[[2]][1],
	" and ",x$data.name[[1]][2]," with ",x$data.name[[2]][2],"\n", sep="")
      } else {
	cat("data:  ",x$data.name[[1]][1]," with ",x$data.name[[2]][1],
	" and ",x$data.name[[1]][2],"\n", sep="")
      }
    }
  }
  if(x$nsamples == 1){
    if(length(x$data.name) == 1){
      cat("data:  ",x$data.name[[1]][1],"\n", sep="")
    }
    if(length(x$data.name) == 2){
      cat("data:  ",x$data.name[[1]][1]," with ",x$data.name[[2]][1],"\n", sep="")
    }
  }
  
  out <- character()

  if(!is.null(x$null.value)){
    if(length(x$null.value) == 1){
      if(!is.null(x$parameter)){
	if(x$nsamples == 2){
	  Nparam <- dim(x$param)[2]
	  cat("Sample 1 parameters : ")
	  for(i in 1:Nparam){
	    out <- c(out,paste(colnames(x$parameter)[i], "=", format(signif(x$parameter[1,i], 
	      max(1, digits - 2)))))
	  }
	  cat(strwrap(paste(out, collapse=", ")), sep="\n")
	  out <- character()
	  cat("Sample 2 parameters : ")
	  for(i in 1:Nparam){
	    out <- c(out,paste(colnames(x$parameter)[i], "=", format(signif(x$parameter[2,i], 
	      max(1, digits - 2)))))
	  }
	} else {
	  Nparam <- length(x$parameter)
	  for(i in 1:Nparam){
	    if(names(x$parameter[i])=="mfn"){
	      out <- c(out,paste(names(x$parameter[i]), "=", format(x$parameter[i], 
	      scientific=TRUE,digit=5)))
	    } else {
	      out <- c(out,paste(names(x$parameter[i]), "=", format(signif(x$parameter[i], 
	      max(1, digits - 2)))))
	    }
	  }
	}
      }
      cat(strwrap(paste(out, collapse=", ")), sep="\n")
      cat("---------------------------------- Statistics --------------------------------\n")
      if(!is.null(x$Tstat)){
	cat("Tstat =", format(signif(x$Tstat, 
	  max(1, digits - 2))),"\n")
      }
      if (!is.null(x$p.value)) {
	fp <- format.pval(x$p.value, digits= max(1, digits - 
	    3))
	cat("p-value for",names(x$null.value[1]), if (substr(fp, 1, 1) == 
	  "<") fp else paste("=", fp),"\n")
      }
      
      if (!is.null(x$alternative)) {
	cat("Alternative hypothesis: ")
	alt.char <- switch(x$alternative, two.sided="not equal to", 
		less="less than", greater="greater than")
	cat("true ", names(x$null.value), " is ", alt.char, 
	  " ", x$null.value, "\n", sep="")
      }
      
      if (!is.null(x$conf.int)) {
	cat(format(100 * attr(x$conf.int, "conf.level")), " percent confidence interval:\n",
	paste(format(x$conf.int),collapse="  "),
	"\n")
      }
      if (!is.null(x$estimate)) {
	if(x$nsamples == 1){
	  cat("Sample estimates : \n")
	  print(x$estimate)
	}
	if(x$nsamples == 2){
	  cat("Sample 1 estimates : \n")
	  print(x$estimate[1,])
	  cat("Sample 2 estimates : \n")
	  print(x$estimate[2,])
	}
      }
    } else if(length(x$null.value)==2){
      
      if(!is.null(x$parameter)){
	if(x$nsamples == 2){
	  cat("Sample 1 parameters : ")
	  Nparam <- dim(x$param)[2]
	  for(i in 1:Nparam){
	    out <- c(out,paste(colnames(x$parameter)[i], "=", format(signif(x$parameter[1,i], 
	      max(1, digits - 2)))))
	  }
	  cat(strwrap(paste(out, collapse=", ")), sep="\n")
	  out <- character()
	  cat("Sample 2 parameters : ")
	  for(i in 1:Nparam){
	    if(colnames(x$parameter)[i] == "mfn"){
	      out <- c(out,paste(colnames(x$parameter)[i], "=", format(x$parameter[2,i], 
	      scientific=TRUE,digit=5)))
	    } else {
	      out <- c(out,paste(colnames(x$parameter)[i], "=", format(signif(x$parameter[2,i], 
	      max(1, digits - 2)))))
	    }
	  }
	} 
	if(x$nsamples == 1){
	  cat("Sample parameters : ")
	  Nparam <- length(x$parameter)

	  for(i in 1:Nparam){
	    out <- c(out,paste(names(x$parameter[i]), "=", format(signif(x$parameter[i], 
	      max(1, digits - 2)))))
	  }
	}
      }
      cat(strwrap(paste(out, collapse=", ")), sep="\n")
      cat("---------------------------------- Statistics --------------------------------\n")
      if(!is.null(x$Tstat)){
	cat("Tstat = (",format(signif(x$Tstat[1], 
	  max(1, digits - 2)))," , ",format(signif(x$Tstat[2],
	  max(1, digits - 2))),")","\n",sep="")
      }
      if (!is.null(x$p.value)) {
	fp1 <- format.pval(x$p.value[1], digits=max(1, digits - 
	    3))
	fp2 <- format.pval(x$p.value[2], digits=max(1, digits - 
	    3))
	cat("p-value for",names(x$null.value[1]), if (substr(fp1, 1, 1) == 
	  "<") fp1 else paste("=", fp1),"\np-value for",names(x$null.value[2]),if (substr(fp2, 1, 1) == 
	  "<") fp2 else paste("=", fp2),"\n")
      }
      if (!is.null(x$alternative)) {
	cat("Alternative hypotheses: ")
	alt.char1 <- switch(x$alternative[1], two.sided="not equal to", 
		less="less than", greater="greater than")
	alt.char2 <- switch(x$alternative[2], two.sided="not equal to", 
		less="less than", greater="greater than")
	cat("true ", names(x$null.value[1]), " is ", alt.char1, 
	  " ", x$null.value[1], "\n", sep="")
	cat("                        ")
	cat("true ", names(x$null.value[2]), " is ", alt.char2, 
	  " ", x$null.value[2], "\n", sep="")
      }
      if (!is.null(x$conf.int)) {
	cat(format(100 * attr(x$conf.int, "conf.level")), " percent confidence interval for ",names(x$null.value[1]),": \n", 
	  paste(format(x$conf.int[,1]),collapse="   "),
	  "\n",
	  format(100 * attr(x$conf.int, "conf.level")), " percent confidence interval for ",names(x$null.value[2]),": \n",
	  paste(format(x$conf.int[,2]),collapse="   "),
	  "\n",sep="")
      }
      if (!is.null(x$estimate)) {
	if(x$nsamples == 1){
	  cat("Sample estimates : \n")
	  print(x$estimate[1,])
	}
	if(x$nsamples == 2){
	  cat("Sample 1 estimates : \n")
	  print(x$estimate[1,])
	  cat("Sample 2 estimates : \n")
	  print(x$estimate[2,])
	}
      }
    } 
    cat("\n")
  }    
}


draw.clone <- function(t,mutprob=1e-2,fitness=1,death=0.,
	     dist=list("lnorm",meanlog=-0.3795851,sdlog=0.3016223)){
  #   Simulates a clone up to time t, represents the clone, and
  #   returns the vector of split times.
  #   At each division, the probability of having a mutant is mutprob.
  #   Variable death is the probability of death of a cell. 
  #   The division times of have distribution dist, normalized
  #   to unit growth rate for mutant, to fitness for normal cells.
  #   The distribution dist is given as a list of a character chain 
  #   followed by the list of parameters. 
  #
  #   usage: draw.clone(5)
  #          draw.clone(t=9,fitness=0.5,death=0.1)
  #          draw.clone(t=3,mutprob=0.1,fitness=2,dist=dexp,death=0.2)
  #
  
  
    names(dist)[1] <- "name"
    if(dist$name == "exp"){
      names(dist)[2] <- "rate"
    } else if(dist$name == "dirac"){
      names(dist)[2] <- "location"
    } else if(dist$name == "lnorm"){
      names(dist)[2] <- "meanlog"
      names(dist)[3] <- "sdlog"
    } else if(dist$name == "gamma"){
      names(dist)[2] <- "shape"
      names(dist)[3] <- "scale"
    } else stop("'dist[[1]]' must be a character chain among 'exp', 'dirac', 'lnorm', 'gamma'")

    
		    # initialization
  p <- mutprob                            # mutation probability
  distn <- adjust.rate(dist,fitness=fitness,death=death)      # division times of normal cells
  distm <- adjust.rate(dist,fitness=1,death=death)            # division times of mutants
  g <- 0                                  # current generation
  ce <- 1                                 # cells of current generation
  mu <- FALSE                                 # mutants in current generation
  bd <- 0                                 # birth dates of current generation
  dt <- rdt(1,distn)                      # life times of current generation
  de <- bd+dt                             # death dates of current generation
  or <- 1                                 # location of cells in current generation
  div <- {which((de<t)&                   # indices of cells dead before t
	      (runif(length(de))>death))} # that will divide
  de <- min(de,t)                         # truncate death dates
  nd <- length(div)                       # number of dividing cells
  nc <- nd                                # total number of cells
  CE <- ce                                # all cells
  MU <- mu                                # all mutants
  BD <- bd                                # all birth dates
  DE <- de                                # all death dates
  OR <- or                                # all locations
		    # main loop
  while((nd>0)&&(nc<1e+4)){               # while divisions still happen  
      g <- g+1                            # next generation
      ce <- 2*ce[div]; ce <- c(ce,ce+1)   # cells in next generation
	  mu <- mu[div]                   # mutants that divide
	  mu1 <- {sample(c(TRUE,FALSE),length(div),                      
	      replace=TRUE,prob=c(p,1-p))}# new mutants in next generation
  # 	if (){mu1}
	  mu1 <- mu|mu1                   # mutants in next generation
	  mu <- c(mu1,mu)                 # all mutants in next generation
      nd <- nd*2                          # double the number of cells
      bd <- c(de[div],de[div])            # birth dates of daughters
	  imu <- which(mu)                # indices of mutants
	  nmu <- length(imu)              # number of mutants
	  ino <- which(!mu)               # indices of normal cells
	  nno <- length(ino)              # number of normal cells
      dtmu <- rdt(nmu,distm)              # division times of mutants
      dtno <- rdt(nno,distn)              # division times of normal cells
	  dt <- rep(0,nd)                 # division times of daughters
      dt[imu] <- dtmu                     # insert division times of mutants
	  dt[ino] <- dtno                 # insert division times of normal cells
      de <- bd+dt                         # death dates of daughters
      or<-c(or[div]-2^(-g),or[div]+2^(-g))# location of daughters 
      div <- {which((de<t)&               # indices of cells dead before t
	      (runif(length(de))>death))} # that will divide
      de <- mapply(min,de,rep(t,nd))      # truncate death dates
      nd <- length(div)                   # number of dividing daughters
      nc <- nc+nd                         # total number of cells
      CE <- c(CE,ce)                      # stack new cells
      MU <- c(MU,mu)                      # stack mutants
      BD <- c(BD,bd)                      # stack birth dates
      DE <- c(DE,de)                      # stack death dates
      OR <- c(OR,or)                      # stack locations
  }                                       # end while
  if (nc > 1e+4){
    warning("too many cells to plot")
    }else{
		    # treatment for graphics
  if(mutprob == 0){colnor <- "green4"}
  else{
  colnor <- "green4"                      # color for normal cells
  colmut <- "orangered"                   # color for mutants
  }
  nc <- length(CE)                        # number of cells
  if (nc>1){                              # something to plot
  OR <- sort(OR,index.return=TRUE)$ix     # order for plotting
  CE <- CE[OR]                            # reorder cells
  MU <- MU[OR]                            # reorder mutants
  BD <- BD[OR]                            # reorder births
  DE <- DE[OR]                            # reorder deaths
  nmu <- length(which(MU))                # number of mutants
  nno <- nc - nmu                         # number of normal
  BDN <- BD[!MU]                          # births of normal
  DEN <- DE[!MU]                          # deaths of normal
  XN <- rbind(which(!MU),which(!MU))      # abscissas for vertical lines normal
  YN <- rbind(BDN,DEN)                    # ordinates for vertical lines normal
  nlc <- nc - length(which(DE<t))         # living cells at time t
  nlm <- nmu - length(which((DE<t)&MU))   # living mutants at time t
  para <- {substitute(list(               # list of parameters
	  "mutation probability" == t1, fitness == t2, death==t3,
	  cells == t4, mutants ==t5),    # names
	  list(t1=p,t2=fitness,t3=death,
	  t4=nlc,t5=nlm))}               # values
#   if(!Lab){para = NULL}
  name <- switch(dist[[1]],               # distribution of division times
	  dirac="constant lifetimes",
	  exp="exponential lifetimes",
	  gamma="gamma lifetimes",
	  lnorm="log-normal lifetimes")
  {matplot(XN,YN,                         # plot vertical lines
	  main=para,                     # add title
	  xlim=c(1,nc),             
	  ylim = c(t,0),                 # reverse y axis
	  type="l",col=colnor,           # color for normal
	  lty=rep(1,nno),                # solid lines
	  axes=FALSE,                    # remove axes
	  xlab=name,ylab="time")}        # axes labels
  axis(side=2)                            # draw time axis
  if(nmu>0){                              # if any mutant
    BDM <- BD[MU]                        # births of mutant
    DEM <- DE[MU]                        # deaths of mutant
    XM <- rbind(which(MU),which(MU))     # abscissas for vertical lines mutant
    YM <- rbind(BDM,DEM)                 # ordinates for vertical lines mutant
    {matlines(XM,YM,                     # plot horizontal lines
	  type="l",col=colmut,           # color for mutant
	  lty=rep(1,nmu))}               # solid lines
  }                                       # end if any mutants
  div <- which((DE<t)&(!MU))              # indices dividing cells normal
  DIN <- CE[div]                          # identities dividing cells normal
  DLN <- 2*DIN                            # identities of left daughters
  IO <- sort(CE,index.return=TRUE)        # numeral order of cells
  IOx <- IO$x; IOix <- IO$ix              # get cells and order
  indn <- which(IOx %in% DLN)             # identify left daughters normal 
  XLN <- IOix[indn]                       # indices of left daughters normal
  XRN <- IOix[indn+1]                     # indices of right daughters normal
  XN <- rbind(XLN,XRN)                    # abscissas for horizontal lines
  YN <- rbind(BD[XLN],BD[XRN])            # ordinates for horizontal lines
  {matlines(XN,YN,                        # plot horizontal lines
	  type="l",col=colnor,           # color normal
	  lty=rep(1,nno))}               # solid lines
  div <- which((DE<t)&MU)                 # indices dividing cells mutant
  if (length(div>0)){                     # if any dividing mutant
    DIM <- CE[div]                       # identities dividing cells mutant
    DLM <- 2*DIM                         # identities of left daughters
    indm <- which(IOx %in% DLM)          # identify left daughters mutant
    XLM <- IOix[indm]                    # indices of left daughters mutant
    XRM <- IOix[indm+1]                  # indices of right daughters mutant
    XM <- rbind(XLM,XRM)                 # abscissas for horizontal lines
    YM <- rbind(BD[XLM],BD[XRM])         # ordinates for horizontal lines
    {matlines(XM,YM,                     # plot horizontal lines
	  type="l",col=colmut,           # color mutant
	  lty=rep(1,nmu))}               # solid lines
  }                                       # end if any dividing mutant
  }                                       # end if something to plot
  }                                       # end if graphics
#   ST <- sort(BD)                          # split times
#   ST <- append(ST,t)                      # final instant
#   list(ST=ST,nc=nc)                              # return birth dates
}


rdt <- function(n,dist){
#   Returns a sample of size n of distribution dist.
#   The distribution dist is given as a list of a character chain 
#   followed by the list of parameters. 
#
#   usage: dist <- list("lnorm",meanlog=0,sdlog=1)
#          rdt(100,dist)
#          rdt(1000,goflnormBA)
#
switch(dist[[1]],
    dirac={                             # constant division times
        l <- dist$location              # location parameter
        rep(l,n)
          },                            # end Dirac distribution
    exp={                               # exponential distribution
        r <- dist$rate                  # rate parameter
        rexp(n,rate=r)
          },                            # end exponential distribution
    gamma={                             # gamma distribution
        a <- dist$shape                 # shape parameter
        l <- dist$scale                 # scale parameter
        rgamma(n,shape=a,scale=l)
          },                            # end gamma distribution
    lnorm={                             # Log Normal distribution
        a <- dist$sdlog                 # sdlog parameter
        l <- dist$meanlog               # meanlog parameter
        rlnorm(n,meanlog=l,sdlog=a)
          }                             # end log normal distribution
)                                       # end switch
}                                       # end function rdt