tests/example2nd.R

In TIMP: Fitting Separable Nonlinear Models in Spectroscopy and Microscopy

## Example in modeling second order kinetics, by
## David Nicolaides.

## On simulated data.

##############################
## load TIMP
##############################

library("TIMP")

##############################
## SIMULATE DATA
##############################

## set up the Example problem, a la in-situ UV-Vis spectroscopy of a simple
## reaction.
## A + 2B -> C + D, 2C -> E

cstart <- c(A = 1.0, B = 0.8, C = 0.0, D = 0.0, E = 0.0)
times <- c(seq(0, 2, length = 21), seq(3, 10, length = 8))
k <- c(kA = 0.5, k2C = 1)

## stoichiometry matrices as per
## Puxty, G., Maeder, M., and Hungerbuhler, K. (2006) Tutorial on the fitting
## of kinetics models to multivariate spectroscopic measurements with
## non-linear least-squares regression, Chemometrics and Intelligent
## Laboratory Systems 81, 149-164.

rsmatrix <- c(1, 2, 0, 0, 0, 0, 0, 2, 0, 0)
smatrix <- c(-1, -2, 1, 1, 0, 0, 0, -2, 0, 1)
concentrations <- calcD(k, times, cstart, rsmatrix, smatrix)

wavelengths <- seq(500, 700, by = 2)
spectra <- matrix(nrow = length(wavelengths), ncol = length(cstart))
location <- c(550, 575, 625, 650, 675)
delta <- c(10, 10, 10, 10, 10)
spectra[, 1] <- exp(-log(2) * (2 * (wavelengths - location[1]) / delta[1])^2)
spectra[, 2] <- exp(-log(2) * (2 * (wavelengths - location[2]) / delta[2])^2)
spectra[, 3] <- exp(-log(2) * (2 * (wavelengths - location[3]) / delta[3])^2)
spectra[, 4] <- exp(-log(2) * (2 * (wavelengths - location[4]) / delta[4])^2)
spectra[, 5] <- exp(-log(2) * (2 * (wavelengths - location[5]) / delta[5])^2)

sigma <- .001
Psi_q <- concentrations %*% t(spectra) + sigma *
  rnorm(dim(concentrations)[1] * dim(spectra)[1])

## store the simulated data in an object of class "dat"
kinetic_data <- dat(
  psi.df = Psi_q, x = times, nt = length(times),
  x2 = wavelengths, nl = length(wavelengths)
)

##############################
## DEFINE MODEL
##############################

## starting values
kstart <- c(kA = 1, k2C = 0.5)

## model definition for 2nd order kinetics
kinetic_model <- initModel(
  mod_type = "kin", seqmod = FALSE, kinpar = kstart,
  numericalintegration = TRUE,
  initialvals = cstart,
  reactantstoichiometrymatrix = rsmatrix,
  stoichiometrymatrix = smatrix
)

##############################
## FIT INITIAL MODEL
## adding constraints to non-negativity of the
## spectra via the opt option nnls=TRUE
##############################

kinetic_fit <- fitModel(
  data = list(kinetic_data), modspec = list(kinetic_model),
  opt = kinopt(
    nnls = TRUE, iter = 80,
    selectedtraces = seq(1, kinetic_data@nl, by = 2)
  )
)

## look at estimated parameters

parEst(kinetic_fit)

## make a png of various results

## concentrations

conRes <- getX(kinetic_fit)

matplot(times, conRes,
  type = "b", col = 1, pch = 21, bg = 1:5, xlab = "time (sec)",
  ylab = "concentrations", main = "Concentrations (2nd order kinetics)"
)


## spectra

specRes <- getCLP(kinetic_fit)

matplot(wavelengths, specRes,
  type = "b", col = 1, pch = 21, bg = 1:5,
  xlab = "wavelength (nm)",
  ylab = "amplitude", main = "Spectra"
)


## see help(getResults) for how to get more results information from
## kinetic_fit

##############################
## CLEANUP GENERATED FILES
##############################
# This removes the files that were generated in this example
# (do not run this code if you wish to inspect the output)
file_list_cleanup <- c(Sys.glob("*paramEst.txt"), Sys.glob("*.ps"), Sys.glob("Rplots*.pdf"))

# Iterate over the files and delete them if they exist
for (f in file_list_cleanup) {
  if (file.exists(f)) {
    unlink(f)
  }
}