R/set_plotter_kin.R

In glotaran/TIMP: Fitting Separable Nonlinear Models in Spectroscopy and Microscopy

setMethod("plotter", signature(model = "kin"), function(model,
                                                        multimodel, multitheta, plotoptions) {
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  if (!plotoptions@notraces) {
    plotTracesSuper(multimodel, multitheta, plotoptions)
  }
  if (plotoptions@residplot) {
    if (!plotoptions@FLIM) {
      plotResids(multimodel, multitheta, plotoptions)
    }
  }
  if (plotoptions@residtraces) {
    plotTracesResids(multimodel, multitheta, plotoptions)
  }
  if (plotoptions@writefit || plotoptions@writefitivo) {
    writeFit(multimodel, multitheta, plotoptions)
  }
  if (length(plotoptions@breakdown) > 0) {
    plotKinBreakDown(multimodel, multitheta, plotoptions)
  }
  if (plotoptions@FLIM) {
    plotFLIM(multimodel, multitheta, plotoptions)
    return()
  }
  if (dev.cur() != 1) {
    dev.new()
  }
  par(
    mgp = c(2, 1, 0), mar = c(3, 3, 3, 2), oma = c(1, 0, 4, 0),
    mfrow = c(plotoptions@summaryplotrow, plotoptions@summaryplotcol)
  )
  resultlist <- multimodel@fit@resultlist
  m <- multimodel@modellist
  t <- multitheta
  groups <- multimodel@groups
  allx2 <- allx <- vector()
  for (i in 1:length(m)) {
    allx2 <- append(allx2, m[[i]]@x2)
    allx <- append(allx, m[[i]]@x)
  }
  xmax <- max(allx)
  xmin <- min(allx)
  x2max <- max(allx2)
  x2min <- min(allx2)
  conmax <- tauList <- muList <- contoplotList <- list()

  if (plotoptions@adddataimage) {
    for (i in 1:length(m)) {
      datarev <- t(apply(m[[i]]@psi.df, 2, rev))
      xr <- rev(m[[i]]@x)
      x2lab <- xr[seq(1, length(xr), length = 10)]
      x2at <- m[[i]]@x[seq(1, length(xr), length = 10)]
      image(
        x = m[[i]]@x2,
        y = m[[i]]@x,
        z = datarev,
        ylab = plotoptions@xlab,
        yaxt = "n",
        main = "Data",
        xlab = plotoptions@ylab,
        col = tim.colors()
      )
      axis(2, at = x2at, labels = x2lab)
    }
  }

  if (m[[i]]@anispec$useparperp) {
    calcAniSignal(m, plotoptions)
  }
  f1 <- function(x) {
    x[[1]]
  }
  f2 <- function(x) {
    x[[2]]
  }
  ## will plot the first concentration from each dataset
  grtoplot <- vector("list", length(m))
  for (i in 1:length(m)) {
    cnt <- 1
    notfound <- TRUE
    while (notfound) {
      for (j in 1:length(groups[[cnt]])) {
        if (groups[[cnt]][[j]][2] == i) {
          grtoplot[[i]] <- list(groups[[cnt]], j)
          notfound <- FALSE
        }
      }
      cnt <- cnt + 1
    }
  }
  for (i in 1:length(m)) {
    group <- grtoplot[[i]][[1]]
    place <- grtoplot[[i]][[2]]
    dset <- group[[place]][2]
    irfmu <- unlist(lapply(resultlist[[i]]@irfvec, f1))
    irftau <- unlist(lapply(resultlist[[i]]@irfvec, f2))
    muList[[i]] <- irfmu
    tauList[[i]] <- irftau
    if (plotoptions@writerawcon ||
      length(plotoptions@writeplaincon) > 0) {
      C <- getKinConcen(group, multimodel, t, doConstr = FALSE, oneDS = place)
      if (plotoptions@writerawcon) {
        write.table(C,
          file = paste(plotoptions@makeps,
            "_rawconcen_dataset_", dset, ".txt",
            sep = ""
          ),
          quote = FALSE,
          row.names = m[[dset]]@x
        )
      }
      if (length(plotoptions@writeplaincon) > 0) {
        xplot <- plotoptions@writeplaincon$x
        CWRITE <- C
        write.table(CWRITE,
          file = paste(plotoptions@makeps,
            "_plaincon_dataset_", dset, ".txt",
            sep = ""
          ),
          quote = FALSE,
          row.names = m[[dset]]@x
        )
      }
    }
    contoplotList[[length(contoplotList) + 1]] <- getConToPlot(getKinConcen(
      group, multimodel, t,
      oneDS = place
    ), m[[i]]@cohspec, m[[i]]@cohcol)
    conmax[[i]] <- attributes(contoplotList[[length(contoplotList)]])$max
  }
  minc <- do.call("min", contoplotList)
  maxc <- do.call("max", contoplotList)
  for (i in 1:length(m)) {
    matlinlogplot(m[[i]]@x, contoplotList[[i]], muList[[i]][1],
      plotoptions@linrange,
      type = "l", add = !(i == 1), lty = i, ylab = "concentration",
      xlab = plotoptions@xlab, main = "Concentrations",
      xlim = c(xmin, xmax), ylim = c(minc, maxc)
    )
    if (plotoptions@writecon) {
      write.table(contoplotList[[i]],
        file = paste(plotoptions@makeps,
          "_concen_dataset_", i, ".txt",
          sep = ""
        ), quote = FALSE,
        row.names = m[[i]]@x
      )
    }
  }
  if (length(model@kin2scal) != 0) {
    perA <- vector()
    for (i in 1:length(m)) {
      perA <- append(perA, t[[i]]@kin2scal[2])
    }
    matplot(1:length(m), perA,
      type = "l",
      main = "% Concentration Photoconverted",
      xlab = "Dataset number", ylab = "percent"
    )
  }
  spectralist <- getSpecList(multimodel, t)
  specList <- list()
  for (i in 1:length(m)) {
    if (length(conmax) > 0) {
      spec <- getSpecToPlot(
        spectralist[[i]], conmax[[i]],
        m[[i]]@cohcol, plotoptions@plotcohcolspec
      )
    } else {
      spec <- spectralist[[i]]
    }

    specList[[length(specList) + 1]] <- spec
  }
  maxs <- do.call("max", specList)
  mins <- do.call("min", specList)
  for (i in 1:length(m)) {
    matplot(m[[i]]@x2, specList[[i]],
      type = "l", main = "Spectra",
      xlab = plotoptions@ylab, ylab = "amplitude", lty = i,
      add = !(i == 1), xlim = c(x2min, x2max), ylim = c(mins, maxs)
    )
    if (plotoptions@writespec) {
      write.table(specList[[i]],
        file = paste(plotoptions@makeps,
          "_spec_dataset_", i, ".txt",
          sep = ""
        ),
        row.names = m[[i]]@x2, quote = FALSE
      )
    }
  }
  abline(0, 0)
  for (i in 1:length(m)) {
    matplot(m[[i]]@x2, normdat(specList[[i]]),
      type = "l",
      main = "Normalized spectra", xlim = c(x2min, x2max),
      xlab = plotoptions@ylab, ylab = "amplitude", lty = i,
      add = !(i == 1)
    )
    abline(0, 0)
    if (plotoptions@writenormspec) {
      write.table(normdat(specList[[i]]),
        file = paste(plotoptions@makeps,
          "_normspec_dataset_", i, ".txt",
          sep = ""
        ),
        row.names = m[[i]]@x2, quote = FALSE
      )
    }
  }
  abline(0, 0)
  notplotted <- TRUE
  minmu <- do.call("min", muList)
  maxmu <- do.call("max", muList)
  for (i in 1:length(m)) {
    if (m[[i]]@dispmu) {
      matplot(m[[i]]@x2, muList[[i]],
        type = "l",
        main = "IRF location", xlim = c(x2min, x2max),
        add = !notplotted, ylim = c(minmu, maxmu), col = i,
        xlab = plotoptions@ylab, ylab = "IRF location"
      )
      notplotted <- FALSE
    }
  }
  notplotted <- TRUE
  mintau <- do.call("min", tauList)
  maxtau <- do.call("max", tauList)
  for (i in 1:length(m)) {
    if (m[[i]]@disptau) {
      matplot(m[[i]]@x2, tauList[[i]],
        type = "l",
        main = "IRF width", add = !notplotted, col = i,
        xlim = c(x2min, x2max), ylim = c(mintau, maxtau),
        xlab = plotoptions@ylab, ylab = "IRF width"
      )
      notplotted <- FALSE
    }
  }
  svddatalist <- list()
  residlist <- svdresidlist <- list()
  for (i in 1:length(m)) {
    residuals <- matrix(nrow = m[[i]]@nt, ncol = m[[i]]@nl)
    for (j in 1:length(resultlist[[i]]@resid)) {
      residuals[, j] <- resultlist[[i]]@resid[[j]]
    }
    svdresidlist[[i]] <- doSVD(residuals, 2, 2)
    residlist[[i]] <- residuals
    svddatalist[[i]] <- doSVD(multimodel@data[[i]]@psi.df, 2, 2)
  }

  rangeleftd <- rangerightd <- rangevald <- 0
  rangeleftr <- rangerightr <- rangevalr <- 0
  for (i in 1:length(m)) {
    rangeleftr <- range(svdresidlist[[i]]$right[1, ], rangeleftr)
    rangerightr <- range(svdresidlist[[i]]$right[1, ], rangerightr)
    rangevalr <- range(svdresidlist[[i]]$values, rangevalr)

    rangeleftd <- range(svddatalist[[i]]$right[1, ], rangeleftd)
    rangerightd <- range(svddatalist[[i]]$right[1, ], rangerightd)
    rangevald <- range(svddatalist[[i]]$values, rangevald)
  }
  ## START RESID PLOTTING
  for (i in 1:length(m)) {
    limd <- max(max(residlist[[i]]), abs(min(residlist[[i]])))
    if (!(any(diff(m[[i]]@x) <= 0) || any(diff(m[[i]]@x2) <= 0))) {
      image.plot(m[[i]]@x, m[[i]]@x2,
        residlist[[i]],
        xlab = plotoptions@xlab,
        ylab = plotoptions@ylab,
        main = paste("Residuals Dataset", i),
        zlim = c(-limd, limd),
        col = diverge_hcl(40,
          h = c(0, 120), c = 60,
          l = c(45, 90), power = 1.2
        )
      )
    }
  }
  for (i in 1:length(m)) {
    if (m[[i]]@nt > 1 && m[[i]]@nl > 1) {
      xpos <- m[[i]]@x
      xpos[which(xpos <= 0)] <- NA
      matplot(xpos, svdresidlist[[i]]$left[, 1],
        type = "l", ylim = rangeleftr,
        main = "LSV Residuals", add = !(i == 1),
        log = "x", xlab = plotoptions@xlab, col = i,
        ylab = plotoptions@ylab
      )
    }
  }
  for (i in 1:length(m)) {
    if (m[[i]]@nt > 1 && m[[i]]@nl > 1) {
      matplot(m[[i]]@x2, svdresidlist[[i]]$right[1, ],
        type = "l", xlim = c(x2min, x2max),
        ylim = rangerightr,
        main = "RSV Residuals",
        xlab = plotoptions@xlab, add = !(i == 1),
        col = i, ylab = plotoptions@ylab
      )
    }
  }
  for (i in 1:length(m)) {
    if (i == 1) {
      plot(1:length(svdresidlist[[i]]$values),
        log10(svdresidlist[[1]]$values),
        xlab = "",
        ylab = plotoptions@ylab, col = i, ylim = rangevalr,
        main = "SVal Residuals", type = "b"
      )
    } else {
      lines(1:length(svdresidlist[[i]]$values),
        log10(svdresidlist[[i]]$values),
        type = "b", col = i
      )
    }
  }
  ## START DATA PLOTTING
  for (i in 1:length(m)) {
    if (m[[i]]@nt > 1 && m[[i]]@nl > 1) {
      xpos <- m[[i]]@x
      xpos[which(xpos <= 0)] <- NA
      matplot(xpos, svddatalist[[i]]$left[, 1],
        type = "l",
        main = "LSV Data", add = !(i == 1),
        log = "x", xlab = plotoptions@xlab, col = i,
        ylim = rangeleftd,
        ylab = plotoptions@ylab
      )
    }
  }
  for (i in 1:length(m)) {
    if (m[[i]]@nt > 1 && m[[i]]@nl > 1) {
      matplot(m[[i]]@x2, svddatalist[[i]]$right[1, ],
        type = "l",
        main = "RSV Data", xlim = c(x2min, x2max),
        ylim = rangerightd,
        xlab = plotoptions@xlab, add = !(i == 1),
        col = i, ylab = plotoptions@ylab
      )
    }
  }
  for (i in 1:length(m)) {
    if (i == 1) {
      plot(1:length(svddatalist[[i]]$values),
        log10(svddatalist[[1]]$values),
        xlab = "",
        ylab = plotoptions@ylab, col = i, ylim = rangevald,
        main = "SingVal Data", type = "b"
      )
    } else {
      lines(1:length(svddatalist[[i]]$values),
        log10(svddatalist[[i]]$values),
        type = "b", col = i
      )
    }
  }
  ## OTHER PLOTS
  for (i in 1:length(m)) {
    pl <- FALSE
    if (length(m[[i]]@dscalspec$perclp) != 0) {
      if (m[[i]]@dscalspec$perclp) {
        if (!pl) {
          plot(m[[i]]@x2, t[[i]]@drel,
            main = "Dataset scaling per clp",
            xlab = plotoptions@ylab, ylab = "", type = "l"
          )
          pl <- TRUE
        } else {
          lines(m[[i]]@x2, t[[i]]@drel,
            type = "l",
            col = i
          )
        }
      }
    }
  }
  if (length(plotoptions@title) != 0) {
    tit <- plotoptions@title
    if (plotoptions@addfilename) tit <- paste(tit, m[[i]]@datafile)
  } else {
    tit <- ""
    if (plotoptions@addfilename) tit <- paste(tit, m[[i]]@datafile)
  }
  mtext(tit, side = 3, outer = TRUE, line = 1)

  writeEst(multimodel, multitheta, plotoptions)
  # TODO: replace functionality provided by displayEst
  displayEst(plotoptions)

  if (dev.interactive() && length(plotoptions@makeps) != 0) {
    if (plotoptions@output == "pdf") {
      pdev <- pdf
    } else {
      pdev <- postscript
    }
    dev.print(
      device = pdev,
      file = paste(plotoptions@makeps,
        "_summary.", plotoptions@output,
        sep = ""
      )
    )
  }
  if (plotoptions@plotkinspec) {
    plotClp(multimodel, t, plotoptions)
  }
  if (plotoptions@kinspecest) {
    plotKinSpecEst(t, plotoptions, multimodel)
  }
})