R/getResults.R

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

## all of these functions act on the output of the fitModel function to
## return various results
getC <- function(result, dataset = 1, file = "") {
  C <- getX(result, dataset = dataset, file = "", lreturnC = TRUE)
  if (file != "") {
    write.table(C,
      file = paste(file,
        "_C_dataset_", dataset, ".txt",
        sep = ""
      ),
      row.names = result$currModel@modellist[[dataset]]@x,
      quote = FALSE
    )
  }
  C
}
getCLPList <- function(result, getclperr = FALSE, file = "") {
  specList <- getSpecList(result$currModel, result$currTheta, getclperr)
  if (file != "") {
    for (i in 1:length(specList)) {
      write.table(specList[[i]],
        file = paste(file,
          "_spec_dataset_", i, ".txt",
          sep = ""
        ),
        row.names = result$currModel@modellist[[i]]@x2,
        quote = FALSE
      )
    }
  }
  specList
}
getCLP <- function(result, getclperr = FALSE, dataset = 1, file = "") {
  spec <- getSpecList(result$currModel, result$currTheta, getclperr)[[dataset]]
  if (file != "") {
    write.table(spec,
      file = paste(file,
        "_spec_dataset_", dataset, ".txt",
        sep = ""
      ),
      row.names = result$currModel@modellist[[dataset]]@x2,
      quote = FALSE
    )
  }
  spec
}
getData <- function(result, dataset = 1, weighted = FALSE) {
  if (weighted) {
    datamat <- result$currModel@data[[dataset]]@psi.weight
  } else {
    datamat <- result$currModel@data[[dataset]]@psi.df
  }
  datamat
}
getDAS <- function(result, getclperr = FALSE, dataset = 1, file = "") {
  spec <- getSpecList(result$currModel, result$currTheta, getclperr)[[dataset]]
  if ((result$currModel@modellist[[dataset]]@fullk == TRUE) ||
    (result$currModel@modellist[[dataset]]@seqmod == TRUE)) {
    A <- getX(result, dataset = dataset, file = "", lreturnA = TRUE)
    ncomp <- nrow(A)
    specA <- spec[, 1:ncomp] %*% t(A)
    spec[, 1:ncomp] <- specA
  }
  if (file != "") {
    write.table(spec,
      file = paste(file,
        "_DAS_dataset_", dataset, ".txt",
        sep = ""
      ),
      row.names = result$currModel@modellist[[dataset]]@x2,
      quote = FALSE
    )
  }
  spec
}
getSVDData <- function(result, numsing = 2, dataset = 1) {
  datamat <- getData(result, dataset)
  doSVD(datamat, numsing, numsing)
}
getResiduals <- function(result, dataset = 1) {
  residlist <- result$currModel@fit@resultlist[[dataset]]@resid
  residmat <- unlist(residlist)
  dim(residmat) <- c(length(residlist[[1]]), length(residlist))
  residmat
}
getSVDResiduals <- function(result, numsing = 2, dataset = 1) {
  residmat <- getResiduals(result, dataset)
  doSVD(residmat, numsing, numsing)
}
getTraces <- function(result, dataset = 1, file = "") {
  fitted <- result$currModel@fit@resultlist[[dataset]]@fitted
  tracesmat <- unlist(fitted)
  dim(tracesmat) <- c(length(fitted[[1]]), length(fitted))
  if (file != "") {
    write.table(tracesmat,
      file = paste(file,
        "fit.txt",
        sep = ""
      ),
      col.names = result$currModel@modellist[[dataset]]@x2,
      row.names = result$currModel@modellist[[dataset]]@x,
      quote = FALSE
    )
  }
  tracesmat
}
getdim1 <- function(result, dataset = 1) {
  result$currModel@modellist[[dataset]]@x
}
getdim2 <- function(result, dataset = 1) {
  result$currModel@modellist[[dataset]]@x2
}
parEst <- function(result, param = "", dataset = NA, verbose = TRUE,
                   file = "", stderr = TRUE) {
  currTheta <- result$currTheta
  currModel <- result$currModel
  if (stderr && currModel@optlist[[1]]@sumnls) {
    stderr <- TRUE
    if (is(currModel@fit@nlsres$onls, "timp.nls.lm")) {
      currErr <- getThetaCl(sumnls(result)$coefficients[, 2], currModel)
    } else {
      currErr <- getThetaCl(sumnls(result)$parameters[, 2], currModel)
    }
  } else {
    stderr <- FALSE
  }
  reslist <- stderrlist <- list()
  if (param == "") {
    param <- slotNames(theta())
  }
  if (is.na(dataset)) {
    dataset <- 1:length(currTheta)
  }

  for (nm in param) {
    for (j in dataset) {
      if (length(slot(currTheta[[j]], nm)) > 0) {
        if (is.null(reslist[[nm]])) {
          reslist[[nm]] <- stderrlist[[nm]] <- list()
          if (verbose) cat("Parameters:", nm, "\n", file = file)
        }
        reslist[[nm]][[length(reslist[[nm]]) + 1]] <- slot(currTheta[[j]], nm)
        if (stderr) {
          if (nm %in% currModel@modellist[[j]]@positivepar) {
            std <- stdErrTransform(
              slot(currTheta[[j]], nm),
              slot(currErr[[j]], nm)
            )
          } else {
            std <- slot(currErr[[j]], nm)
          }
          stderrlist[[nm]][[length(stderrlist[[nm]]) + 1]] <- std
        }
        if (verbose) {
          cat("dataset ", j, ": ", toString(slot(currTheta[[j]], nm)),
            "\n",
            sep = "", file = file
          )
          if (stderr) {
            cat("    standard errors: ",
              toString(stderrlist[[nm]][[length(stderrlist[[nm]])]]),
              "\n",
              sep = "", file = file
            )
          }
        }
      }
    }
  }
  invisible(list(reslist = reslist, stderrlist = stderrlist))
}
stdErrTransform <- function(k, err) {
  x <- rep(0, length(err))
  kl <- log(k)
  errl <- log(err)
  for (i in 1:length(k)) {
    b1 <- abs(exp(kl[i] + errl[i]) - exp(kl[i]))
    b2 <- abs(exp(kl[i] - errl[i]) - exp(kl[i]))
    x[i] <- ifelse(b1 > b2, b1, b2)
  }
  x
}


sumnls <- function(result) {
  result$currModel@fit@nlsres[[2]]
}
onls <- function(result) {
  result$currModel@fit@nlsres[[1]]
}