#---------------------------------
#Run this part to generate the data frame to be processed by the function
dosetimes <- c(0,12)
sampletimes <- sort(unique(c(seq(from=0,to=100,by=1))))
#number of subjects
ID <- 1:10
#Make dataframe
df <- expand.grid("ID"=ID,"TIME"=sampletimes, "AMT"=0)
doserows <- subset(df, TIME%in%dosetimes)
doserows$AMT <- 100
#Add back dose information
df <- rbind(df,doserows)
df <- df[order(df$ID,df$TIME),]
df <- subset(df, (TIME==0 & AMT==0)==F)
df$CL[df$ID<=1] <- 2
df$CL[df$ID>=2] <- 2.5
df$V2[df$ID<=1] <- 10
df$V2[df$ID>=2] <- 15
df$Q [df$ID<=1] <- 20
df$Q [df$ID>=2] <- 20
df$V3[df$ID<=1] <- 30
df$V3[df$ID>=2]  <- 30
df$KA[df$ID<=1] <- 0.2
df$KA[df$ID>=2]  <- 0.25
df$F1[df$ID<=1] <- 1
df$F1[df$ID>=2]  <- 0.5

#calculate micro-parameters
df$k20 <- df$CL/df$V2
df$k23 <- df$Q/df$V2
df$k32 <- df$Q/df$V3

#-----------------------------------
# This is the function to provess the data frame
TwoCompOral <- function(d){
  #Accepts a NONMEM style data frame for 1 subject with columns for TIME, AMT,MDV,DV, CL, V2, Q, V3, KA & F1
  #Returns a dataframe with populated columns for A1, A2, A3 and DV

  #set initial values in the compartments
  d$A1[d$TIME==0] <- d$AMT[d$TIME==0]*d$F1[1]  # Amount in the absorption compartment at time zero.
  d$A2[d$TIME==0] <- 0                      # Amount in the central compartment at time zero.
  d$A3[d$TIME==0] <- 0                      # Amount in the peripheral compartment at time zero.

  #This loop calculates micro-rate constants based on individual's PK parameter values.
  #It uses these values to calculate macro-rate constants (Lambda1/lambda2).
  #Rate constants(micro- and macro), along other parameters, are used in the equations to calculate drug amounts in each compartment.
  #The loop advances the solution from one time interval to the next.
  k20 <- d$k20[1]
  k23 <- d$k23[1]
  k32 <- d$k32[1]
  KA <- d$KA[1]
  k30 <- 0
  E2 <- k20+k23
  E3 <- k32+k30

  #calculate hybrid rate constants
  lambda1 = 0.5*((E2+E3)+sqrt((E2+E3)^2-4*(E2*E3-k23*k32)))
  lambda2 = 0.5*((E2+E3)-sqrt((E2+E3)^2-4*(E2*E3-k23*k32)))


  for(i in 2:nrow(d))
  {

    t <- d$TIME[i]-d$TIME[i-1]
    A2last <- d$A2[i-1]
    A3last <- d$A3[i-1]
    A1last <- d$A1[i-1]

    A2term1 = (((A2last*E3+A3last*k32)-A2last*lambda1)*exp(-t*lambda1)-((A2last*E3+A3last*k32)-A2last*lambda2)*exp(-t*lambda2))/(lambda2-lambda1)
    A2term2 = A1last*KA*(exp(-t*KA)*(E3-KA)/((lambda1-KA)*(lambda2-KA))+exp(-t*lambda1)*(E3-lambda1)/((lambda2-lambda1)*(KA-lambda1))+exp(-t*lambda2)*(E3-lambda2)/((lambda1-lambda2)*(KA-lambda2)))
    d$A2[i] = A2term1+A2term2  #Amount in the central compartment

    A3term1 = (((A3last*E2+A2last*k23)-A3last*lambda1)*exp(-t*lambda1)-((A3last*E2+A2last*k23)-A3last*lambda2)*exp(-t*lambda2))/(lambda2-lambda1)
    A3term2 = A1last*KA*k23*(exp(-t*KA)/((lambda1-KA)*(lambda2-KA))+exp(-t*lambda1)/((lambda2-lambda1)*(KA-lambda1))+exp(-t*lambda2)/((lambda1-lambda2)*(KA-lambda2)))
    d$A3[i] = A3term1+A3term2  #Amount in the peripheral compartment

    A1last = A1last*exp(-t*KA)
    d$A1[i] = A1last + d$AMT[i]*d$F1[1]  #Amount in the absorption compartment

    d$DV[i] <- d$A2[i]/d$V2[i]        #Concentration in the central compartment

  }
  d
}

#Apply the function to each ID
system.time(simdf <- ddply(df, .(ID), TwoCompOral))