R/bsd_package.R

Defines functions de.trans de.shift de.birth de.death de.particleT solve.trans solve.shift solve.birth solve.death solve.particleT makeGrid.trans makeGrid.shift.r1 makeGrid.birth.r1 makeGrid.death.r1 makeGrid.particleT.r0 makeGrid.shift.partial makeGrid.birth.partial makeGrid.death.partial makeGrid.particleT.partial getTrans.timeList getParticleT.timeList getDeathMeans.timeList getShiftMeans.timeList getBirthMeans.timeList getTrans.initList getBirthMeans.initList getShiftMeans.initList getDeathMeans.initList getParticleT.initList simulate.true sim.one.true sim.N.true sim.one.ran ransim.N.true sim.eventcount sim.one.eventcount sim.N.eventcount makedata.simple PatientDesignExample MakePatientRates MakePatientData getTrans.MC logFFT.patients FFT.pat.optim ESTEP ESTEP.slow MSTEP EM.run getStandardErrors getTrans.FreqMon logFM FM.data FM.optim getCover getEst FM.run FM.replicate logFFT FFT.optim FFT.run FFT.replicate ESTEP.realdata

Documented in de.birth de.death de.particleT de.shift de.trans EM.run ESTEP ESTEP.slow FFT.optim FFT.pat.optim FFT.replicate FFT.run FM.data FM.optim FM.replicate FM.run getBirthMeans.initList getBirthMeans.timeList getCover getDeathMeans.initList getDeathMeans.timeList getEst getParticleT.initList getParticleT.timeList getShiftMeans.initList getShiftMeans.timeList getStandardErrors getTrans.FreqMon getTrans.initList getTrans.MC getTrans.timeList logFFT logFFT.patients logFM makedata.simple makeGrid.birth.partial makeGrid.birth.r1 makeGrid.death.partial makeGrid.death.r1 makeGrid.particleT.partial makeGrid.particleT.r0 makeGrid.shift.partial makeGrid.shift.r1 makeGrid.trans MakePatientData MakePatientRates MSTEP PatientDesignExample ransim.N.true sim.eventcount sim.N.eventcount sim.N.true sim.one.eventcount sim.one.ran sim.one.true simulate.true solve.birth solve.death solve.particleT solve.shift solve.trans

######################################
##### functions for solving ODEs #####
######################################
#this is the ODE for our probability generating function (pgf) used for calculating transition probabilities

# arguments: 
# t is a number or vector of times to be evaluated
# param is a vector containing b/s/d rates lam,v,mu and argument s2
# state contains the state vectors at their initial conditions
# returns the rate of change in list form
# for further info see deSolve documentation/vignette

#' Transition probability ODE
#' 
#' Evaluates the ODE for generating function corresponding to transition probabilities. 
#' This is in a format to be solved using zvode from package deSolve; see deSolve documentation/vignette
#' for further details
#' 
#' @param t A number or vector of numbers: evaluation times of ODE
#' @param state A named vector containing initial value of the state variable G, complex valued
#' @param param A named vector of numbers containing the other arguments lam, v, mu, and complex number s2
#' @return The rate of change of G the generating function in list form
#' @examples
#' t = .5; state = c(G=exp(2*1i)); param = c(lam = .2, v = .05, mu = .1, s2 = exp(3*1i))
#' de.trans(t,state,param)
de.trans <- function(t, state, param){  
  with(as.list(c(state,param)), {    
    C <- 1/(s2-1+.5*.Machine$double.eps) + lam/(lam-mu)
    f <- 1 + 1/(lam/(mu-lam) + C*exp((mu-lam)*t))
    dG <- G*(lam*f - lam - v - mu) + v*f + mu
    #return rates of change
    list(dG)
  }) #end with as.list
}

#the following are the ODEs to be used for mean 
#births, mean deaths, mean shifts, and particle times (restricted moments)

#' Expected shifts ODE
#' 
#' Evaluates the ODE for generating function corresponding to transition probabilities. 
#' This is in a format to be solved using zvode from package deSolve; see deSolve documentation/vignette
#' for further details
#' 
#' @param t A number or vector of numbers: evaluation times of ODE
#' @param state A named vector containing initial value of the state variable G, complex valued
#' @param param A named vector of numbers containing the other arguments lam, v, mu, r, and complex number s2
#' @return The rate of change of G, the generating function for expected shifts, in list form
#' @examples
#' t = .5; state = c(G=exp(2*1i)); param = c(lam = .2, v = .05, mu = .1, r = 3, s2 = exp(3*1i))
#' de.shift(t,state,param)
de.shift <- function(t, state, param){  
  with(as.list(c(state,param)), {    
    C <- 1/(s2-1+.5*.Machine$double.eps) + lam/(lam-mu)
    f <- 1 + 1/(lam/(mu-lam) + C*exp((mu-lam)*t))
    dG <- G*(lam*f - lam - v - mu) + v*r*f + mu
    #return rates of change
    list(dG)
  }) #end with as.list
}

#' Expected births ODE
#' 
#' Evaluates the ODE for generating function corresponding to transition probabilities. 
#' This is in a format to be solved using zvode from package deSolve; see deSolve documentation/vignette
#' for further details
#' 
#' @param t A number or vector of numbers: evaluation times of ODE
#' @param state A named vector containing initial value of the state variable G, complex valued
#' @param param A named vector of numbers containing the other arguments lam, v, mu, r, and complex number s2
#' @return The rate of change of G, the generating function for expected births, in list form
#' @examples
#' t = .5; state = c(G=exp(2*1i)); param = c(lam = .2, v = .05, mu = .1, r = 3, s2 = exp(3*1i))
#' de.birth(t,state,param)
de.birth <- function(t, state, param){
  with(as.list(c(state,param)), {    
    yb <- (lam + mu + sqrt(lam^2 + 2*lam*mu + mu^2 - 4*lam*mu*r + .5*.Machine$double.eps)) / (2*lam*r)
    C <- 1/(s2-yb) + lam*r/(2*lam*r*yb - lam - mu)
    f <- yb + 1/( -lam*r/(2*lam*r*yb - lam - mu) + C*exp(-(2*yb*lam*r - lam - mu)*t) )
    dG <- G*(lam*r*f - lam - v - mu) + v*f + mu
    #return rates of change
    list(dG)
  }) #end with as.list
}

#' Expected deaths ODE
#' 
#' Evaluates the ODE for generating function corresponding to transition probabilities. 
#' This is in a format to be solved using zvode from package deSolve; see deSolve documentation/vignette
#' for further details
#' 
#' @param t A number or vector of numbers: evaluation times of ODE
#' @param state A named vector containing initial value of the state variable G, complex valued
#' @param param A named vector of numbers containing the other arguments lam, v, mu, r, and complex number s2
#' @return The rate of change of G, the generating function for expected deaths, in list form
#' @examples
#' t = .5; state = c(G=exp(2*1i)); param = c(lam = .2, v = .05, mu = .1, r = 3, s2 = exp(3*1i))
#' de.death(t,state,param)
de.death <- function(t, state, param){
  with(as.list(c(state,param)), {    
    yd <- (lam + mu + sqrt(lam^2 + 2*lam*mu + mu^2 - 4*lam*mu*r + .5*.Machine$double.eps)) / (2*lam)
    C <- 1/(s2-yd) + lam/(2*lam*yd - lam - mu)
    f <- yd + 1/( -lam/(2*lam*yd - lam - mu) + C*exp(-(2*yd*lam - lam - mu)*t) )
    dG <- G*(lam*f - lam - v - mu) + v*f + mu*r
    #return rates of change
    list(dG)
  }) #end with as.list
}

#' Expected particle time ODE
#' 
#' Evaluates the ODE for generating function corresponding to transition probabilities. 
#' This is in a format to be solved using zvode from package deSolve; see deSolve documentation/vignette
#' for further details
#' 
#' @param t A number or vector of numbers: evaluation times of ODE
#' @param state A named vector containing initial value of the state variable G, complex valued
#' @param param A named vector of numbers containing the other arguments lam, v, mu, r, and complex number s2
#' @return The rate of change of G, the generating function for expected particle time, in list form
#' @examples
#' t = .5; state = c(G=exp(2*1i)); param = c(lam = .2, v = .05, mu = .1, r = 3, s2 = exp(3*1i))
#' de.particleT(t,state,param)
de.particleT <- function(t, state, param){
  with(as.list(c(state,param)), {    
    y <- (lam + mu + r +sqrt( (lam + mu + r)^2 - 4*lam*mu + .5*.Machine$double.eps)) / (2*lam)
    C <- 1/(s2-y+.Machine$double.eps) + lam/(2*lam*y - lam - mu - r)
    f <- y + 1/( -lam/(2*lam*y - lam - mu - r) + C*exp(-(2*y*lam - lam - mu - r)*t) )
    dG <- G*(lam*f - lam - v - mu - r) + v*f + mu
    #return rates of change
    list(dG)
  }) #end with as.list
}

####################################
###The following solve these pgfs###
####################################

#solves the pgf for transitiions: param needs to be a vector of (lam,v,mu), time is the interval length
#numsteps is steps taken by deSolve, s1 s2 are arguments

#' Transition ODE solver
#' 
#' Numerically solves the ODE defined in \code{\link{de.trans}} using \code{zvode} from package 
#' \code{deSolve} given evaluation time, initial state values s1 and s2, and birth/shift/death rates
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1 A complex number giving the initial value of the ODE G
#' @param s2 A complex number
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return The function value of the transition probability generating function 
#' @examples 
#' time = 5;  dt = 5; s1 = exp(2*1i); s2 = exp(3*1i); lam = .5; v = .2; mu = .4
#' solve.trans(time,dt,s1,s2,lam,v,mu)
solve.trans <- function(time, dt, s1, s2, lam,v,mu){
  t <- seq(0, time, by=dt)
  state <- c(G = s1)
  param <- c(lam=lam, v=v, mu=mu, s2=s2)
  out <- zvode(y = state, times = t, func = de.trans, parms = param, atol = 1e-10, rtol = 1e-10)
  tail(out, n=1)[2]
}


#' Expected shifts ODE solver
#' 
#' Numerically solves the ODE defined in \code{\link{de.shift}} using \code{zvode} from package 
#' \code{deSolve} given evaluation time, initial state values s1 and s2, and birth/shift/death rates
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1 A complex number giving the initial value of the ODE G
#' @param s2 A complex number
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @param r A real number
#' @return The function value of the shifts generating function 
#' @examples 
#' time = 5;  dt = 1; s1 = exp(2*1i); s2 = exp(3*1i); lam = .5; v = .2; mu = .4; r = 2
#' solve.shift(time,dt,s1,s2, r, lam,v,mu)
solve.shift <- function(time, dt, s1, s2, r, lam, v, mu){
  t <- seq(0, time, by=dt)
  state <- c(G = s1)
  param <- c(lam=lam, v=v, mu=mu, s2=s2, r = r)
  out <- zvode(y = state, times = t, func = de.shift, parms = param, atol = 1e-10, rtol = 1e-10)
  tail(out, n=1)[2]
}


