#This file contains the generator functions for recruiting extra subjects 

#into a predict from data simulation. It also contains the code for

#estimating k in the accrual cdf G(t) t^k/B^k from recruitment data.

##' @include eventData.R fromDataInternal.R

NULL

##' A class to allow additional subjects 

##' to be recruited onto the study

##' @slot f A function with 1 argument (the number of subjects to be recruited)

##' which returns a vector of recruitment times (of the required length)

##' @slot model A string, name of the accrual procedure used  

##' @slot text The summary text to be output when printing the accrual procedure

##' for example "a poisson process with rate..."

##' @seealso \code{\link{show,AccrualGenerator-method}}

##' @export

setClass("AccrualGenerator",

         slots=list(f="function",

                    model="character",

                    text="character"

         )

)

##' @name show

##' @rdname show-methods

##' @aliases show,AccrualGenerator-method

##' @export

setMethod("show","AccrualGenerator",function(object){

  cat(paste("Accrual model:",object@model,"\n"))

  cat(paste("Model Description:\n"))

  cat(object@text)  

})

##' Create A Poisson process accrual generator

##' @param start.date The start date for the recruitment process

##' @param rate The rate of the required Poisson process.

##' @return An \code{AccrualGenerator} object which can 

##' be used when recruiting subjects with a Possion process

##' with given \code{start.date} and \code{rate}

##' @export 

Generate.PoissonAccrual <- function(start.date,rate){

  if(rate <= 0) stop("rate must be positive")

  start.date <- FixDates(start.date)

  f <- function(N){

    start.date + ceiling( cumsum( rexp( N, rate=rate ) ) )

  }

  r <- if(rate < 0.01) r else round(rate,2)

  text <- paste("a poisson process with rate=", 

                 r, " and start date=", start.date,".",sep="")

  new("AccrualGenerator",f=f,

      model="Poisson process",text=text)    

}

##' Create An AccrualGenerator using a power law recruitment

##' 

##' Subjects are accrued according to the c.d.f 

##' \code{G(t)=t^k/B^k} where \code{k} is a parameter, 

##' \code{t} is the time and \code{B} is the recruitment period. 

##'  

##' See the predict from data vignette for further details  

##' @param start.date The start of the subject accrual period

##' @param end.date The date the last subject is accrued so \code{B} = 

##' \code{end.date} - \code{start.date} unless \code{rec.start.date} is used see 

##' below

##' @param k The non-uniformity accrual parameter

##' @param deterministic Logical, if FALSE then the recruitment times 

##' are non-stochastically chosen so that their cumulative distribution function is \code{G(t)}

##' otherwise they are generated by sampling random variables with a cdf \code{G(t)}

##' @param rec.start.date If this argument is used the subjects are still recruited between 

##' start.date and end.date but they follow the cdf \code{G(t)=(t^k-L^k)/(B^k-L^k)} where

##' \code{t} is in \code{[L,B]} and \code{B = end.date - rec.start.date} and \code{L = start.date- rec.start.date}.

##' @return An AccrualGenerator object which generates the required subject accruals

##' @export

Generate.Accrual <- function(start.date,end.date,k,deterministic=FALSE,rec.start.date=NULL){

  start.date <- FixDates(start.date)

  end.date <- FixDates(end.date)

  rec.start.date <- if(is.null(rec.start.date)) start.date else FixDates(rec.start.date)

  # Length of recruitment period

  B <- as.integer(end.date - rec.start.date)

  L <- as.integer(start.date - rec.start.date)

  if(L < 0 || k <= 0|| B <= L){

    stop("Invalid arguments")

  }

  if(deterministic){

    f <- function(N){

      as.Date(((1:N/N)*(B^k-L^k) + L^k)^(1/k)-L,origin=start.date)

    }

    stext <- "non-stochastic"

  }

  else{

    f <- function(N){

      ans <- ((B^k-L^k)*runif(N) + L^k)^(1/k)-L

      ans[order(ans)]

      as.Date(ans,origin=start.date)

    }

    stext <- "stochastic"

  }  

  text <- paste("a ",stext  ," allocation following G(t)=t^k/B^k, where k=", 

  round(k,2), " and B=",B," days is the recruitment period [", rec.start.date, ", ", end.date, "].",sep="")

  if(L!=0){

    text <- paste(text," New subjects will be recruited after ",start.date,".",sep="") 

  }

  new("AccrualGenerator",f=f,

      model="Power law allocation",text=text) 

} 

##' Generic function to estimate non-uniformity accrual 

##' parameter \code{k}

##' 

##' Note The estimate produced assumes recruitment has completed.

##' @docType methods

##' @param object An \code{EventData} object

##' @param ... Additional arguments to be passed to specific method

##' @rdname estimateAccrualParameter-methods  

##' @export

setGeneric("estimateAccrualParameter",function(object,...) standardGeneric("estimateAccrualParameter"))

##' @param start.date The startdate for recruitment, by default the date 

##' the first subject is recruited. If subjects are recruited on start.date, it is

##' assumed they are recruited on the following day to avoid log(0) in the 

##' maximum likelihood estimate.

##' @param end.date The end date for recruitment, by default the date 

##' the last subject is recruited. 

##' @return A maximum likelihood estimate for k, see vignette 

##' for further details    

##' @rdname estimateAccrualParameter-methods 

##' @aliases estimateAccrualParameter,EventData-method

##' @export

setMethod("estimateAccrualParameter", "EventData",

  function(object, start.date=min(object@subject.data$rand.date),

           end.date=max(object@subject.data$rand.date)){

    message("Note this estimate assumes recruitment has completed")          

    start.date <- FixDates(start.date)

    end.date <- FixDates(end.date)

    B <- as.numeric(end.date-start.date)

    if(B <= 0) stop("end.date must be after start.date")

    recs <- (object@subject.data$rand.date-start.date)

    recs <- ifelse(recs==0,0.5,recs)

    if(any(recs< 0 | recs > B)) stop("Subjects recruited outside of recruitment period")

    return(MLestimateK(B,as.numeric(recs)))

  }

)