## model
dM1 <- function(dt, M1, M2, e01, k12, k21, k10)
  dt*(e01 - k12*M1 + k21*M2 - k10*M1)
dM2 <- function(dt, M1, M2, k12, k21)
  dt*(    + k12*M1 - k21*M2)

## `sol2C` implements the Euler method of
## numerical integration, where h is the step size
## in time, and f(t, y) is the differential equation dy/dt:
## y_n = y_n-1 + h*f(t_n-1, y_n-1)
sol2C <- function(e01, k12=2, k21=2, k10=0.7, dt=1/60) {
  M1 <- M2 <- time <- rep(0, length(e01))
  for(t in 1:(length(e01)-1)) {
    d1 <- dM1(dt, M1[t], M2[t], e01[t], k12, k21, k10)
    d2 <- dM2(dt, M1[t], M2[t], k12, k21)
    M1[t+1] <- M1[t] + d1
    M2[t+1] <- M2[t] + d2
    time[t+1] <- time[t] + dt
  }
  ret <- data.frame(M1=M1, M2=M2, e01=e01, time=time)
  class(ret) <- c("sol2C", class(ret))
  return(ret)
}


## plot the results of `sol2C`
plot.sol2C <- function(obj) {
  par(mfrow=c(2,1), mai=rep(0,4), omi=rep(0.75,4))
  plot(obj$time, obj$M1, type="l",
       ylim=range(obj$M1, obj$M2),
       ylab="Mass",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("Mass", 2, 1, xpd=NA)
  lines(obj$time, obj$M2, lty=2)
  abline(h=0, lty=2, col="gray")
  legend("topright", legend=c("M1", "M2"),
         lty=c(1,2), cex=0.7, bty="n")
  plot(obj$time, obj$e01, type="l",
       ylab="e01",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("e01", 2, 1, xpd=NA)
  abline(h=0, lty=2, col="gray")
}

## stimulus functions
dt   <- 1/60 
time <- seq(0,8,dt)
## single IV infusion
e01_1iv <- ifelse(time <= 30/60, 1, 0)
## double IV infusion
e01_2iv <- ifelse(time <= 30/60 | 
                 (time > 60/60 & time <= 90/60), 1, 0)
## single IV followed by removal
e01_1iv1rm <- ifelse(time <= 30/60, 1, 0)
e01_1iv1rm <- ifelse(time > 60/60 & time <= 90/60,
                     -0.45, e01_1iv1rm)
## single IV followed by removal then additional IV
e01_2iv1rm <- ifelse(time <= 30/60, 1, 0)
e01_2iv1rm <- ifelse(time > 60/60 & time <= 90/60,
                     -0.45, e01_2iv1rm)
e01_2iv1rm <- ifelse(time > 120/60 & time <= 150/60,
                     1.0, e01_2iv1rm)
## sawtooth
e01_saw <- ifelse(time <= 30/60, time, 0)
e01_saw <- ifelse(time > 60/60 & time <= 90/60,
                  time-60/60, e01_saw)


## compute the Jacobian using a numerical derivative
sen2C <- function(e01, k12=2, k21=2, k10=0.7, dt=1/60) {
  sol0 <- sol2C(e01, k12, k21, k10, dt)
  delk12 <- sqrt(.Machine$double.eps)*k12
  solk12 <- sol2C(e01, k12+delk12, k21, k10, dt)
  delk21 <- sqrt(.Machine$double.eps)*k21
  solk21 <- sol2C(e01, k12, k21+delk21, k10, dt)
  delk10 <- sqrt(.Machine$double.eps)*k10
  solk10 <- sol2C(e01, k12, k21, k10+delk10, dt)
  ret <- data.frame(M1=sol0$M1, M2=sol0$M2,
                    k12=(solk12$M1-sol0$M1)/delk12,
                    k21=(solk21$M1-sol0$M1)/delk21,
                    k10=(solk10$M1-sol0$M1)/delk10,
                    e01=e01,
                    time=sol0$time)
  class(ret) <- c("sen2C", class(ret))
  return(ret)
}

## compute the parameter correlations
cor2C <- function(sen, times=c(seq(0.25, 3, 0.25),4,6,8)) {
  senm <- as.matrix(sen[,c("k12","k21","k10")])
  covm <- t(senm)%*%senm
  return(cov2cor(covm))
}

## plot the model solutions and Jacobian
plot.sen2C <- function(obj) {
  par(mfrow=c(5,1), mai=rep(0,4), omi=rep(0.5,4))
  plot(obj$time, obj$M1, type="l",
       ylim=range(obj$M1, obj$M2),
       ylab="Mass",
       xaxt="n", yaxt="n")
  mtext("Mass", 2, 1, xpd=NA, cex=0.8)
  lines(obj$time, obj$M2, lty=2)
  abline(h=0, lty=2, col="gray")
  legend("topright", legend=c("M1", "M2"),
         lty=c(1,2), bty="n")
  plot(obj$time, obj$e01, type="l",
       ylab="e01",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("e01", 2, 1, xpd=NA, cex=0.8)
  abline(h=0, lty=2, col="gray")
  plot(obj$time, obj$k12, type="l",
       ylab="dM1/dk12",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("dM1/dk12", 2, 1, xpd=NA, cex=0.8)
  abline(h=0, lty=2, col="gray")
  plot(obj$time, obj$k21, type="l",
       ylab="dM1/dk21",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("dM1/dk21", 2, 1, xpd=NA, cex=0.8)
  abline(h=0, lty=2, col="gray")
  plot(obj$time, obj$k10, type="l",
       ylab="dM1/dk10",
       xlab="Time",
       xaxt="n", yaxt="n")
  mtext("dM1/dk10", 2, 1, xpd=NA, cex=0.8)
  mtext("Time", 1, 1, xpd=NA, cex=0.8)
  abline(h=0, lty=2, col="gray")
}

## plot the solutions
plot(sol2C(e01_1iv))
plot(sol2C(e01_2iv))
plot(sol2C(e01_saw))
plot(sol2C(e01_1iv1rm))
plot(sol2C(e01_2iv1rm))

## plot the solutions with Jacobian
plot(sen2C(e01_1iv))
plot(sen2C(e01_2iv))
plot(sen2C(e01_saw))
plot(sen2C(e01_1iv1rm))
plot(sen2C(e01_2iv1rm))
     
## show the correlations
cor2C(sen2C(e01_1iv))
cor2C(sen2C(e01_2iv))
cor2C(sen2C(e01_saw))
cor2C(sen2C(e01_1iv1rm))
cor2C(sen2C(e01_2iv1rm))