# INTERNAL
Nonpareil.__init_globals <- function(
    ### Internal Nonpareil function.
    soft=TRUE){
  if(soft && exists('Nonpareil.LastFactor')) return();
  Nonpareil.LastFactor <<- NULL;
  Nonpareil.LastXmax <<- NULL;
  Nonpareil.OpenColors <<- c();
  Nonpareil.OpenNames <<- c();
}

# MODEL FUNCTIONS
#Nonpareil.f <- function(x,a,b){ return( ( 1 - exp(-x/(2^a)) )^b ) }
#Nonpareil.antif <- function(y,a,b){ return( -(2^a)*log(1 - y^(1/b)) ) }
Nonpareil.f <- function(
    ### Function of the projected model.
    x,
    ### Values of sequencing effort (typically in bp).
    a,
    ### Parameter alpha of the gamma CDF.
    b
    ### Parameter beta of the Gamma CDF.
    ){
  return(pgamma(log(x+1),a, b))
  ### Predicted values of abundance-weighted average coverage.
}
Nonpareil.antif <- function(
    ### Complement function of Nonpareil.f.
    y,
    ### Values of abundance-weighted average coverage
    a,
    ### Parameter alpha of the gamma CDF.
    b
    ### Parameter beta of the gamma CDF.
    ){
  return(exp(qgamma(y,a,b))-1)
  ### Estimated sequencing effort.
}

# NONPAREIL CURVES
Nonpareil.curve.batch <- function(
    ### Generates a single plot (and a single table of estimates) with multiple
    ### Nonpareil curves, from different .npo files.
    files,
    ### Vector of characters with the paths to the .npo files.
    overlap,
    ### Value of the '-L' parameter (in Nonpareil, the default is 50). It can be
    ### a number (if all the curves were generated with the same value) or a
    ### vector (in the same order of 'files'). See the 'overlap' value of
    ### 'Nonpareil.curve()'. Use only with Nonpareil < v2.0.
    r=NA,
    ### Values of 'r' in 'Nonpareil.curve()' (red value).
    g=NA,
    ### Values of 'g' in 'Nonpareil.curve()' (green value).
    b=NA,
    ### Values of 'b' in 'Nonpareil.curve()' (blue value).
    libnames=NA,
    ### Vector of names of the libraries (corresponding to 'libname' in
    ### 'Nonpareil.curve()'). It must be characters, not factors.
    read.lengths=NA,
    ### A vector of numbers indicating the length of the reads (corresponding to
    ### 'read.length' in 'Nonpareil.curve()'). Use only with Nonpareil < v2.0.
    col=NA,
    ### Values of 'col' in 'Nonpareil.curve()'.
    ...
    ### Any other parameter accepted by 'Nonpareil.curve()' is supported.
    ){
  if(!is.vector(files)) files = as.vector(files);
  if(missing(overlap)) overlap=rep(NULL, length(files));
  new=TRUE;
  for(i in 1:length(files)){
    o = Nonpareil.curve(as.character(files[i]), overlap, r=r[i], g=g[i], b=b[i],
            col=col[i], new=new,
            libname=ifelse(is.na(libnames[1]), NA, as.character(libnames[i])),
            read.length=ifelse(is.na(read.lengths[1]), NA, read.lengths[i]),
            ...);
    if(new){
      out.m = matrix(NA, ncol=length(o), nrow=length(files));
      colnames(out.m) <- rownames(as.matrix(o));
      if(!is.na(libnames[1])) rownames(out.m) <- libnames;
    }
    out.m[i, ] <- as.numeric(o);
    if(!is.null(Nonpareil.LastXmax)) new <- FALSE;
  }
  return(out.m);
  ### A dataframe containing the values generated by 'Nonpareil.curve()'.
}