#' Expected births ODE solver
#' 
#' Numerically solves the ODE defined in \code{\link{de.birth}} using \code{zvode} from package 
#' \code{deSolve} given evaluation time, initial state values s1 and s2, and birth/shift/death rates
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1 A complex number giving the initial value of the ODE G
#' @param s2 A complex number
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @param r a real number
#' @return The function value of the births generating function 
#' @examples 
#' time = 5;  dt = 1; s1 = exp(2*1i); s2 = exp(3*1i); lam = .5; v = .2; mu = .4; r = 3
#' solve.birth(time,dt,s1,s2,r,lam,v,mu)
solve.birth <- function(time, dt, s1, s2, r, lam, v, mu){
  t <- seq(0, time, by=dt)
  state <- c(G = s1)
  param <- c(lam=lam, v=v, mu=mu, s2=s2, r = r)
  out <- zvode(y = state, times = t, func = de.birth, parms = param, atol = 1e-10, rtol = 1e-10)
  tail(out, n=1)[2]
}

#' Expected deaths ODE solver
#' 
#' Numerically solves the ODE defined in \code{\link{de.death}} using \code{zvode} from package 
#' \code{deSolve} given evaluation time, initial state values s1 and s2, and birth/shift/death rates
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1 A complex number giving the initial value of the ODE G
#' @param s2 A complex number
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @param r a real number
#' @return The function value of the deaths generating function 
#' @examples 
#' time = 5;  dt = 1; s1 = exp(2*1i); s2 = exp(3*1i); lam = .5; v = .2; mu = .4; r = 3
#' solve.death(time,dt,s1,s2,r,lam,v,mu)
solve.death <- function(time, dt, s1, s2, r, lam, v, mu){
  t <- seq(0, time, by=dt)
  state <- c(G = s1)
  param <- c(lam=lam, v=v, mu=mu, s2=s2, r = r)
  out <- zvode(y = state, times = t, func = de.death, parms = param, atol = 1e-10, rtol = 1e-10)
  tail(out, n=1)[2]
}

#' Expected births ODE solver
#' 
#' Numerically solves the ODE defined in \code{\link{de.particleT}} using \code{zvode} from package 
#' \code{deSolve} given evaluation time, initial state values s1 and s2, and birth/shift/death rates
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1 A complex number giving the initial value of the ODE G
#' @param s2 A complex number
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @param r a real number
#' @return The function value of the particle time generating function 
#' @examples 
#' time = 5;  dt = 1; s1 = exp(2*1i); s2 = exp(3*1i); lam = .5; v = .2; mu = .4; r = 3
#' solve.particleT(time,dt,s1,s2,r,lam,v,mu)
solve.particleT <- function(time, dt, s1, s2, r, lam, v, mu){
  t <- seq(0, time, by=dt)
  state <- c(G = s1)
  param <- c(lam=lam, v=v, mu=mu, s2=s2, r = r)
  out <- zvode(y = state, times = t, func = de.particleT, parms = param, atol = 1e-10, rtol = 1e-10)
  tail(out, n=1)[2]
}


#' Evaluates transition probability ODE over 2D grid of arguments
#' 
#' Applies the function \code{\link{solve.trans}} to a grid of inputs s1, s2 for a fixed time.
#' 
#' @param time A number corresponding to the desired evaluation time of ODE
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.trans(time,dt,s1.seq,s2.seq,lam,v,mu)

makeGrid.trans <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- solve.trans(time, dt, s1.seq[i], s2.seq[j], lam,v,mu) #put the result in the matrix
    }
  }
  return(result)
}

#' Evaluates expected shifts ODE over 2D grid of arguments
#' 
#' Applies the function \code{\link{solve.shift}} to a grid of inputs s1, s2 at one fixed time and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.shift.r1(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.shift.r1 <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- solve.shift(time, dt, s1.seq[i], s2.seq[j], 1.0, lam,v,mu) #put the result in the matrix
    }
  }
  return(result)
}

#' Evaluates expected births ODE over 2D grid of arguments
#' 
#' Applies the function \code{\link{solve.birth}} to a grid of inputs s1, s2 at one fixed time and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.birth.r1(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.birth.r1 <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- solve.birth(time, dt, s1.seq[i], s2.seq[j], 1.0, lam,v,mu) #put the result in the matrix
    }
  }
  return(result)
}

#' Evaluates expected deaths ODE over 2D grid of arguments
#' 
#' Applies the function \code{\link{solve.death}} to a grid of inputs s1, s2 at one fixed time and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.death.r1(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.death.r1 <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- solve.death(time, dt, s1.seq[i], s2.seq[j], 1.0, lam,v,mu) #put the result in the matrix
    }
  }
  return(result)
}

#' Evaluates expected particle time ODE over 2D grid of arguments
#' 
#' Applies \code{\link{solve.particleT}} to a grid of inputs s1, s2 at one fixed time and r=0
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.particleT.r0(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.particleT.r0 <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- solve.particleT(time, dt, s1.seq[i], s2.seq[j], 0.0, lam,v,mu) #put the result in the matrix
    }
  }
  return(result)
}

#' Evaluates partial derivative of expected shifts ODE over 2D grid
#' 
#' Numerically partially differentiates the solution
#' given in \code{\link{solve.shift}} over a grid of input values s1, s2 at one fixed time, and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.shift.partial(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.shift.partial <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- grad(solve.shift, 1.0, time = time, dt = dt, s1 = s1.seq[i], s2 = s2.seq[j], lam = lam, v = v, mu=mu, method = "simple")
    }
  }
  return(result)
}

#' Evaluates partial derivative of expected births ODE over 2D grid
#' 
#' Numerically partially differentiates the solution
#' given in \code{\link{solve.birth}} over a grid of input values s1, s2 at one fixed time, and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.birth.partial(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.birth.partial <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- grad(solve.birth, 1.0, time = time, dt = dt, s1 = s1.seq[i], s2 = s2.seq[j], lam = lam, v = v, mu=mu, method = "simple")
    }
  }
  return(result)
}

#' Evaluates partial derivative of expected deaths ODE over 2D grid
#' 
#' Numerically partially differentiates the solution
#' given in \code{\link{solve.death}} over a grid of input values s1, s2 at one fixed time, and r=1
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.death.partial(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.death.partial <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- grad(solve.death, 1.0, time = time, dt = dt, s1 = s1.seq[i], s2 = s2.seq[j], lam = lam, v = v, mu=mu, method = "simple")
    }
  }
  return(result)
}

#' Evaluates partial derivative of expected particle time ODE over 2D grid
#' 
#' Numerically partially differentiates the solution
#' given in \code{\link{solve.particleT}} over a grid of input values s1, s2 at one fixed time, and r=0
#' 
#' @param time A number corresponding to the desired evaluation time of ODEs
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param lam Birth rate
#' @param v Shift rate
#' @param mu Death rate
#' @return A matrix of dimension length(s1.seq) by length(s2.seq) of the function values
#' 
#' @examples 
#' time = 5;  dt = 5; lam = .5; v = .2; mu = .4
#' gridLength = 32
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' makeGrid.particleT.partial(time,dt,s1.seq,s2.seq,lam,v,mu)
makeGrid.particleT.partial <- function(time, dt, s1.seq, s2.seq, lam, v, mu){
  result <- matrix(nrow = length(s1.seq), ncol = length(s2.seq)) #this will store the grid
  for(i in 1:length(s1.seq)){
    for(j in 1:length(s2.seq)){    
      result[i,j] <- grad(solve.particleT, 0.0, time = time, dt = dt, s1 = s1.seq[i], s2 = s2.seq[j], lam = lam, v = v, mu=mu, method = "simple")
    }
  }
  return(result)
}

############################################################################
### The following functions create transition probability approximations ###
### by transforming solutions from makeGrid functions via FFT            ###
############################################################################

#' Compute transition probabilities over a grid at a list of evaluation times
#' 
#' Transforms the output matrix from \code{\link{makeGrid.trans}} to interpretable 
#' transition probabilities using \code{fft}. Does this at several input times given in tList.
#' Returns a list of matrices, where each list entry corresponds to a time in tList. 
#' The i,j entry of each matrix corresponds to the probability of transitioning from initNum type 1 particles 
#' and 0 type 2 particles to i type 1 particles, j type 2 particles, over the corresponding time length.
#' 
#' @param tList A vector of numbers corresponding to the desired evaluation times
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param initNum An integer giving the number of initial particles
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq), each list entry corresponds to an evaluation time from tList, each matrix is a matrix of transition probabilities
#' 
#' @examples 
#' tList = c(1,2);  dt = 1; lam = .5; v = .2; mu = .4; initNum = 10
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getTrans.timeList(tList, lam, v, mu, initNum, s1.seq, s2.seq, dt)
getTrans.timeList <- function(tList, lam, v, mu, initNum, s1.seq,s2.seq, dt){
  tpm.list <- vector("list", length(tList))      #store the trans prob matrices
  for(i in 1:length(tList)){
    grid.de <- makeGrid.trans(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid.de <- grid.de^initNum
    fourier.de <- fft(grid.de)/length(grid.de)
    tpm.list[[i]] <- Re(fourier.de)
  }
  return(tpm.list)
}

#' Compute expected particle time over a grid at a list of evaluation times
#' 
#' Transforms the output matrices from \code{\link{makeGrid.particleT.r0}} and \code{\link{makeGrid.particleT.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several input times given in tList.
#' Returns a list of matrices, where each list entry corresponds to a time in tList. 
#' The i,j entry of each matrix corresponds to the expected particle time spent with i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, over the corresponding time length.
#' 
#' @param tList A vector of numbers corresponding to the desired evaluation times
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param initNum An integer giving the number of initial particles
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an evaluation time from tList
#' @examples 
#' tList = c(1,2);  dt = 1; lam = .5; v = .2; mu = .4; initNum = 10
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getParticleT.timeList(tList, lam, v, mu, initNum, s1.seq, s2.seq, dt)
getParticleT.timeList <- function(tList, lam, v, mu, initNum, s1.seq,s2.seq, dt){
  tpm.list <- vector("list", length(tList))      #store the trans prob matrices
  for(i in 1:length(tList)){
    grid.partial <- makeGrid.particleT.partial(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid <- makeGrid.particleT.r0(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid.de <- initNum*grid.partial*grid^(initNum-1)
    fourier.de <- fft(grid.de)/length(grid.de)
    tpm.list[[i]] <- Re(fourier.de)
  }
  return(tpm.list)
}

#' Compute expected deaths over a grid at a list of evaluation times
#' 
#' Transforms the output matrices from \code{\link{makeGrid.death.r1}} and \code{\link{makeGrid.death.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several input times given in tList.
#' Returns a list of matrices, where each list entry corresponds to a time in tList. 
#' The i,j entry of each matrix corresponds to the number of expected deaths when process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, at the end of corresponding time interval.
#' 
#' @param tList A vector of numbers corresponding to the desired evaluation times
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param initNum An integer giving the number of initial particles
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an evaluation time from tList
#' 
#' @examples 
#' tList = c(1,2);  dt = 1; lam = .5; v = .2; mu = .4; initNum = 10
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getDeathMeans.timeList(tList, lam, v, mu, initNum, s1.seq, s2.seq, dt)
getDeathMeans.timeList <- function(tList, lam, v, mu, initNum, s1.seq,s2.seq, dt){
  tpm.list <- vector("list", length(tList))      #store the trans prob matrices
  for(i in 1:length(tList)){
    grid.partial <- makeGrid.death.partial(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid <- makeGrid.death.r1(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid.de <- initNum*grid.partial*grid^(initNum-1)
    fourier.de <- fft(grid.de)/length(grid.de)
    tpm.list[[i]] <- Re(fourier.de)
  }
  return(tpm.list)
}

#' Compute expected shifts over a grid at a list of evaluation times
#' 
#' Transforms the output matrices from \code{\link{makeGrid.shift.r1}} and \code{\link{makeGrid.shift.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several input times given in tList.
#' Returns a list of matrices, where each list entry corresponds to a time in tList. 
#' The i,j entry of each matrix corresponds to the number of expected shifts when process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, at the end of corresponding time interval.
#' 
#' @param tList A vector of numbers corresponding to the desired evaluation times
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param initNum An integer giving the number of initial particles
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an evaluation time from tList
#' 
#' @examples 
#' tList = c(1,2);  dt = 1; lam = .5; v = .2; mu = .4; initNum = 10
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getShiftMeans.timeList(tList, lam, v, mu, initNum, s1.seq, s2.seq, dt)
getShiftMeans.timeList <- function(tList, lam, v, mu, initNum, s1.seq,s2.seq, dt){
  tpm.list <- vector("list", length(tList))      #store the trans prob matrices
  for(i in 1:length(tList)){
    grid.partial <- makeGrid.shift.partial(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid <- makeGrid.shift.r1(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid.de <- initNum*grid.partial*grid^(initNum-1)
    fourier.de <- fft(grid.de)/length(grid.de)
    tpm.list[[i]] <- Re(fourier.de)
  }
  return(tpm.list)
}

#' Compute expected birth over a grid at a list of evaluation times
#' 
#' Transforms the output matrices from \code{\link{makeGrid.birth.r1}} and \code{\link{makeGrid.birth.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several input times given in tList.
#' Returns a list of matrices, where each list entry corresponds to a time in tList. 
#' The i,j entry of each matrix corresponds to the number of expected births when process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, at the end of corresponding time interval.
#' 
#' @param tList A vector of numbers corresponding to the desired evaluation times
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param initNum An integer giving the number of initial particles
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an evaluation time from tList
#' 
#' @examples 
#' tList = c(1,2);  dt = 1; lam = .5; v = .2; mu = .4; initNum = 10
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getBirthMeans.timeList(tList, lam, v, mu, initNum, s1.seq, s2.seq, dt)
getBirthMeans.timeList <- function(tList, lam, v, mu, initNum, s1.seq,s2.seq, dt){
  tpm.list <- vector("list", length(tList))      #store the trans prob matrices
  for(i in 1:length(tList)){
    grid.partial <- makeGrid.birth.partial(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid <- makeGrid.birth.r1(time = tList[i], dt, s1.seq, s2.seq, lam,v,mu)
    grid.de <- initNum*grid.partial*grid^(initNum-1)
    fourier.de <- fft(grid.de)/length(grid.de)
    tpm.list[[i]] <- Re(fourier.de)
  }
  return(tpm.list)
}

#' Compute transition probabilites over a grid at a list of initial particle counts
#' 
#' Transforms the output matrices from \code{\link{makeGrid.trans}} 
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several initial particle counts from initList
#' Returns a list of matrices, where each list entry corresponds to a number of initial particles
#' The i,j entry of each matrix corresponds to the probability of transitioning to i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, by the end of corresponding time interval.
#' 
#' @param initList A vector of integers corresponding to the desired initial particle counts
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param u A number giving the observation interval length, equivalently the time to evaluate the ODEs
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an initial number of particles from initList
#' 
#' @examples 
#' initList = c(10,11) 
#' #gives matrices of transition probabilities corresponding to 10 initial particles and 11 initial particles
#' u = 1; dt = 1; lam = .5; v = .2; mu = .4
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getTrans.initList(u, initList, lam, v, mu, s1.seq, s2.seq, dt)

getTrans.initList <- function(u, initList, lam, v, mu, s1.seq, s2.seq, dt){  #initList list of different initial sizes, and u = dt
  tpm.list <- vector("list", length(initList))      #store the trans prob matrices
  grid.de <- makeGrid.trans(u, dt, s1.seq, s2.seq, lam,v,mu) 
  for(i in 1:length(initList)){
    grid.temp <- grid.de^initList[i]
    tpm.list[[i]] <- Re(fft(grid.temp)/length(grid.temp))
  }
  return(tpm.list)
}

#below, we use the fact that d/dr f^k evaluated at r=1 is equal to k* d/dr(f at 1)*f(1)^(k-1)

#' Compute expected births over a grid at a list of initial particle counts
#' 
#' Transforms the output matrices from \code{\link{makeGrid.birth.r1}} and \code{\link{makeGrid.birth.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several initial particle counts from initList
#' Returns a list of matrices, where each list entry corresponds to a number of initial particles
#' The i,j entry of each matrix corresponds to the expected births when the process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, by the end of corresponding time interval.
#' 
#' @param initList A vector of integers corresponding to the desired initial particle counts
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param u A number giving the observation interval length, equivalently the time to evaluate the ODEs
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an initial number of particles from initList
#' 
#' @examples 
#' initList = c(10,11) 
#' #gives matrices of expected sufficient statistics corresponding to 10 initial particles and 11 initial particles
#' u = 1; dt = 1; lam = .5; v = .2; mu = .4
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getBirthMeans.initList(u, initList, lam, v, mu, s1.seq, s2.seq, dt)
getBirthMeans.initList <- function(u, initList, lam, v, mu, s1.seq, s2.seq, dt){  #initList list of different initial sizes, and u = dt
  tpm.list <- vector("list", length(initList))      #store the trans prob matrices
  grid.partial <- makeGrid.birth.partial(u, dt, s1.seq, s2.seq, lam,v,mu)
  grid <- makeGrid.birth.r1(u, dt, s1.seq, s2.seq, lam,v,mu)
  for(i in 1:length(initList)){
    grid.temp <- initList[i]*grid.partial*grid^(initList[i]-1)
    tpm.list[[i]] <- Re(fft(grid.temp)/length(grid.temp))
  }
  return(tpm.list)
}

#' Compute expected shifts over a grid at a list of initial particle counts
#' 
#' Transforms the output matrices from \code{\link{makeGrid.shift.r1}} and \code{\link{makeGrid.shift.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several initial particle counts from initList
#' Returns a list of matrices, where each list entry corresponds to a number of initial particles
#' The i,j entry of each matrix corresponds to the expected shifts when the process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, by the end of corresponding time interval.
#' 
#' @param initList A vector of integers corresponding to the desired initial particle counts
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param u A number giving the observation interval length, equivalently the time to evaluate the ODEs
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an initial number of particles from initList
#' 
#' @examples 
#' initList = c(10,11) 
#' #gives matrices of expected sufficient statistics corresponding to 10 initial particles and 11 initial particles
#' u = 1; dt = 1; lam = .5; v = .2; mu = .4
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getShiftMeans.initList(u, initList, lam, v, mu, s1.seq, s2.seq, dt)
getShiftMeans.initList <- function(u, initList, lam, v, mu, s1.seq, s2.seq, dt){  #initList list of different initial sizes, and u = dt
  tpm.list <- vector("list", length(initList))      #store the trans prob matrices
  grid.partial <- makeGrid.shift.partial(u, dt, s1.seq, s2.seq, lam,v,mu)
  grid <- makeGrid.shift.r1(u, dt, s1.seq, s2.seq, lam,v,mu)
  for(i in 1:length(initList)){
    grid.temp <- initList[i]*grid.partial*grid^(initList[i]-1)
    tpm.list[[i]] <- Re(fft(grid.temp)/length(grid.temp))
  }
  return(tpm.list)
}

#' Compute expected deaths over a grid at a list of initial particle counts
#' 
#' Transforms the output matrices from \code{\link{makeGrid.death.r1}} and \code{\link{makeGrid.death.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several initial particle counts from initList
#' Returns a list of matrices, where each list entry corresponds to a number of initial particles
#' The i,j entry of each matrix corresponds to the expected deaths when the process has i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, by the end of corresponding time interval.
#' 
#' @param initList A vector of integers corresponding to the desired initial particle counts
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param u A number giving the observation interval length, equivalently the time to evaluate the ODEs
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an initial number of particles from initList
#' 
#' @examples 
#' initList = c(10,11) 
#' #gives matrices of expected sufficient statistics corresponding to 10 initial particles and 11 initial particles
#' u = 1; dt = 1; lam = .5; v = .2; mu = .4
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getDeathMeans.initList(u, initList, lam, v, mu, s1.seq, s2.seq, dt)
getDeathMeans.initList <- function(u, initList, lam, v, mu, s1.seq, s2.seq, dt){  #initList list of different initial sizes, and u = dt
  tpm.list <- vector("list", length(initList))      #store the trans prob matrices
  grid.partial <- makeGrid.death.partial(u, dt, s1.seq, s2.seq, lam,v,mu)
  grid <- makeGrid.death.r1(u, dt, s1.seq, s2.seq, lam,v,mu)
  for(i in 1:length(initList)){
    grid.temp <- initList[i]*grid.partial*grid^(initList[i]-1)
    tpm.list[[i]] <- Re(fft(grid.temp)/length(grid.temp))
  }
  return(tpm.list)
}

#' Compute expected particle time over a grid at a list of initial particle counts
#' 
#' Transforms the output matrices from \code{\link{makeGrid.particleT.r0}} and \code{\link{makeGrid.particleT.partial}}
#' to a list of expected sufficient statistics matrices using \code{fft}. Does this at several initial particle counts from initList
#' Returns a list of matrices, where each list entry corresponds to a number of initial particles
#' The i,j entry of each matrix corresponds to the expected particle time spent with i type 1 particles and j type 2 particles,
#' beginning with initNum type 1 particles, by the end of corresponding time interval.
#' 
#' @param initList A vector of integers corresponding to the desired initial particle counts
#' @param dt A number giving the increment length used in solving the ODE
#' @param s1.seq A vector of complex numbers; initial values of the ODE G
#' @param s2.seq A vector of complex numbers as inputs of s2.seq
#' @param u A number giving the observation interval length, equivalently the time to evaluate the ODEs
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @return A list of matrices of dimension length(s1.seq) by length(s2.seq); each list entry corresponds to an initial number of particles from initList
#' 
#' @examples 
#' initList = c(10,11) 
#' #gives matrices of expected sufficient statistics corresponding to 10 initial particles and 11 initial particles
#' u = 1; dt = 1; lam = .5; v = .2; mu = .4
#' gridLength = 16
#' s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
#' getParticleT.initList(u, initList, lam, v, mu, s1.seq, s2.seq, dt)
getParticleT.initList <- function(u, initList, lam, v, mu, s1.seq, s2.seq, dt){  #initList list of different initial sizes, and u = dt
  tpm.list <- vector("list", length(initList))      #store the trans prob matrices
  grid.partial <- makeGrid.particleT.partial(u, dt, s1.seq, s2.seq, lam,v,mu)
  grid <- makeGrid.particleT.r0(u, dt, s1.seq, s2.seq, lam,v,mu)
  for(i in 1:length(initList)){
    grid.temp <- initList[i]*grid.partial*grid^(initList[i]-1)
    tpm.list[[i]] <- -Re(fft(grid.temp)/length(grid.temp))
    #notice the negative sign
  }
  return(tpm.list)
}


#############################
### Simulation Functions  ###
#############################
#the following function simulates from the birth-shift-death process:
#the "indices" correspond to locations; we just let S=100,000
#returns information at end of time interval as ordered pair (new,old)

#' Simulate a birth-shift-death-process
#' 
#' This function simulates one realization of a birth-shift-death CTMC with per-particle birth, death, and shift rates
#' lambda, nu, and mu. 
#' 
#' The process begins with initNum particles, each assigned a random location indexed by a number between 1 and 100,000.
#' Simulation occurs until time t.end. The process returns an ordered pair giving the number of initial locations (indices)
#' still occupied, followed by the number of new indices present, by t.end.
#' 
#' @param t.end A number giving the length of time for the simulation
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @param initNum An integer, initial total number of particles
#' @return A pair of integers giving the number of initial indices still present followed by the number of
#' new indices present in the population. Otherwise, returns '999' as an error code
#' 
#' @examples
#' simulate.true(2,.2,.12,.15,10)
simulate.true <- function(t.end, lam, v, mu, initNum){
  initInd <- sample(1:100000,initNum)            #locations of initial copy generation
  K <- length(initInd)                          #number of particles at current time
  ind <- initInd                                #indices/locations of current collection
  t.cur <- copies <- shifts <- deaths <- 0
  for(i in 1:99999999){                          
    if(K == 0){
      return(c(0,0))  #entire process dies
    }
    
    copy <- shift <- death <- 0   #rates
    copy <-lam*K; shift <- v*K; death <- mu*K; rates <- c(copy, shift, death)
    t.next <- rexp(1, sum(rates)) #time of next event
    t.cur <- t.cur+t.next
    
    if(t.cur > t.end){ #end and return the population
      overlap <- length(intersect(ind,initInd))
      return( c(overlap, length(ind)-overlap) )  #common elements are number of old/type 1 particles; rest are new type 2
    }
    
    decision <- sample( c(1,2,3), 1, prob = rates/sum(rates))  #ratios determine type of event that occurs
    if(decision == 1){ ind <- c(ind, sample(1:100000,1))       #add a new location or index
                       copies <- copies + 1}
    else if(decision == 2){ ind[sample(1:length(ind),1)] <- sample(1:100000,1) #replace an index randomly
                            shifts <- shifts + 1}
    else if(decision == 3){ ind <- ind[-sample(1:length(ind),1)]  #delete an index randomly
                            deaths <- deaths + 1}
    K <- length(ind)  #update the overall rate
    
  }
  return(999) #for catching errors
}


#' Simulate a BSD process with error catch
#'
#' Wrapper for \code{\link{simulate.true}} that catches any possible errors
#' @inheritParams simulate.true
#' @return A pair of integers giving the number of initial indices still present followed by the number of
#' new indices present in the population
#' @examples
#' sim.one.true(2,.2,.12,.15,10)
sim.one.true <- function(t.end, lam, v, mu, initNum){
  res = 999 #this is a check so that nothing weird happens
  while(res[1] == 999){
    res <- simulate.true(t.end, lam, v, mu, initNum)  }
  return(res)
}

#' Simulate N realizations from the BSD process
#'
#' Repeats the function \code{\link{sim.one.true}} N times using \code{replicate}
#' 
#' @inheritParams simulate.true
#' @param N An integer, the number of realizations to simulate
#' @return A 2 by N matrix; each \emph{i}th column is a pair of integers giving the number of initial indices still present followed by the number of
#' new indices present in the population in the \emph{i}th realization
#' @examples
#' sim.N.true(10,2,.2,.12,.15,10)
sim.N.true<- function(N, t.end,lam,v,mu,initNum){
  replicate(N, sim.one.true(t.end, lam, v, mu, initNum))
}

#this simulation is for use with random initial starting sizes across replications;
#same as sim.one.true but keeps track of initial number as well in return statement

#' Simulate a BSD process, returning initial number as well as end state
#'
#' Wrapper for \code{\link{simulate.true}} that catches any possible errors and returns the initial number.
#' 
#' @inheritParams simulate.true
#' @return A vector containing the initial number of particles, followed by a pair of integers giving the number of initial indices still present followed by the number of
#' new indices present in the population
#' @examples
#' sim.one.ran(2,.2,.12,.15,10)
sim.one.ran <- function(t.end, lam, v, mu, initNum){  #also stores initial number
  res = 999 #this is a check so that nothing weird happens
  while(res[1] == 999){
    res <- simulate.true(t.end, lam, v, mu, initNum)  }
  return(c(initNum,res))
}

#' Simulate N realizations from the BSD process, each with random initial number each time
#'
#' Repeats \code{\link{sim.one.ran}} N times using \code{replicate}. Each replication 
#' begins with a random number of initial particles uniformly generated from a specified range.
#' 
#' @inheritParams simulate.true
#' @param N An integer, the number of realizations to simulate
#' @param range A vector containing possible initial populations, typically a sequence of integers
#' @return A 3 by N matrix; each \emph{i}th column corresponds to the \emph{i}th realization. The first row contains
#' initial numbers, the second row contains the number of original indices still present in the population by t.end,
#' and the third row contains the number of new indices present.
#' @examples
#' ransim.N.true(5,2,.2,.12,.15, seq(7,15))
ransim.N.true<- function(N, t.end,lam,v,mu,range){
  replicate(N, sim.one.ran(t.end, lam, v, mu, sample(range,1)))
}

##############################################################################
### The following functions simulate from the same true process, but return the number of shifts, 
### births, deaths as well as the ending state, necessary for verifying restricted moment calculations
### Not necessary for any inference in itself  
##############################################################################

#' Simulate a BSD process and return sufficient statistics
#' 
#' Simulates a birth-shift-death process and records the number of births, deaths, and shifts, as well as total
#' particle time, over the observation interval until time t.end. Used toward simulation functions that check
#' accuracy of expected restricted moment calculations
#' 
#' @inheritParams simulate.true
#' @return A vector of four numbers giving total copies, shifts, deaths, and particle time, respectively. Returns 
#' the integer '999' as an error code if error occurs.
#' @examples
#' sim.eventcount(2,.2,.12,.15,10)
sim.eventcount <- function(t.end, lam, v, mu, initNum){
  initInd <- sample(1:100000,initNum)      	    #locations of initial copy generation, S =10,000 for now
  K <- length(initInd)                          #number of particles at current time
  ind <- initInd                                #indices/locations of current collection
  t.cur <- copies <- shifts <- deaths <- particleT <- 0 # number of copies/shifts/deaths, and total particle time
  for(i in 1:9999999){
    if(K == 0){
      return(c(copies,shifts,deaths, particleT))	#entire process dies
    }
    
    copy <- shift <- death <- 0   #rates
    copy <-lam*K; shift <- v*K; death <- mu*K; rates <- c(copy, shift, death)
    t.next <- rexp(1, sum(rates))
    t.cur <- t.cur+t.next
    
    if(t.cur > t.end){ #end and return the population
      particleT <- particleT + K*(t.end - t.cur + t.next) #add the remaining chunk of time
      return(c(copies,shifts,deaths,particleT))
    }
    
    decision <- sample( c(1,2,3), 1, prob = rates/sum(rates))	#ratios determine type of event that occurs
    if(decision == 1){ ind <- c(ind, sample(1:100000,1))       #add a new location or index
                       copies <- copies + 1}
    else if(decision == 2){ ind[sample(1:length(ind),1)] <- sample(1:100000,1) #replace an index randomly
                            shifts <- shifts + 1}
    else if(decision == 3){ ind <- ind[-sample(1:length(ind),1)]  #delete an index randomly
                            deaths <- deaths + 1}
    particleT <- particleT + K*t.next
    K <- length(ind)  #update the overall rate
  }
  return(999)
}

#' Simulate a BSD process and return sufficient statistics with error catch
#' 
#' Simulates a birth-shift-death process and records the number of births, deaths, and shifts, as well as total
#' particle time, over the observation interval until time t.end. Used toward simulation functions that check
#' accuracy of expected restricted moment calculations
#' 
#' @inheritParams simulate.true
#' @return A vector of four numbers giving total copies, shifts, deaths, and particle time, respectively.
#' @examples 
#' sim.one.eventcount(2,.2,.12,.15,10)
sim.one.eventcount <- function(t.end, lam, v, mu, initNum){
  res = 999 #this is a check so that nothing weird happens
  while(res[1] == 999){
    res <- sim.eventcount(t.end, lam, v, mu, initNum)  }
  return(res)
}

#' Simulate N BSD processes and return sufficient statistics
#' 
#' Replicates \code{\link{sim.one.eventcount}} N times. Each replication 
#' begins with initNum particles.
#' 
#' @inheritParams simulate.true
#' @param N The number of replications, an integer
#' @return A 4 by N matrix, where rows correspond to total copies, shifts, deaths, and particle time per realization, respectively.
#' Each column corresponds to one realization of the process.
#' @examples
#' sim.N.eventcount(3,2,.2,.12,.15,10)
sim.N.eventcount<- function(N, t.end,lam,v,mu,initNum){
  replicate(N, sim.one.eventcount(t.end, lam, v, mu, initNum))
}


#generates a list of simulated data matrices, each entry of the list corresponds to a time in tList
#each column of each matrix stores one "interval", for N total intervals
#row one is initial number of [old] particles, row 2 is number of type 1 [old] at end of interval, 
#row 3 is type 2 [new] particles

#' Generate synthetic dataset for BSD process with simple rates
#' 
#' \code{makedata.simple} creates a list of synthetic observed datasets, each with N observation intervals
#' from the BSD process. Observation interval lengths are entries in tList. This is used in simulation studies checking
#' inferential procedures. Simulates observations using \code{\link{ransim.N.true}}.
#' 
#' @param N An integer specifying the number of observation intervals/realization from the simple BSD process.
#' @param tList A list of observation interval lengths. The number of datasets returned is equal to the length of tList
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @param initList A vector containing possible initial population sizes
#' @return A list of matrices. Each entry in the list is in the format returned by \code{\link{ransim.N.true}}, and 
#' corresponds to equal observation times specified in the corresponding entry in tList
#' @examples
#' makedata.simple(10,c(.2,.4,.6),.2,.12,.15,seq(8,12))
makedata.simple <- function(N, tList, lam, v, mu, initList){
  simDataList <- vector("list", length(tList))
  for(i in 1:length(tList)){
    simDataList[[i]] <- ransim.N.true(N, tList[i],lam,v,mu,initList)
  }
  return(simDataList)
}

#' Create example design matrix in the format required for \code{\link{MakePatientRates}}
#' 
#' This function creates one possible design matrix as an input for the function \code{\link{MakePatientRates}}
#' This particular design matrix features one covariate uniformly sampled between 6 and 10, as well as a row 
#' of 1's corresponding to an intercept term. The design matrix is thus 2 by number of patients
#' 
#' @param num.patients An integer giving the number of patients in the design matrix
#' @return A 2 by num.patients matrix: the first row is constant 1's corresponding to intercepts, and 
#' second row is a patient-specific covariate sampled uniformly between 6 and 10
#' 
PatientDesignExample <- function(num.patients){
  x1 <- runif(num.patients,6,10)
  patients.design <- rbind(1, x1) 
  return(patients.design)
}

#' Makes matrix of patient-specific birth, shift, and death rates
#' 
#' This function returns a matrix containing the birth, shift, and death rates of each patient. The rates are calculated
#' based on the log-linear relationship between regression coefficients 
#' \eqn{\mathbf{\beta} = (\mathbf{\beta^\lambda}, \mathbf{\beta}^\nu, \mathbf{\beta}^\mu)} 
#' and covariates \eqn{\mathbf{z_p}} in the \emph{p}th column of the design matrix given by
#' \eqn{\log(\lambda_p) = \mathbf\beta^\lambda \cdot \mathbf{z_p},  \log(\nu_p) = \mathbf\beta^\nu \cdot \mathbf{z_p}, \log(\mu_p) = \mathbf\beta^\mu \cdot \mathbf{z_p} }
#' 
#' @param patients.design A \eqn{n \times m} matrix, where n is number of covariates (including intercept) in the model and m is number of patients
#' @param betaVec The vector \eqn{\mathbf\beta} of regression coefficients. Note this function assumes 
#' that \eqn{ \mathbf\beta^\lambda, \mathbf\beta^\nu, \mathbf\beta^\mu} are equal in length. That is,
#' each rate depends on the same number of covariates
#' 
#' @return A \eqn{3 \times m} matrix where m is the number of patients, and rows correspond to birth, shift, and death rates respectively. 
#' 
#' @examples
#' num.patients = 10
#' patients.design <- PatientDesignExample(num.patients)
#' beta.lam <- c(log(8), log(.6)); beta.v <- c( log(.5), log(.7)); beta.mu <- c(log(.8), log(.8))
#' betas <- c(beta.lam,beta.v,beta.mu)
#' MakePatientRates(patients.design, betas)
MakePatientRates <- function(patients.design, betaVec){ 
  beta.lam <- betaVec[1:(length(betaVec)/3)]; beta.v = betaVec[(length(betaVec)/3+1):(2*length(betaVec)/3)]
  beta.mu = betaVec[(2*length(betaVec)/3 + 1) : length(betaVec)]
  return(exp( rbind( colSums(beta.lam*patients.design), colSums(beta.v*patients.design), colSums(beta.mu*patients.design)) ))
}

#' Generates synthetic patient data for inference in discretely observed BSD process with covariates
#' 
#' \code{MakePatientData} is the main function for generating a synthetic dataset with covariates. It simulates observations from
#' a birth-shift-death process for a number of synthetic "patients", using \code{\link{ransim.N.true}}. This provides data
#' for simulation studies toward assessing inference in a panel data setting, with rates depending on multiple covariates.
#' 
#' Each observation interval has a common fixed length given by the argument t. For each patient, between 2 and 6 observation intervals
#' are generated, and each begins with between 2 and 12 initial particles. These numbers are initialized uniformly at random, and
#' passed in as arguments to \code{\link{ransim.N.true}}, with the true rates set to the corresponding column of patients.rates.
#' 
#' @param t A number giving the observation interval length
#' @param num.patients An integer, the number of synthetic patients
#' @param patients.rates A matrix in the format returned by \code{\link{MakePatientRates}}
#' @return A \eqn{5 \times m} matrix where m is the number of total observation intervals generated. The first row
#' corresponds to the patient ID, the second row gives the initial number of particles for that observation interval,
#' the third through fifth column are in the format returned by \code{\link{ransim.N.true}}
#' 
#' @examples
#' num.patients = 10; t = .4
#' patients.design <- PatientDesignExample(num.patients)
#' beta.lam <- c(log(8), log(.6)); beta.v <- c( log(.5), log(.7)); beta.mu <- c(log(.8), log(.8))
#' betas <- c(beta.lam,beta.v,beta.mu)
#' patients.rates <- MakePatientRates(patients.design, betas)
#' MakePatientData(t, num.patients, patients.rates)
MakePatientData <- function(t, num.patients, patients.rates){ 
  res <- c()
  #first create design matrix  
  for(i in 1:num.patients){
    obs <- sample(seq(2,6),1) #how many observations we have for this patient, uniform between 2 and 6
    init.patient <- sample(seq(2,12),1)
    #let's say the possible values per patient are only + or - 1 from the init.patient, which is
    #kind of the "average" starting value for a given patient
    #this is more realistic and easier computationally
    res <- cbind(res, rbind(i, init.patient, ransim.N.true(obs, t, patients.rates[1,i], patients.rates[2,i], patients.rates[3,i], seq(init.patient-1, init.patient+1)) ) )
  }
  return(res)
}

#' Calculate empirical transition probabilities
#' 
#' Calculates transition probabilities empirically using Monte Carlo simulation
#' 
#' Note: function can be modified to initialize matrices to be larger sized in the case where rates are large
#' 
#' @param N An integer, the number of MC simulations per element of tList
#' @param tList A list of observation interval lengths to simulate
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @param initNum Integer giving the number of initial particles
#' @return A list of matrices, where each entry of the list corresponds to an element of tList. The i,j entry of
#' each matrix in the list give the probability of the process ending with i type 1 particles and j type 2 particles, beginning
#' with initNum type 1 particles, by the end of the corresponding observation length.
#' 
#' @examples
#' N = 20; tList = c(.5,1); lam = .2; v = .1; mu = .15; initNum = 10
#' getTrans.MC(N,tList,lam,v,mu,initNum)
getTrans.MC <- function(N, tList, lam,v,mu, initNum){
  tpm.list <- vector("list", length(tList))     #this will store 2 by 2 transitions
  for(t in 1:length(tList)){
    result <- sim.N.true(N, tList[t], lam, v, mu, initNum)
    trans.count <- matrix(0, (initNum+1),(initNum+20))  #make big enough to account for everything that may happen: count end states
    for(i in 1:N){
      id <- result[,i]+1 #indices in the resulting transition count: ie if you end at (1,1), you add a count to the (2,2) entry of the count matrix, etc
      trans.count[id[1], id[2]] = trans.count[id[1], id[2]] + 1
    }
    tpm.list[[t]] <- trans.count/sum(trans.count)
  }
  return(tpm.list)
}


#########################################################################################
### The following functions are relevant to EM algorithm and other simulation studies ###
#########################################################################################

##NOTE may include the real data version of functions where times can vary, etc

#' Log likelihood for BSD process with covariates
#' 
#' \code{logFFT.patients} evaluates the log likelihood of a dataset with observations corresponding to "patients" in the setting
#' where rates of the process depend on patient-specific covariates. The transition probabilities given the states of the process at 
#' endpoints of each observation interval are computed using the FFT/generating function method, relying on \code{\link{getTrans.initList}}. 
#' The log likelihood is then the sum of these log transition probabilities.
#' 
#' Note: this function is used so that MLE estimation of the coefficient vector can be accomplished, i.e. using \code{optim}, and 
#' is also used in numerically computing standard errors at the MLE. 
#' 
#' Vectors s1.seq and s2.seq should be of length greater than the total number of particles of either type at any observation interval
#' 
#' @param betas A vector of numbers \eqn{\mathbf\beta = (\mathbf\beta^\lambda, \mathbf\beta^\nu, \mathbf\beta^\mu)}
#' @param t.pat A number, the observation interval length
#' @param num.patients An integer, number of unique patients
#' @param PATIENTDATA A matrix in the form returned by \code{\link{MakePatientData}} containing the set of observation intervals
#' @param patients.design A design matrix in the same form as returned by \code{\link{PatientDesignExample}}
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @return The negative log likelihood of the observations in PATIENTDATA, given rates determined by coefficient vector betas and 
#' covariate values in patients.design
logFFT.patients <- function(betas, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq){
  b.lam <- betas[1:(length(betas)/3)]; b.v = betas[(length(betas)/3+1):(2*length(betas)/3)]; b.mu = betas[(2*length(betas)/3 + 1) : length(betas)]
  logresult <- 0
  for(i in 1:num.patients){
    lam <- exp( sum(b.lam*patients.design[,i])); v <- exp( sum(b.v*patients.design[,i])); mu <- exp(sum(b.mu*patients.design[,i]))
    #this takes the subset of the fake data corresponding to patient i
    dat.i <- PATIENTDATA[,which(PATIENTDATA[1,]==i)]
    tpm.list <- getTrans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq,t.pat) #stores trans matrices for all possible init size for this patient
    for(j in 1:dim(dat.i)[2]){
      initial <- dat.i[3,j] - dat.i[2,1] + 2  #subtract the init.patient to get the proper index in the tpm list
      n.old <- dat.i[4,j]
      n.new <- dat.i[5,j]
      logresult <- logresult + log(tpm.list[[initial]][n.old+1, n.new+1]) #initial is a list index, and then pick the entries in the matrix chosen
    }
  }
  return(-logresult)
}

#' Optimizes the function \code{\link{logFFT.patients}} using \code{optim} package
#' 
#' Function uses Nelder-Mead optimization as implemented in \code{optim} to maximize the log likelihood function
#' \code{\link{logFFT.patients}}. 
#' 
#' Example code is not included here because of runtime; see vignette for tutorial on using this function.
#' 
#' @param betaInit A vector, the initial guess for the algorithm
#' @param t.pat A number, the observation interval length
#' @param num.patients An integer, number of unique patients
#' @param PATIENTDATA A matrix in the form returned by \code{\link{MakePatientData}} containing the set of observation intervals
#' @param patients.design A design matrix in the same form as returned by \code{\link{PatientDesignExample}}
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @param tol A number for setting the relative tolerance for the algorithm (the reltol argument in \code{optim})
#' @param max An integer, the max number of iterations before termination
#' @return An \code{optim} type object
#optimizes the function logFFT.patients beginning with betaInit as initial guess
#uses Nelder-Mead optimization in optim package
FFT.pat.optim <- function(betaInit, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq, tol, max=2000, hess=TRUE){
  optim(betaInit, logFFT.patients, t.pat= t.pat, num.patients = num.patients, PATIENTDATA = PATIENTDATA, patients.design = patients.design, 
        s1.seq = s1.seq, s2.seq = s2.seq, hessian=hess, control = list(maxit = max, reltol = tol))
}


#' Perform one accelerated E-step of the EM algorithm
#' 
#' \code{ESTEP} performs one E-step of the EM algorithm, computing expected sufficient statistics given current settings
#' of the parameters. This function is the "accelerated" version, meaning that intervals with no observed changes are computed
#' more efficiently using closed form expressions, bypassing generating function and FFT calculations on these intervals.
#' 
#' @param betaVec A vector, the setting of beta coefficients
#' @param t.pat A number, the observation interval length
#' @param num.patients An integer, number of unique patients
#' @param PATIENTDATA A matrix in the form returned by \code{\link{MakePatientData}} containing the set of observation intervals
#' @param patients.design A design matrix in the same form as returned by \code{\link{PatientDesignExample}}
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @return A list containing a matrix of the expected sufficient statistics as well as the observed log likelihood value
ESTEP <- function(betaVec, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq){
  bLam <- betaVec[1:(length(betaVec)/3)]; bNu = betaVec[(length(betaVec)/3+1):(2*length(betaVec)/3)]; bMu = betaVec[(2*length(betaVec)/3 + 1) : length(betaVec)]
  U <- D <- V <- P <- rep(0, num.patients) #initialize vectors of expected sufficient statistics
  observedLikelihood <- 0
  for(i in 1:num.patients){
    #get lam, v, mu for current patients
    lam <- exp( sum(bLam*patients.design[,i]))
    v <- exp( sum(bNu*patients.design[,i])) 
    mu <- exp(sum(bMu*patients.design[,i]))
    #this takes the subset of the fake data corresponding to patient i
    dat.i <- PATIENTDATA[,which(PATIENTDATA[1,]==i)]
    
    #generate the transition probabilities and restricted moments for all possible initial sizes of patients
    #all information for the E step is calculated in these functions: this is the expensive part/meat of update
    tpm.list <- getTrans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat) #stores trans matrices for all possible init size for this patient
    restrictedBirth.list <- getBirthMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedShift.list <- getShiftMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedDeath.list <- getDeathMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedParticle.list <- getParticleT.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    
    #loop over number of observations associated with patient i
    for(j in 1:dim(dat.i)[2]){
      #these three items are indices to look up terms in our transition lists
      initial <- dat.i[3,j] - dat.i[2,1] + 2  #subtract the init.patient to get the proper index in the tpm list
      n.old <- dat.i[4,j]
      n.new <- dat.i[5,j]
      #initial is a list index, and then pick the entries in the matrix chosen
      #we divide all the restricted moments by the corresponding transition probability
      #we add to the i'th entry for each interval associated with patient i in this loop
      
      #if no changes occur, use quicker 'FM' methods
      if( dat.i[3,j] == dat.i[4,j] & n.new == 0){
        #print('no event')
        #zero expected counts to add, but particle time needs to be accounted for
        P[i] <- P[i] + dat.i[3,j]*t.pat   #particle time is just initial particles multiplied by interval length
        observedLikelihood <- observedLikelihood - dat.i[3,j]*(lam+v+mu)*t.pat
      } else {
        #NOTE we should check for any divide by zero etc that may happen here
        observedLikelihood <- observedLikelihood + log(tpm.list[[initial]][n.old+1, n.new+1])
        U[i] <- U[i] + restrictedBirth.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
        D[i] <- D[i] + restrictedDeath.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
        V[i] <- V[i] + restrictedShift.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
        P[i] <- P[i] + restrictedParticle.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1]
        #note the P has a minus
      }
    }
  }
  result <- rbind(U,D,V,P)
  result[result < 0] = .000001 #numerical stability for negatives
  resultList <- list( 'matrix' = result, 'loglike' = observedLikelihood)
  return(resultList)  #return the vectors of expected sufficient statistics
}


#' Perform one E-step of the EM algorithm
#' 
#' \code{ESTEP.slow} performs one E-step of the EM algorithm, computing expected sufficient statistics given current settings
#' of the parameters. This function is the un-accelerated version, simply using the generating function approach to compute
#' necessary quantities for all observation intervals.
#' 
#' @inheritParams ESTEP
#' @return A list containing a matrix of the expected sufficient statistics as well as the observed log likelihood value
ESTEP.slow <- function(betaVec, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq){
  bLam <- betaVec[1:(length(betaVec)/3)]; bNu = betaVec[(length(betaVec)/3+1):(2*length(betaVec)/3)]; bMu = betaVec[(2*length(betaVec)/3 + 1) : length(betaVec)]
  U <- D <- V <- P <- rep(0, num.patients) #initialize vectors of expected sufficient statistics
  observedLikelihood <- 0
  for(i in 1:num.patients){
    #get lam, v, mu for current patients
    lam <- exp( sum(bLam*patients.design[,i]))
    v <- exp( sum(bNu*patients.design[,i])) 
    mu <- exp(sum(bMu*patients.design[,i]))
    #this takes the subset of the fake data corresponding to patient i
    dat.i <- PATIENTDATA[,which(PATIENTDATA[1,]==i)]
    #generate the transition probabilities and restricted moments for all possible initial sizes of patients
    #all information for the E step is calculated in these functions: this is the expensive part/meat of update
    tpm.list <- getTrans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat) #stores trans matrices for all possible init size for this patient
    restrictedBirth.list <- getBirthMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedShift.list <- getShiftMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedDeath.list <- getDeathMeans.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    restrictedParticle.list <- getParticleT.initList(t.pat, seq(dat.i[2,1]-1, dat.i[2,1]+1), lam, v, mu, s1.seq, s2.seq, t.pat)
    
    #loop over number of observations associated with patient i
    for(j in 1:dim(dat.i)[2]){
      #these three items are indices to look up terms in our transition lists
      initial <- dat.i[3,j] - dat.i[2,1] + 2  #subtract the init.patient to get the proper index in the tpm list
      n.old <- dat.i[4,j]
      n.new <- dat.i[5,j]
      #initial is a list index, and then pick the entries in the matrix chosen
      #we divide all the restricted moments by the corresponding transition probability
      #we add to the i'th entry for each interval associated with patient i in this loop
      
      #NOTE we should check for any divide by zero etc that may happen here
      observedLikelihood <- observedLikelihood + log(tpm.list[[initial]][n.old+1, n.new+1])
      U[i] <- U[i] + restrictedBirth.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
      D[i] <- D[i] + restrictedDeath.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
      V[i] <- V[i] + restrictedShift.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1] 
      P[i] <- P[i] + restrictedParticle.list[[initial]][n.old+1, n.new+1]/tpm.list[[initial]][n.old+1, n.new+1]
      #note the P has a minus
    }
  }
  result <- rbind(U,D,V,P)
  result[result < 0] = .000001 #numerical stability for negatives
  resultList <- list( 'matrix' = result, 'loglike' = observedLikelihood)
  return(resultList)  #return the vectors of expected sufficient statistics
}

#performs one newton raphson step: takes in matrix returned by E step and current estimate of parameters as arguments
#returns the new betas after a newton raphson update

#' Execute a Newton-Raphson step within M-step of the EM algorithm
#' 
#' \code{MSTEP} executes one iteration of a Newton-Raphson algorithm as part of the maximization (M-step) of the EM algorithm. 
#' Given the matrix of expected sufficient statistics returned by \code{\link{ESTEP}}, this function uses closed form gradient and
#' hessian expressions to efficiently optimize the current settings of the coefficients beta. This is called up to 10 times per M-step 
#' within \code{\link{EM.run}}
#' 
#' @param matrix A matrix in the format returned by \code{\link{ESTEP}} or \code{\link{ESTEP.slow}}
#' @param betaVec A vector of regression coefficients
#' @param num.patients An integer, the number of unique patients
#' @param patients.design A design matrix in the format generated by \code{\link{PatientDesignExample}}
#' @return An updated coefficient vector after one Newton-Raphson step
MSTEP <- function(matrix, betaVec, num.patients, patients.design){
  U <- matrix[1,]; D <- matrix[2,]; V <- matrix[3,]; P <- matrix[4,]
  bLam <- betaVec[1:(length(betaVec)/3)]; bNu = betaVec[(length(betaVec)/3+1):(2*length(betaVec)/3)]; bMu = betaVec[(2*length(betaVec)/3 + 1) : length(betaVec)]
  lamList <- vList <- muList <- rep(0, num.patients) #initialize vectors of expected sufficient statistics
  for(i in 1:num.patients){
    #get lam, v, mu for current patients
    lamList[i] <- exp( sum(bLam*patients.design[,i]))
    vList[i] <- exp( sum(bNu*patients.design[,i])) 
    muList[i] <- exp(sum(bMu*patients.design[,i]))
  }
  
  gradient <- c( patients.design %*% ( -diag(P) %*% lamList + U), patients.design %*% ( -diag(P) %*% vList + V), patients.design %*% ( -diag(P) %*% muList + D) )
  lamBlock <- -patients.design %*% diag(P) %*% diag(lamList) %*% t(patients.design)
  nuBlock <- -patients.design %*% diag(P) %*% diag(vList) %*% t(patients.design)
  muBlock <- -patients.design %*% diag(P) %*% diag(muList) %*% t(patients.design)
  hess <- bdiag(lamBlock,nuBlock,muBlock)
  hessInverse <- solve(hess)
  betaNew <- betaVec - as.vector( hessInverse %*% gradient )
  return(betaNew)
}

###runs the EM algorithm, with initial guess betaInit
#returns log likelihood and beta estimate at convergence, and number of iterations

#' Inference via EM algorithm
#' 
#' This function runs the EM algorithm from an initial guess. Infers the coefficient vector in setting where rates
#' depend on patient-specific covariates. The EM algorithm alternates between calling \code{\link{ESTEP}} and \code{\link{MSTEP}}
#' until the change in observed log-likelihood changes less than a specified relative tolerance between iterations
#' 
#' Examples are not included here due to runtime, but see vignette for usage.
#' 
#' @param betaInit A vector, the initial guess for coefficients beta
#' @param t.pat A number, the observation interval length
#' @param num.patients An integer, number of unique patients
#' @param PATIENTDATA A matrix in the form returned by \code{\link{MakePatientData}} containing the set of observation intervals
#' @param patients.design A design matrix in the same form as returned by \code{\link{PatientDesignExample}}
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @param relTol A number, the relative convergence criterion
#' 
#' @return A list containing the log-likelihood value at convergence, the final beta estimate, and the number of iterations
EM.run <- function(betaInit, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq, relTol, verbose=FALSE){
  llist <- c(999,99) #records log likelihoods, ignore these first two entries
  betacur <- bMat <- betaInit #betacur will keep track of initial estimate, bMat records all estimates
  mat <- mat.or.vec(4,num.patients) 
  iter = 2
  
  while(abs(llist[iter]-llist[iter-1]) > relTol*abs(llist[iter])){
    ###E STEP
    mat <- ESTEP(betacur, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq)
    if(verbose){print(betacur)}
    ###M step
    #initialize betaold, a vector to store previous beta setting after each Newton-Raphson step
    betaold <- rep(99,length(betacur)); numMsteps = 0
    #convergence criterion based on relTol for number of N-R steps, also max set at 10      
    while( sum(abs(betaold - betacur)) > relTol*sum(abs(betacur)) & numMsteps < 10){
      betaold <- betacur
      betacur <- MSTEP(mat$matrix,betacur, num.patients, patients.design)
      numMsteps = numMsteps + 1
    }
    #print(paste("M steps:", numMsteps))
    llist <- c(llist, mat$loglike )
    bMat <- rbind(bMat,as.vector(betacur))
    iter <- iter +1
  }
  
  #list of results: converged log like, final beta estimate, and number of iterations
  resList <- list('loglike' = llist[length(llist)], 'betaEstimate' =  bMat[dim(bMat)[1],] , 'iters' = iter-2)
  return(resList)
}

#calculates fisher info at the value of betas, and returns a vector of standard deviations

#' Calculate numerical standard errors of estimated coefficients
#' 
#' \code{getStandardErrors} uses \code{optimHess} to numerically compute the hessian at the value of the estimated MLE, 
#' which can be obtained using \code{\link{EM.run}}. The function inverts the hessian and takes the square root of diagonal entries
#' to yield standard errors for each coefficient estimate.
#' 
#' @inheritParams logFFT.patients
#' @return A vector the length of the coefficient vector, giving standard errors for each estimated coefficient
getStandardErrors <- function(betas, t.pat, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq){
  hessian <- optimHess(betas, logFFT.patients, t.pat = t.pat, num.patients = num.patients, PATIENTDATA = PATIENTDATA,
                       patients.design = patients.design, s1.seq = s1.seq, s2.seq = s2.seq)
  fisherInfo <- solve(hessian)
  return(sqrt(diag(fisherInfo)))
}


######################################
### Frequent Monitoring functions ####
######################################

#This calculates the four transition probabilities available under frequent monitoring, for comparison

#' Calculates the four transition probabilities defined under Frequent Monitoring
#' 
#' \code{getTrans.FreqMon} uses simple closed form expressions to compute transition probabilities available under the
#' Frequent Monitoring assumption, which allows at most one event to occur per observation interval. Computes the four
#' possible transition possibilities for intervals with lengths corresponding to entries in tList
#' 
#' @param tList A vector of observation interval lengths
#' @param lam Per-particle birth rate
#' @param v Per-particle shift rate
#' @param mu Per-particle death rate
#' @param initNum Integer, the number of initial particles
#' @return A list of \eqn{2 \times 2} matrices, each containing the probabilitiy of one birth, one shift, one death, or no event occurring
#' over a corresponding observation interval length from tList
#' 
#' @examples
#' getTrans.FreqMon(c(.5,1,10),.3,.1,.2,10)
getTrans.FreqMon <- function(tList, lam, v, mu, initNum){
  tpm.list <- vector("list", length(tList))     #this will store 2 by 2 transitions
  theta <- v + mu + lam
  for(u in 1:length(tList)){
    a <-  (1 - exp(-initNum*theta*tList[u]))*mu/theta  #top left (N,0) -> (N-1, 0)
    b <-  exp(-initNum*theta*tList[u])        #bottom left (N,0) -> (N,0)
    c <-  (1 - exp(-initNum*theta*tList[u]))*v/theta   #top right (N,0) -> (N-1,1)
    d <-  (1 - exp(-initNum*theta*tList[u]))*lam/theta   #bototm right (N,0) -> (N,1)
    tpm.list[[u]] <- matrix(c(a,b,c,d), 2, 2)      #append the matrix
  }
  return(tpm.list)
}

#the frequent monitoring log likelihood, u is length of interval, param is (lam, v, mu)
#this can be directly maximized via optim()
#takes in FM object

#' Calculate approximate log-likelihood based on Frequent Monitoring computations
#' 
#' \code{logFM} calculates the log-likelihood where transition probabilities are computed under the frequent monitoring 
#' assumption. Used for simulation studies comparing results using frequent monitoring with our methods.
#' 
#' @param param A vector of three numbers containing the birth, shift, and death rate respectively
#' @param u A number, the observation interval length
#' @param FM An object returned by \code{\link{FM.data}} containing the number of each event under FM
#' @return The negative frequent monitoring log-likelihood
logFM <- function(param, u, FM){
  FM.births <- FM$births; FM.deaths <- FM$deaths; FM.shifts <- FM$shifts; FM.nothing <- FM$nothing
  Bi <- length(FM.births)
  De <- length(FM.deaths)
  Sh <- length(FM.shifts)
  No <- length(FM.nothing)
  lam <- param[1]; v = param[2]; mu = param[3]
  theta <- lam+v+mu
  birthsum <- Bi*log(lam/theta) + sum( log(1 - exp(-FM.births*theta*u)) )
  deathsum <- De*log(mu/theta) + sum( log(1 - exp(-FM.deaths*theta*u)) )
  shiftsum <- Sh*log(v/theta) + sum( log(1 - exp(-FM.shifts*theta*u)) )
  nosum <- sum( log(exp(-FM.nothing*theta*u)))
  return( -(birthsum+deathsum+shiftsum+nosum ) )
}


##This function now returns the quantities necessary for logFM computation
#returns the number of births, deaths. shifts, and no event under FM in a vector called FM

#' Finds all intervals compatible under frequent monitoring assumption in a synthetic dataset
#' 
#' \code{FM.data} takes a simulated dataset in the format generated by \code{\link{makedata.simple}} and finds
#' the subset of observation intervals where at most one event occurred. This is necessary for computation of the
#' log likelihood in \code{\link{logFM}}
#' 
#' @param simDataList A list of synthetic observed datasets, returned by \code{\link{makedata.simple}}
#' @param u The index of the desired entry of simDataList
#' @return A list containing information about each type of FM event
FM.data <- function(simDataList, u){
  data <- simDataList[[u]]
  #takes only the intervals consistent under FM assumption:
  FM <- data[, intersect( which(data[1,] - data[2,] < 2), which(data[3,] < 2) )]  #only keep those that differ by at most 1 old particle
  #split up into births, deaths, shifts, or nothing happening, first row (initial number at that interval)
  FM.births <- FM[, intersect(which(FM[3,]== 1), which(FM[1,] == FM[2,])), drop = FALSE ][1,]
  FM.deaths <- FM[, intersect(which(FM[2,] < FM[1,]), which(FM[3,] == 0)), drop = FALSE ][1,]
  FM.shifts <- FM[, intersect(which(FM[2,] < FM[1,]), which(FM[3,] == 1)), drop = FALSE ][1,]
  FM.nothing <- FM[, intersect(which(FM[3,] == 0), which(FM[1,] == FM[2,])), drop = FALSE ][1,]
  resultList <- list('births' = FM.births, 'deaths' = FM.deaths, 'shifts' = FM.shifts, 'nothing' = FM.nothing)
  return(resultList)
}

#' Optimizes the frequent monitoring log-likelihood
#' 
#' This function converts a dataset generated in the format returned by \code{\link{makedata.simple}} 
#' and finds the intervals compatible under FM using \code{\link{FM.data}}, then infers the MLE birth, 
#' shift, and death rates.
#' 
#' @param simDataList A list of synthetic observed datasets, returned by \code{\link{makedata.simple}}
#' @param u The index of the desired entry of simDataList
#' @param initGuess A vector, initial guess for beta
#' @return An optim object corresponding to maximized frequent monitoring log-likelihood
FM.optim <- function(simDataList, u, initGuess){
  optim(initGuess, logFM, u=tList[u], FM = FM.data(simDataList, u), hessian = TRUE, control = list(trace = 6))
}

#takes in optim type object, returns vector of length 3 of coverage counts for each parameter
#for evaluating coverage post-hoc
#' Calculates coverage indicators for estimated parameters
#' 
#' This function takes in an object returned by \code{optim} and returns logicals on whether the corresponding 95%
#' confidence interval of each estimate contains the true value.
#' 
#' @param opt An optim object
#' @param trueParam A vector containing the true parameters
#' @return A vector of indicators, 1 indicating that the true parameter was included in the corresponding confidence interval
getCover <- function(opt, trueParam){
  sd <- sqrt(diag(solve(opt$hessian)))
  est <- opt$par
  as.numeric( est - 1.96*sd < trueParam)*as.numeric(est + 1.96*sd > trueParam )  #multiply numerics so that its in both interval endpoints
}

#' Get coefficient estimate from optim object
#' 
#' @param opt An optim object
#' @return The estimates contained in the optim object
getEst <- function(opt){
  opt$par
}

#main FM function, generates new fake data, returns a list with each entry an optim object
#each list entry corresponds to an element u from tList

#' Main function for generating observations from birth-shift-death process and inferring parameters under frequent monitoring
#' 
#' This function generates observation intervals from the simple birth-shift-death process without covariates, and generates
#' a dataset for each observation time length in tList. 
#' Next, the frequent monitoring log-likelihood is maximized using \code{\link{FM.optim}} for each dataset.
#' 
#' See accompanying vignette for example usage.
#' 
#' @inheritParams makedata.simple
#' @param initGuess Vector of numbers, initial guess for \code{optim}
#' @return A list of optim objects
FM.run <- function(N, tList, lam, v, mu,initList, initGuess){
  simDataList <- makedata.simple(N,tList,lam,v,mu,initList)  #simulate data
  result <- vector("list",length(tList))
  for(u in 1:length(tList)){
    result[[u]] <- FM.optim(simDataList, u, initGuess)
  }
  return(result)
}

#' Replicate \code{\link{FM.run}}
#' 
#' @param numReps The number of desired replications
#' @inheritParams FM.run
#' @return An array with entries of type returned by \code{\link{FM.run}}. Rows correspond to a dt value in tList, and
#' columns correspond to replications
#' 
#' @examples 
#' tList <- c(.2,.4,.6); initList <- c(1:15)
#' lam = .06; v = .02; mu = .11
#' trueParam <- c(lam,v,mu) 
#' N = 50; numReps = 2
#' example <- FM.replicate(numReps,N,tList,lam,v,mu,initList, trueParam)
FM.replicate <- function(numReps,N,tList,lam,v,mu,initList, initGuess){
  replicate(numReps, FM.run(N,tList,lam,v,mu,initList, initGuess))
}

### Repeat these for FFT method ###

#' Calculate log-likelihood using generating function technique 
#' 
#' \code{logFFT} calculates the log-likelihood where transition probabilities are computed using our FFT generating
#' function method
#' 
#' @param param A vector of three numbers containing the birth, shift, and death rate respectively
#' @param u A number, the observation interval length
#' @param simData An matrix in the format of a list entry returned by \code{\link{makedata.simple}}
#' @param initList A vector of possible initial populations
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @return The negative log-likelihood value
logFFT <- function(param, u, simData, initList, s1.seq, s2.seq){
  lam <- param[1]; v = param[2]; mu = param[3]
  tpm.list <- getTrans.initList(u, initList, lam, v, mu, s1.seq, s2.seq,.01) #this stores transition matrices for all initial sizes
  logresult <- 0
  for(i in 1:dim(simData)[2]){ #fake is a global holding the fake data
    initial <- simData[1,i]
    n.old <- simData[2,i]
    n.new <- simData[3,i]
    logresult <- logresult + log(tpm.list[[initial]][n.old+1,n.new+1])
  }
  return(-logresult)
}

#' Maximize the log-likelihood in \code{\link{logFFT}}
#' 
#' \code{FFT.optim} maximizes the log-likelihood using \code{optim}, with hessian = TRUE.
#' 
#' @param initGuess A vector containing the initial guess of birth, shift, and death rates
#' @param param A vector of three numbers containing the birth, shift, and death rate respectively
#' @param u A number, the observation interval length
#' @param simDataList A list of matrices in the format returned by \code{\link{makedata.simple}}
#' @param initList A vector of possible initial populations
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @return An optim type object 
FFT.optim <- function(simDataList, u, initGuess, initList, s1.seq, s2.seq){
  optim(initGuess, logFFT, u=tList[u], simData = simDataList[[u]], initList = initList, 
        s1.seq = s1.seq, s2.seq = s2.seq, hessian=T, control = list(trace = 6))
}

#' Generate synthetic data from simple BSD process and infer rates using generating function method
#' 
#' The main function for simulation studies assessing our generating function approach in the the discretely observed
#' simple birth-shift-death-process without covariates. Generates synthetic datasets using \code{\link{makedata.simple}}, and
#' infers the MLE rates for each using \code{\link{FFT.optim}}.
#' 
#' @inheritParams makedata.simple
#' @inheritParams FFT.optim
#' @return A list of optim objects
FFT.run <- function(N, tList, lam, v, mu, initList, initGuess, s1.seq, s2.seq ){
  simDataList <- makedata.simple(N,tList,lam,v,mu,initList)  #simulate data
  result <- vector("list",length(tList))
  for(u in 1:length(tList)){
    result[[u]] <- FFT.optim(simDataList, u, initGuess, initList, s1.seq, s2.seq)
  }
  return(result)
}

#' Replicate the function \code{\link{FFT.run}}
#' 
#' This function does the same as \code{\link{FM.replicate}}, but optimizes the likelihood based on transition
#' probabilities computed by the generating function method instead of using frequent monitoring. An example
#' of usage is included in the vignette.
#' 
#' @param numReps The number of replications
#' @inheritParams FFT.run
#' @return An array of optim objects in the same layout as \code{\link{FM.replicate}}
FFT.replicate <- function(numReps, N, tList, lam, v, mu, initList, initGuess, s1.seq, s2.seq){
  replicate(numReps, FFT.run(N,tList,lam,v,mu,initList,initGuess, s1.seq, s2.seq))
}


#' Perform one E-step of the EM algorithm for unevenly spaced observations
#' 
#' \code{ESTEP.realdata} performs one E-step of the EM algorithm similarly to \code{\link{ESTEP}},
#' but for datasets with unevenly spaced time intervals between observations. In particular, the second
#' row of the PATIENTDATA argument, a matrix returned by \code{\link{MakePatientData}} for simulated data
#' contains no relevant information for the algorithm. When working with real data, we may create a matrix
#' of the same format, except replace this row to instead contain information about the time between observations.
#' 
#' 
#' @param betaVec A vector, the setting of beta coefficients
#' @param num.patients An integer, number of unique patients
#' @param PATIENTDATA A matrix in the form returned by \code{\link{MakePatientData}} containing the set of observation intervals,
#' but the second row now contains the observation interval lengths 
#' @param patients.design A design matrix in the same form as returned by \code{\link{PatientDesignExample}}
#' @param s1.seq A vector of complex arguments evenly spaced along the unit circle
#' @param s2.seq A vector of complex arguments evenly spaced along the unit circle
#' @return A list containing a matrix of the expected sufficient statistics as well as the observed log likelihood value
ESTEP.realdata <- function(betaVec, num.patients, PATIENTDATA, patients.design, s1.seq, s2.seq){  
  bLam <- betaVec[1:(length(betaVec)/3)]; bNu = betaVec[(length(betaVec)/3+1):(2*length(betaVec)/3)]; bMu = betaVec[(2*length(betaVec)/3 + 1) : length(betaVec)]
  U <- D <- V <- P <- rep(0, num.patients) #initialize vectors of expected sufficient statistics
  observedLikelihood <- 0
  for(i in 1:num.patients){
    #get lam, v, mu for current patient
    lam <- exp( sum(bLam*patients.design[,i]))
    v <- exp( sum(bNu*patients.design[,i])) 
    mu <- exp(sum(bMu*patients.design[,i]))
    #this takes the subset of the fake data corresponding to patient i
    dat.i <- as.matrix(PATIENTDATA[,which(PATIENTDATA[1,]==UIDs[i])])
    #loop over number of observations associated with current patient:
    for(j in 1:dim(dat.i)[2]){
      #these three items are indices to look up terms in our transition lists
      n.initial <- dat.i[3,j]  #now this is actually initial number of particles, not an index
      n.old <- dat.i[4,j]
      n.new <- dat.i[5,j]
      timeInterval <- dat.i[2,j] #interval length
      #set the dt for numerical solver
      dt = timeInterval
      if( n.initial == n.old & n.new == 0){
        #print('no event in this interval')
        #zero expected counts to add, but particle time needs to be accounted for
        P[i] <- P[i] + n.initial*timeInterval  #particle time is just initial particles multiplied by interval length in row 2
        observedLikelihood <- observedLikelihood -n.initial*(lam+v+mu)*timeInterval
      } else {
        tpmMat <- getTrans(timeInterval, n.initial, lam, v, mu, s1.seq, s2.seq, dt)
        restrictedBirthMat <- getBirthMeans(timeInterval, n.initial, lam, v, mu, s1.seq, s2.seq, dt)
        restrictedDeathMat <- getDeathMeans(timeInterval, n.initial, lam, v, mu, s1.seq, s2.seq, dt)
        restrictedShiftMat <- getShiftMeans(timeInterval, n.initial, lam, v, mu, s1.seq, s2.seq, dt)
        restrictedParticleMat <- getParticleT(timeInterval, n.initial, lam, v, mu, s1.seq, s2.seq, dt)       
        #we divide all the restricted moments by the corresponding transition probability
        #we add to the i'th entry for each interval associated with patient i in this loop
        observedLikelihood <- observedLikelihood + log(tpmMat[n.old+1, n.new+1])
        U[i] <- U[i] + restrictedBirthMat[n.old+1, n.new+1]/tpmMat[n.old+1, n.new+1]
        D[i] <- D[i] + restrictedDeathMat[n.old+1, n.new+1]/tpmMat[n.old+1, n.new+1]
        V[i] <- V[i] + restrictedShiftMat[n.old+1, n.new+1]/tpmMat[n.old+1, n.new+1]
        P[i] <- P[i] + restrictedParticleMat[n.old+1, n.new+1]/tpmMat[n.old+1, n.new+1]
      }
    }
  }
  result <- rbind(U,D,V,P)
  result[result < 0] = .000001 #numerical stability for negatives
  resultList <- list( 'matrix' = result, 'loglike' = observedLikelihood)
  return(resultList)  #return the vectors of expected sufficient statistics
}
jasonxu90/bdsem documentation built on May 18, 2019, 5:54 p.m.