Nonpareil.curve <- function(
    ### Generates a Nonpareil curve from a .npo file.
    file,
    ### Path to the .npo file, containing the read redundancy.
    ksize=NULL,
    ### K-mer size (if -T kmer). By default this value is extracted
    ### from the .npo file headers.
    overlap=NULL,
    ### Value of the '-L' parameter (in Nonpareil, the default is 50). If not
    ### set, it tries to find the value in the .npo file (supported in Nonpareil
    ### >= 2.0), or fails with an error message. Use only with Nonpareil < v2.0.
    factor=1,
    ### Multiplier of the sequencing effort. This can be used to express the
    ### sequencing effort in units other than base pairs (bp). For example, to
    ### express sequencing effort as Gbp, use factor=1e-9. This can also affect
    ### the fit of the model, and it's considered EXPERIMENTAL.
    plotDispersion=NA,
    ### Indicates if (and how) dispersion of the replicates should be plotted.
    ### It requires modelOnly=FALSE to take effect. It can be NA, in which case
    ### no dispersion is plotted, or any of the following strings: 'sd' (one
    ### standard deviation around the mean), 'ci95' (95% confidence interval),
    ### 'ci90' (90% confidence interval), 'ci50' (50% confidence interval), 'iq'
    ### (inter-quartile range).
    returnModelValues=FALSE,
    ### If TRUE, returns the coordinates of the model as model.x and model.y.
    returnModelParameters=FALSE,
    ### If TRUE, returns the model itself as model.
    xmax=10e12,
    ### Maximum sequencing effort to plot.
    xmin=1e3,
    ### Minimum sequencing effort to plot
    ymax=1,
    ### Maximum coverage to plot.
    ymin=1e-6,
    ### Minimum coverage to plot.
    xlab=NULL,
    ### Label of the X-axis. If NULL, it's set to sequencing effort and the
    ### units (see factor).
    ylab=NULL,
    ### Label of the Y-axis. If NULL, it's set to Estimated average coverage.
    col=NA,
    ### The color of the curve. If passed, it overrides `r`, `g`, and `b`.
    r=NA,
    ### Red component of the curve's color. If NA, it's randomly set. If <=1,
    ### it's assumed to be in the range [0,1]; if >1, it's assumed to be in the
    ### range [0,256].
    g=NA,
    ### Green component of the curve's color. Same as `r`.
    b=NA,
    ### Blue component of the curve's color. Same as `r`.
    new=TRUE,
    ### If FALSE, it attempts to use a previous (active) canvas to plot the
    ### curve.
    plot=TRUE,
    ### Determines if the plot should be produced. If FALSE, it still computes
    ### the coverage and the model.
    libname=NA,
    ### Name of the library. If NA, it's determined by the file name. This is
    ### useful if you plan to call Nonpareil.legend().
    modelOnly=FALSE,
    ### If TRUE, the rarefied data is not presented, only the fitted model.
    plotModel=TRUE,
    ### If FALSE, the model is not plotted (but it's still computed).
    plotDiversity=FALSE,
    ### If TRUE, the diversity estimate is plotted as a small arrow below the
    ### Nonpareil curve.
    curve.lwd=2,
    ### Line width of the Nonpareil curve.
    curve.alpha=0.4,
    ### Alpha value (from 0 to 1) of the Nonpareil curve.
    model.lwd=1,
    ### Line width of the model.
    model.alpha=1,
    ### Alpha value (from 0 to 1) of the model.
    log='x',
    ### Axis to plot in logarithmic scale. It can be 'x' (sequencing effort,
    ### default), 'y' (coverage), 'xy' (both logarithmic), or '' (both linear).
    data.consistency=TRUE,
    ### If TRUE, it checks the consistency of the data before plotting.
    useValue='mean',
    ### Controls how the replicates are to be summarized at each point of
    ### sequencing effort. It can be any of: 'mean' (average of the replicates),
    ### 'median' (median), 'ub' (upper boundary of the 95% confidence interval),
    ### 'lb' (lower boundary of the 95% confidence interval), 'q1' (quartile 1),
    ### 'q3' (quartile 3). Note that the quartile 2 is also 'median'.
    star=95,
    ### Objective coverage (in percentage). By default: 95, which means the
    ### sequencing effort required to reach 95% average coverage is to be
    ### estimated.
    read.length=NA,
    ### Length of the reads. Use only with Nonpareil < v2.0.
    weights.exp=NA,
    ### Vector of values to be tested (in order) as exponent of the weights
    ### distribution. If the model fails to converge, sometimes manual
    ### modifications in this parameter may help. By default (NA), five
    ### different values are tested in the following order: For linear sampling,
    ### -1.1, -1.2, -0.9, -1.3, -1. For logarithmic sampling (-d option in
    ### Nonpareil), 0, 1, -1, 1.3, -1.1, 1.5, -1.5.
    correction.factor=TRUE,
    ### Should the overlap-dependent correction factor be applied? If FALSE,
    ### redundancy is assumed to equal coverage.
    ...
    ### Any other parameters accepted by plot().
    ){
  # Create environment
  Nonpareil.__init_globals(!new);

  # Examine consistency
  if(is.null(file)) stop("The file argument is mandatory");
  if(!new && (is.null(Nonpareil.LastFactor) || is.null(Nonpareil.LastXmax)))
    stop("No previous plot found, use new=TRUE if this is the first plot.");

  # Read metadata (.npo header)
  log_sampling <- 0;
  meta_data <- gsub('^# @', "", grep("^# @", readLines(file), value=TRUE));
  keys <- gsub(': .*', "", meta_data);
  vals <- gsub('.*: ', "", meta_data);
  type <- "alignment";
  if("ksize" %in% keys) type <- "kmer";
  if(type=="kmer"){
    if(is.null(ksize) & "ksize" %in% keys)
      ksize = as.numeric(vals[keys=="ksize"]);
    if(is.na(read.length) & "AL" %in% keys)
      read.length = as.numeric(vals[keys=="AL"]);
    overlap <- 0;
  }else{
    if(is.null(overlap) & "overlap" %in% keys)
      overlap = as.numeric(vals[keys=="overlap"]);
    if(is.null(overlap))
      stop("The overlap argument is mandatory for Nonpareil v < 2.0.");
    if(is.na(read.length) & "L" %in% keys)
      read.length = as.numeric(vals[keys=="L"]);
    ksize <- 1;
  }
  if(is.na(read.length)) read.length=101-ksize+1;
  num_reads = NULL;
  if("R" %in% keys) num_reads = as.numeric(vals[keys=="R"]);
  if("divide" %in% keys) log_sampling = as.numeric(vals[keys=="divide"]);
  if("logsampling" %in% keys)
    log_sampling = as.numeric(vals[keys=="logsampling"]);
  
  # Read input
  out <- list(kappa=0, C=0, LRstar=0, LR=0, modelR=0, diversity=0);
  a <- read.table(file, sep="\t", h=F);
  out$kappa <- a$V2[nrow(a)];
  if(is.null(num_reads)) num_reads = max(a$V1);
  LR <- exp(log(num_reads) + log(read.length));
  out$LR <- LR;
  cor.f <- 1.0;
  if(correction.factor) cor.f <- Nonpareil.coverageFactor(overlap)
  for(i in 2:6) a[, i] <- a[, i]^cor.f;
  if(useValue=="median"){
    values = a[, 5]
  }else if(useValue=="ub"){
    values = pmin(1, a[, 2] + a[, 3]*1.9)
  }else if(useValue=="lb"){
    values = pmax(0, a[, 2] - a[, 3]*1.9)
  }else if(useValue=="q1"){
    values = a[, 4]
  }else if(useValue=="q3"){
    values = a[, 6]
  }else{
    values = a[, 2]
  }
  a$V1 = exp(max(log(a$V1)) + max(values^0.27)*(log(a$V1) - max(log(a$V1))));
  a$V1 = exp(log(a$V1)*0.61 + 10);
  a$V1 = a$V1 * read.length / 101;
  horiz.diff = LR / max(a$V1);
  a$V1 = a$V1 * horiz.diff;##
  horiz.diff = 1;
  if(is.na(libname)) {
    libname <- basename(file);
    if(substr(libname, nchar(libname)-3, nchar(libname))==".npo")
      libname <- substr(libname, 0, nchar(libname)-4);
  }

  # Check data
  out$C <- max(values);
  if(data.consistency){
    twenty.pc = which.max(a$V1[a$V1<=0.5*max(a$V1)]);
    if(a[twenty.pc, 5]==0){
      warning(paste("Median of the curve is zero at 50% of the reads, check",
        "parameters and re-run (e.g., decrease value of -L in nonpareil)."));
      return(out);
    }
    if(a[nrow(a), 2]<=1e-5){
      warning(paste("Curve too low, check parameters and re-run (e.g.,",
        "decrease value of -L in nonpareil)."));
      return(out);
    }
    if(a[2, 2]>=1-1e-5){
      warning(paste("Curve too steep, check parameters and re-run (e.g.,",
        "decrease value of -i in nonpareil or use -d)."));
      return(out);
    }
  }

  # For future calls (of this or other Nonpareil.* functions)
  Nonpareil.LastFactor <<- factor;

  # Draw it
  if(plot){
    if(new){
      if(!is.null(xlab)){     xlab <- as.character(xlab) }
      else if(factor==1){     xlab <- 'Sequencing effort (bp)' }
      else if(factor==1e-1){  xlab <- 'Sequencing effort (Dbp)' }
      else if(factor==1e-2){  xlab <- 'Sequencing effort (Hbp)' }
      else if(factor==1e-3){  xlab <- 'Sequencing effort (Kbp)' }
      else if(factor==1e-6){  xlab <- 'Sequencing effort (Mbp)' }
      else if(factor==1e-9){  xlab <- 'Sequencing effort (Gbp)' }
      else if(factor==1e-12){ xlab <- 'Sequencing effort (Tbp)' }
      else if(factor==1e-15){ xlab <- 'Sequencing effort (Pbp)' }
      else if(factor==1e-18){ xlab <- 'Sequencing effort (Ebp)' }
      else if(factor==1e-21){ xlab <- 'Sequencing effort (Zbp)' }
      else{ xlab <- paste('Sequencing effort (by ', factor, 'bp)', sep='') }
      if(is.null(ylab))       ylab <- 'Estimated average coverage';

      plot(1, t='n', xlim=c(xmin, xmax), ylim=c(ymin, ymax), bty="l",
        xlab=xlab, ylab=ylab, xaxs="i", yaxs="i", log=log, ...);
      abline(h=c(1,star/100), lty=2, col='red');
      abline(v=10^seq(0,15,by=3), lty=2, col='gray80')
    }

    if(!is.na(col)){
      r <- col2rgb(col)['red',1]/256;
      g <- col2rgb(col)['green',1]/256;
      b <- col2rgb(col)['blue',1]/256;
    }
    if(is.na(r)) r <- sample(200,1);
    if(is.na(g)) g <- sample(200,1);
    if(is.na(b)) b <- sample(200,1);
    if(r>1) r <- r/256;
    if(g>1) g <- g/256;
    if(b>1) b <- b/256;

    if(!modelOnly){
      if(!is.na(plotDispersion)){
        if(plotDispersion == 'sd'){
          err.y <- c(values+a$V3, rev(values-a$V3));
        }else if(plotDispersion == 'ci95'){
          err.y <- c(values+a$V3*1.9, rev(values-a$V3*1.9));
        }else if(plotDispersion == 'ci90'){
          err.y <- c(values+a$V3*1.64, rev(values-a$V3*1.64));
        }else if(plotDispersion == 'ci50'){
          err.y <- c(values+a$V3*.67, rev(values-a$V3*.67));
        }else if(plotDispersion == 'iq'){
          err.y <- c(a$V4, rev(a$V6));
        }
        polygon(horiz.diff*c(a$V1, rev(a$V1))*factor, border=NA,
          ifelse(err.y<=ymin*0.1, ymin*0.1, err.y), col=rgb(r,g,b,.2));
      }
      lines(horiz.diff*a$V1*factor, values, col=rgb(r,g,b,curve.alpha),
        lwd=curve.lwd);
    }

    # Save some info
    Nonpareil.LastXmax   <<- xmax;
    Nonpareil.OpenColors <<- c(Nonpareil.OpenColors, rgb(r,g,b));
    Nonpareil.OpenNames  <<- c(Nonpareil.OpenNames, libname);
  }

  # Model it
  sel  <- values>0 & values<0.9;
  x <- a$V1[sel];
  if((length(x)>10) | ! data.consistency){
    y <- values[sel];
    data <- list(x=x, y=y)
    if(is.na(weights.exp[1])){
      if(log_sampling==0){ weights.exp <- c(-1.1,-1.2,-0.9,-1.3,-1) }
      else{ weights.exp <- c(0, 1, -1, 1.3, -1.1, 1.5, -1.5) }
    }
    estModel <- TRUE;
    weights.i <- 0;
    while(estModel & !is.na(weights.exp[weights.i+1])){
      weights.i <- weights.i+1;
      suppressWarnings(model <- nls(y ~ Nonpareil.f(x, a, b), data=data,
        weights=(a$V3[sel]^weights.exp[weights.i]), start=list(a=1, b=0.1),
        lower=c(a=0, b=0), algorithm='port', control=nls.control(
          minFactor=1e-25000, tol=1e-15, maxiter=1024, warnOnly=T)));
      if(summary(model)$convInfo$isConv) estModel <- FALSE;
    }
    if(summary(model)$convInfo$isConv){
      model.lty=2;
      if(modelOnly) model.lty=1;
      if(plot & plotModel){
        model.x <- exp(seq(log(xmin), log(xmax), length.out=1e3));
        model.y <- predict(model, list(x=model.x));
        model.x <- model.x*horiz.diff*factor;
        lines(model.x, model.y, col=rgb(r,g,b,model.alpha), lty=model.lty,
          lwd=model.lwd);
        if(modelOnly)
          points(max(a$V1), predict(model, list(x=max(a$V1))), col=rgb(r,g,b),
            pch=21, bg='white');
      }
      if(returnModelValues) { out$model.x <- model.x ; out$model.y <- model.y; }
      if(returnModelParameters) { out$model=model }
      pa <- coef(model)[1];
      pb <- coef(model)[2];
      if(pa > 1) out$diversity <- (pa-1)/pb;
      if(plot & plotDiversity & out$diversity>0)
        arrows(x0=exp(out$diversity),
          y0=ifelse(log=='y' | log=='xy' | log=='yx', ymin*0.1, ymin-1),
          y1=ymin, length=0.1, col=rgb(r,g,b,model.alpha));
      out$LRstar <- Nonpareil.antif(star/100, pa, pb);
      out$modelR <- cor(y, predict(model, list(x=x)));
    }else{
      warning("Model didn't converge.");
    }
  }else{
    warning("Insufficient resolution below 90% coverage, skipping model");
  }
  return(out);
  ### A list with the following entries:
  ###
  ### kappa: "Redundancy" value of the entire dataset.
  ###
  ### C: Average coverage of the entire dataset.
  ###
  ### LRstar: Estimated sequencing effort required to reach the objective
  ###   average coverage (star, 95% by default).
  ###
  ### LR: Actual sequencing effort of the dataset.
  ###
  ### modelR: Pearson's R coefficient betweeen the rarefied data and the
  ###   projected model.
  ###
  ### diversity: Nonpareil-based diversity index. This value's units are the
  ###   natural logarithm of the units of sequencing effort (by default,
  ###   log-bp), and indicates the inflection point of the fitted model for the
  ###   Nonpareil curve. If the fit doesn't converge, or the model is not
  ###   estimated, the value is zero (0).
  ###
  ### model.x (if retturnModelValues=TRUE): Values of sequencing effort in the
  ###   projected model.
  ###
  ### model.y (if retturnModelValues=TRUE): Values of average coverage in the
  ###   projected model.
  ###
  ### model (if returnModelParameters=TRUE): Fitted model.
  ###
}

Nonpareil.legend <- function(
    ### Generates a legend for Nonpareil plots, after calling Nonpareil.curve()
    ### or Nonpareil.curve.batch().
    x=NULL,
    ### X coordinate, or any character string accepted by legend (e.g.,
    ### 'bottomright').
    y=.3,
    ### Y coordinate.
    ...
    ### Any other parameters supported by legend().
    ){
  if(is.null(Nonpareil.LastFactor) || is.null(Nonpareil.LastXmax))
    stop(paste("There must be at least one active plot to draw a legend. Use",
      "Nonpareil.curve()."));
  if(is.null(x)) x <- 0.75*Nonpareil.LastXmax;
  legend(x=x, y=y, legend=Nonpareil.OpenNames, fill=Nonpareil.OpenColors, ...);
}

## PREDICTIONS
Nonpareil.coverageFactor <- function(
    ### Factor to transform redundancy into coverage (internal function).
    overlap
    ### Value of overlap (-L).
    ){
  if(overlap==0) overlap <- 50 # for type kmer
  return(1-exp(2.23E-2*overlap-3.5698))

  ### A numeric scalar.
}