R/kinfitr_miscfuncs.R

In mathesong/kinfitr: Kinetic Modelling of PET Time Activity Curves

#' Convolution Function
#'
#' Function to convolve two vectors of equal length for kinetic modelling.
#'
#' @param a A numeric vector of TAC function at equally spaced intervals
#' @param b A numeric vector of TAC function at equally spaced intervals of equal length to \code{a}
#' @param stepsize The time interval between each value of \code{a} and \code{b}
#'
#' @return A numeric vector of the convolved function of the same length as \code{a} and \code{b}
#'
#' @examples
#' a <- rnorm(20)
#' b <- rnorm(20)
#' stepsize <- 1
#'
#' kinfit_convolve(a, b, stepsize)
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

kinfit_convolve <- function(a, b, stepsize) {
  if (length(a) != length(b)) {
    stop("The vectors are not of the same length")
  }

  out <- pracma::conv(a, b)[1:length(a)] * stepsize
  return(out)
}

kinfit_uninterpbydur <- function(interptime, interppred, times, durtimes) {
  starttimes <- times - durtimes / 2
  interpstart <- sapply(starttimes, function(x) which.min(abs(interptime - x)))[-1]
  interpend <- c((interpstart - 1)[-1], length(interptime))

  startstop <- data.frame(start = interpstart, stop = interpend)

  predvals <- c(0, mapply(function(start, stop) mean(interppred[start:stop]), startstop$start, startstop$stop))

  return(predvals)
}

zerotrim_input <- function(t_inp, inp) {

  # Create a zero point by drawing a regression between the peak and the peak second differential

  input <- data.frame(Time = t_inp, Value = inp)

  input$Time <- input$Time

  input_start <- subset(input, input$Time < which.max(input$Value))
  # qplot(input_start$Time, input_start$Value)
  input_start$dif1 <- c(0, diff(input_start$Value))
  input_start$dif2 <- c(0, diff(input_start$dif1))
  lineseg <- input_start[which.max(input_start$dif2):which.max(input_start$Value), ]

  abc <- lm(lineseg$Value ~ lineseg$Time)

  # qplot(input_start$Time, input_start$Value) + geom_point() + geom_abline(slope = abc$coefficients[2], intercept=abc$coefficients[1])

  starttime <- -1 * abc$coefficients[1] / abc$coefficients[2]

  input <- subset(input, input$Time > starttime)
  input <- rbind(c(starttime, 0), input)

  input$Time <- input$Time - input$Time[1]

  return(input)
}


plot_fitpoints <- function(ymeasured, yfitted, xmeasured, xfitted = NULL) {
  if (is.null(xfitted)) {
    xfitted <- xmeasured
  }

  # plot the values
  plot(xmeasured, ymeasured)

  # draw the curve
  lines(xfitted, yfitted, col = "blue")
}

plot_fitfunc <- function(ymeasured, xmeasured, fitfunction, parameters) {

  # plot the values
  plot(xmeasured, ymeasured)

  # generate a range of values for xmeasured in small increments to create a smooth line
  xRange <- seq(min(xmeasured), max(xmeasured), length.out = 1024)

  # generate the predicted y values
  yValues <- predict(fitfunction, par = parameters, x = xRange)

  # draw the curve
  lines(xRange, yValues, col = "blue")
}


#' Plot Residuals of a Model
#'
#' Function to visualise the residuals of a model fit.  This function only works for model outputs for
#' which a fit object is available called \code{outputobject$fit}.
#'
#' @param outputobject The output of a kinetic model, including a fit object \code{outputobject$fit}
#'
#' @return A ggplot2 object of the residual plot
#'
#' @examples
#' # Note: Reference region models, and irreversible binding models, should not
#' # be used for PBR28 - this is just to demonstrate function
#'
#' data(pbr28)
#'
#' t_tac <- pbr28$tacs[[2]]$Times / 60
#' reftac <- pbr28$tacs[[2]]$CBL
#' roitac <- pbr28$tacs[[2]]$STR
#' weights <- pbr28$tacs[[2]]$Weights
#'
#' srtmout <- srtm(t_tac, reftac, roitac)
#'
#' plot_residuals(srtmout)
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

plot_residuals <- function(outputobject) {
  Time <- outputobject$tacs$Time
  Residuals <- residuals(outputobject$fit)
  Weights <- weights(outputobject$fit)

  if (length(Time) > length(Residuals)) {
    Time <- tail(Time, length(Residuals))
  }

  plotdf <- data.frame(Time = Time, Residuals = Residuals, Weights = Weights)

  ylimits <- c(-max(abs(plotdf$Residuals[plotdf$Weights > 0])), max(abs(plotdf$Residuals[plotdf$Weights > 0])))

  outplot <- ggplot(plotdf, aes(x = Time, y = Residuals)) +
    geom_point(aes(size = Weights)) +
    geom_hline(aes(yintercept = 0), linetype = "dashed") + ylim(ylimits) +
    scale_size(range = c(1, 3)) + geom_line()

  return(outplot)
}


#' Tidy Up for Reference Region Methods
#'
#' Function to tidy up the input argument vectors for reference region models.
#'
#' @param t_tac Numeric vector of times for each frame in minutes.
#' @param reftac Numeric vector of radioactivity concentrations in the reference tissue for each frame.
#' @param roitac Numeric vector of radioactivity concentrations in the target tissue for each frame.
#' @param weights Optional. Numeric vector of the weights assigned to each frame in the fitting.
#' @param frameStartEnd Optional. This allows one to specify the beginning and final frame to use for modelling, e.g. c(1,20).
#'
#' @return A dataframe containing the tidied up t_tac, reftac, roitac and weights ready for modelling.
#'
#' @details This function i) adds uniform weights if weights are not specified, ii) checks that the
#' lengths of t_tac, reftac, roitac and weights are of the same length, iii) shortens the vectors
#' if a frameStartEnd is specified, iv) adds a zero frame if there is not one, and v) checks that
#' times are in minutes and not in seconds.
#'
#' @examples
#'
#' #' # Note: Reference region models, and irreversible binding models, should not
#' # be used for PBR28 - this is just to demonstrate function
#'
#' data(pbr28)
#'
#' t_tac <- pbr28$tacs[[2]]$Times / 60
#' reftac <- pbr28$tacs[[2]]$CBL
#' roitac <- pbr28$tacs[[2]]$STR
#' weights <- pbr28$tacs[[2]]$Weights
#'
#' tidyinput_ref(t_tac, reftac, roitac, weights, frameStartEnd = c(1, 10))
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

tidyinput_ref <- function(t_tac, reftac, roitac, weights, frameStartEnd) {
  if (is.null(weights) == T) {
    weights <- rep(1, length(reftac))
  }

  lengths <- c(length(t_tac), length(reftac), length(roitac), length(weights))
  if (!all(lengths == lengths[1])) {
    stop("The lengths of t_tac, reftac, roitac and/or weights are not equal")
  }

  if (!is.null(frameStartEnd)) {
    t_tac <- t_tac[ frameStartEnd[1]:frameStartEnd[2] ]
    reftac <- reftac[ frameStartEnd[1]:frameStartEnd[2] ]
    roitac <- roitac[ frameStartEnd[1]:frameStartEnd[2] ]
    weights <- weights[ frameStartEnd[1]:frameStartEnd[2] ]
  }

  if (min(t_tac) > 0) {
    t_tac <- c(0, t_tac)
    roitac <- c(0, roitac)
    reftac <- c(0, reftac)
    weights <- c(0, weights)
  }

  if (max(t_tac) > 300) {
    warning("\n      ******************************************************************************\n          It looks like you have included seconds instead of minutes for time:\n          this can cause wrong/weird results, and should be avoided. Ignore this\n          warning if you just have really long measurements (over 300 minutes).\n      ******************************************************************************")
  }

  out <- data.frame(
    t_tac = t_tac, reftac = reftac, roitac = roitac,
    weights = weights
  )

  return(out)
}


#' Tidy Up for Models with Arterial Input
#'
#' Function to tidy up the input argument vectors for arterial input models.
#'
#' @param t_tac Numeric vector of times for each frame in minutes.
#' @param tac Numeric vector of radioactivity concentrations in the target tissue for each frame.
#' @param weights Optional. Numeric vector of the weights assigned to each frame in the fitting.
#' @param frameStartEnd Optional. This allows one to specify the beginning and final frame to use for modelling, e.g. c(1,20).
#'
#' @return A list containing the tidied up t_tac, tac and weights ready for modelling.
#'
#' @details This function i) adds uniform weights if weights are not specified, ii) checks that the
#' lengths of t_tac, tac and weights are of the same length, iii) shortens the vectors
#' if a frameStartEnd is specified, iv) Checks whether there are negative values in t_tac due to
#' adjusting for delay iv) adds a zero frame if there is not one, and v) checks that
#' times are in minutes and not in seconds.
#'
#' @examples
#'
#' data(pbr28)
#'
#' t_tac <- pbr28$tacs[[2]]$Times / 60
#' tac <- pbr28$tacs[[2]]$STR
#' weights <- pbr28$tacs[[2]]$Weights
#'
#' tidyinput_art(t_tac, tac, weights, frameStartEnd = c(1, 10))
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

tidyinput_art <- function(t_tac, tac, weights, frameStartEnd) {
  if (is.null(weights)) {
    weights <- rep(1, length(tac))
  }

  lengths <- c(length(t_tac), length(tac), length(weights))
  if (!all(lengths == lengths[1])) {
    stop("The lengths of t_tac, tac and/or weights are not equal")
  }

  if (!is.null(frameStartEnd)) {
    t_tac <- t_tac[ frameStartEnd[1]:frameStartEnd[2] ]
    tac <- tac[ frameStartEnd[1]:frameStartEnd[2] ]
    weights <- weights[ frameStartEnd[1]:frameStartEnd[2] ]
  }

  if (min(t_tac) < 0) {
    stop("There are negative times in the TAC")
  }

  if (min(t_tac) > 0) {
    t_tac <- c(0, t_tac)
    tac <- c(0, tac)
    weights <- c(0, weights)
  }

  if (max(t_tac) > 300) {
    warning("\n      ******************************************************************************\n          It looks like you have included seconds instead of minutes for time:\n          this can cause wrong/weird results, and should be avoided. Ignore this\n          warning if you just have really long measurements (over 300 minutes).\n      ******************************************************************************")
  }

  out <- data.frame(
    t_tac = t_tac, tac = tac,
    weights = weights
  )

  return(out)
}

#' Maximum Percentage Residual
#'
#' Function to determine the maximum percentage residual from the fitted value, i.e. the maximum
#' value of a residual as a percentage of its fitted value.  This is useful for obtaining estimates
#' of t*.
#'
#' @param outputobject The output of a kinetic model, including a fit object \code{outputobject$fit}
#'
#' @return The maximum percentage of the fitted value of a residual in the fit.
#'
#' @examples
#'
#' # Note: Reference region models, and irreversible binding models, should not
#' # be used for PBR28 - this is just to demonstrate function
#'
#' data(pbr28)
#'
#' t_tac <- pbr28$tacs[[2]]$Times / 60
#' reftac <- pbr28$tacs[[2]]$CBL
#' roitac <- pbr28$tacs[[2]]$STR
#' weights <- pbr28$tacs[[2]]$Weights
#'
#' refloganout <- refLogan(t_tac, reftac, roitac, 0.1, tstarIncludedFrames = 9)
#'
#' maxpercres(refloganout)
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

maxpercres <- function(outputobject) {
  Residuals <- abs(residuals(outputobject$fit))
  Fitted <- abs(fitted(outputobject$fit))
  PercRes <- 100 * Residuals / Fitted
  MaxPercRes <- max(PercRes)
  return(MaxPercRes)
}

#' Plot Kinetic Model Fit: Generic Function
#'
#' Function to plot the output of a kinetic model. This function calls the specific plotting functions
#' for each model based on the output of the model.
#'
#' @param modelout The output object of the model fitting procedure.
#' @param ... Additional optional arguments.
#'
#' @return A ggplot2 object of the plot.
#'
#' @details This function uses the \code{out$model} name to call the correct function to plot the model fit.
#'
#' @examples
#'
#' data(pbr28)
#'
#' t_tac <- pbr28$tacs[[2]]$Times / 60
#' tac <- pbr28$tacs[[2]]$FC
#' weights <- pbr28$tacs[[2]]$Weights
#'
#' input <- blood_interp(
#'   pbr28$procblood[[2]]$Time / 60, pbr28$procblood[[2]]$Cbl_dispcorr,
#'   pbr28$procblood[[2]]$Time / 60, pbr28$procblood[[2]]$Cpl_metabcorr,
#'   t_parentfrac = 1, parentfrac = 1
#' )
#'
#' fit <- ma1(t_tac, tac, input, 10, weights)
#'
#' plot_kinfit(fit, roiname = "FC")
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export

plot_kinfit <- function(modelout, ...) {
  modelname <- modelout$model
  input_list <- as.list(substitute(list(...)))
  do.call(what = paste0("plot_", modelname, "fit"), args = list(modelout, ...))
  # eval(expr = paste0('plot_', modelname, 'fit'))
}

#' Check and fix multstart parameters
#'
#' Checks and fixes for inputted multstart parameters. Checks missing multstart boundaries,
#' checks if some parameters aren't to be tried with different starting parameters, and checks
#' the lengths of the multstart bounds
#'
#' @param start Original starting values for without multstart
#' @param lower Lower fitting bounds
#' @param upper Upper fitting bounds
#' @param multstart_iter Iterations
#' @param multstart_lower Lower starting bounds
#' @param multstart_upper Upper starting bounds
#'
#' @return A list containing the parameters
#' @export
#'
#' @examples
#' start <- list(a = 1, b = 2, c = 3)
#' upper <- lapply(start, function(x) x + 1)
#' lower <- lapply(start, function(x) x - 1)
#' multstart_iter <- 5
#'
#' multstart_lower <- NULL
#' multstart_upper <- NULL
#'
#' fix_multstartpars(start, lower, upper, multstart_iter, multstart_lower, multstart_upper)
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
fix_multstartpars <- function(start, lower, upper, multstart_iter, multstart_lower, multstart_upper) {
  if (length(multstart_iter) == length(start) || length(multstart_iter) == 1) {
    multstart_l <- as.list(lower)
    multstart_u <- as.list(upper)

    parnames <- names(start)

    ### Adding multstart boundaries
    if (!is.null(multstart_lower)) {
      if (class(multstart_lower) != "list" ||
        sum(!(names(multstart_lower) %in% parnames)) > 0) { # Are there any names not in start?
        stop("multstart_lower should be a named list whose names match the parameters to be fitted")
      }

      multstart_lower <- modifyList(multstart_l, multstart_lower)
    } else {
      multstart_lower <- multstart_l
    }

    if (!is.null(multstart_upper)) {
      if (class(multstart_upper) != "list" ||
        sum(!(names(multstart_upper) %in% parnames)) > 0) { # Are there any names not in start?
        stop("multstart_upper should be a named list whose names match the parameters to be fitted")
      }

      multstart_upper <- modifyList(multstart_u, multstart_upper)
    } else {
      multstart_upper <- multstart_u
    }

    ### No multstart for some variables ###
    if (any(multstart_iter == 1)) {
      non_iterable <- which(multstart_iter == 1)
      multstart_lower[non_iterable] <- start[non_iterable]
      multstart_upper[non_iterable] <- start[non_iterable]
    }

    # Check order
    multstart_lower <- multstart_lower[names(start)]
    multstart_upper <- multstart_upper[names(start)]


    # Turn into named vectors again

    multstart_lower <- as.numeric(multstart_lower)
    multstart_upper <- as.numeric(multstart_upper)

    names(multstart_lower) <- names(start)
    names(multstart_upper) <- names(start)
  } else {
    stop("multstart_iter should be of length 1 or of the same length as the number of parameters")
  }

  out <- list(
    multstart_lower = multstart_lower,
    multstart_upper = multstart_upper
  )

  return(out)
}

#' Estimate the standard error of parameters using the delta method
#'
#' This function returns the SE of parameters, as a proportion of their
#' parameter estimates, calculated using the delta method as implemented in the
#' car package. If the parameter does not exist, this function will return NA.
#'
#' @param fit A fit object for a particular model - in this case the actual lm
#'   or nls object
#' @param expression The expression for which the SE is being calculated
#'
#' @return The SE as a proportion of the parameter estimate.
#' @export
#'
#' @examples
#' data(mtcars)
#' a <- lm(mpg ~ cyl + disp, data=mtcars)
#'
#' get_se(a, "cyl")
#' get_se(a, "cyl/disp")
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
get_se <- function(fit, expression) {

  secov <- function(fit, expression) {
    carout <- car::deltaMethod(fit, expression)
    abs( carout$SE / carout$Estimate )
  }

  possibly_se <- purrr::possibly(secov, otherwise = NA)

  possibly_se(fit, expression)

}

#' Convert between different units of radioactivity
#'
#' @param values Radioactivity values which should be converted. This can be a single value or a vector of values.
#' @param from_units The units from which the radioactivity should be converted.
#' @param to_units The units to which the radioactivity should be converted.
#'
#' @return The converted radioactivity values in the new units
#' @export
#'
#' @examples
#' unit_convert(1, "nCi", "Bq")
#' unit_convert(1:5, "nCi", "Bq")
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
unit_convert <- function(values,
                         from_units = c("GBq", "MBq", "kBq","Bq",
                                       "mBq", "uBq", "nBq","pBq",
                                       "GCi", "MCi", "kCi","Ci",
                                       "mCi", "uCi", "nCi", "pCi"),
                         to_units = c("GBq", "MBq", "kBq","Bq",
                                     "mBq", "uBq", "nBq","pBq",
                                     "GCi", "MCi", "kCi","Ci",
                                     "mCi", "uCi", "nCi", "pCi")) {

  if(length(from_units) > 1 | length(to_units) > 1) {
    stop("Units must be selected to convert between")
  }

  # Match the units

  from <- match.arg(from_units, c("GBq", "MBq", "kBq", "Bq",
                                    "mBq", "uBq", "nBq", "pBq",
                                    "GCi", "MCi", "kCi", "Ci",
                                    "mCi", "uCi", "nCi", "pCi"))

  to <- match.arg(to_units, c("GBq", "MBq", "kBq","Bq",
                                  "mBq", "uBq", "nBq","pBq",
                                  "GCi", "MCi", "kCi","Ci",
                                  "mCi", "uCi", "nCi", "pCi"))

  # Extract the meaning from the units

  from_bqci <- ifelse( grepl("Bq", from), "Bq", "Ci")
  from_mult <- strsplit(x = from, split = from_bqci)[[1]]

  to_bqci <- ifelse( grepl("Bq", to), "Bq", "Ci")
  to_mult <- strsplit(x = to, split = to_bqci)[[1]]

  from_multno <- dplyr::case_when(
    from_mult == "G" ~ 1e9,
    from_mult == "M" ~ 1e6,
    from_mult == "k" ~ 1e3,
    from_mult == ""  ~ 1,
    from_mult == "m" ~ 1e-3,
    from_mult == "u" ~ 1e-6,
    from_mult == "n" ~ 1e-9,
    from_mult == "n" ~ 1e-12,
    TRUE ~ 1
  )

  to_multno <- dplyr::case_when(
    to_mult == "G" ~ 1e9,
    to_mult == "M" ~ 1e6,
    to_mult == "k" ~ 1e3,
    to_mult == ""  ~ 1,
    to_mult == "m" ~ 1e-3,
    to_mult == "u" ~ 1e-6,
    to_mult == "n" ~ 1e-9,
    to_mult == "p" ~ 1e-12,
  )

  # Removing the from multiplier

  values <- values * from_multno

  # Converting between Bq and Ci

  values <- dplyr::case_when(
    from_bqci == to_bqci ~ values,
    from_bqci == "Ci" & to_bqci == "Bq" ~ values * 37000000000,
    from_bqci == "Bq" & to_bqci == "Ci" ~ values / 37000000000,
  )

  # Adding the to multiplier

  values <- values / to_multno

  # Return

  return(values)

}

#' Reverse Decay Correction of TAC values
#'
#' This function uses the equation described here to perform the calculations: http://www.turkupetcentre.net/petanalysis/decay.html
#'
#' @param t_start The starting times of the frames in minutes.
#' @param t_end The end times of the frames in minutes.
#' @param tac The uncorrected radioactivity values.
#' @param radioisotope The radioisotope.
#'
#' @return The tac before decay correction
#' @export
#'
#' @examples
#' data(pbr28)
#' s1 <- pbr28$tacs[[1]]
#'
#' # Assuming the data were not decay-corrected (they are)
#' decay_uncorrect(
#'   s1$StartTime/60,
#'  (s1$StartTime + s1$Duration)/60,
#'  tac = s1$FC)
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
decay_uncorrect <- function(t_start, t_end, tac, radioisotope = c("C11", "O15", "F18")) {

  radioisotope <- match.arg(radioisotope, c("C11", "O15", "F18"))

  hl <- dplyr::case_when(
    radioisotope == "C11" ~ 20.4,
    radioisotope == "O15" ~ 2.05,
    radioisotope == "F18" ~ 109.8
  )

  lambda <- log(2) / hl

  term1 <- exp(lambda*t_start) * lambda * (t_end - t_start)
  term2 <- 1 - exp( -1 * lambda * (t_end - t_start) )

  tac_uncor <- tac * ( term2 / term1 )

  tac_uncor[is.nan(tac_uncor)] <- 0

  return(tac_uncor)

}

#' Perform Decay Correction of TAC values
#'
#' This function uses the equation described here to perform the calculations: http://www.turkupetcentre.net/petanalysis/decay.html
#'
#' @param t_start The starting times of the frames in minutes.
#' @param t_end The end times of the frames in minutes.
#' @param tac_uncor The uncorrected radioactivity values.
#' @param radioisotope The radioisotope.
#'
#' @return The tac after decay correction
#' @export
#'
#' @examples
#' data(pbr28)
#' s1 <- pbr28$tacs[[1]]
#'
#' # Assuming the data were not decay-corrected (they are)
#' decay_correct(
#'   s1$StartTime/60,
#'  (s1$StartTime + s1$Duration)/60,
#'  tac_uncor = s1$FC)
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
decay_correct <- function(t_start, t_end, tac_uncor,
                          radioisotope = c("C11", "O15", "F18")) {

  radioisotope <- match.arg(radioisotope, c("C11", "O15", "F18"))

  hl <- dplyr::case_when(
    radioisotope == "C11" ~ 20.4,
    radioisotope == "O15" ~ 2.05,
    radioisotope == "F18" ~ 109.8
  )

  lambda <- log(2) / hl

  term1 <- exp(lambda*t_start) * lambda * (t_end - t_start)
  term2 <- 1 - exp( -1 * lambda * (t_end - t_start) )

  tac <- tac_uncor * ( term1 / term2 )

  tac[is.nan(tac)] <- 0

  return(tac)

}

#' Cumulative integral over frames
#'
#' This function uses the frame radioactivity values and durations to create a
#' cumulative integral. It sums the multiple of radioactivity and duration
#' for all frames preceding the current frame, and halves that of the current
#' frame. This is an alternative to using the trapezoidal method.
#'
#' @param dur Frame durations in minutes.
#' @param tac Radioactivity concentrations for each frame.
#'
#' @return The cumulative integral of radioactivity.
#' @export
#'
#' @examples
#' data(pbr28)
#'
#' tac <- pbr28$tacs[[2]]$WB
#' dur <- pbr28$tacs[[2]]$Duration / 60
#'
#' frame_integ <- frame_cumsum(dur, tac)
#'
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
frame_cumsum <- function(dur, tac) {

  if( max(dur) > 60 ) {
    warning("The frame durations are very long. Make sure they are in minutes and not seconds")
  }

  durtac <- dur*tac
  cumdurtac <- cumsum(durtac)

  before_cumdurtac <- c(0, cumdurtac[-(length(cumdurtac))])
  halfway_durtac <- 0.5 * durtac

  out <- before_cumdurtac + halfway_durtac

  return(out)

}

camel <- function(x){ #function for camel case
  capit <- function(x) paste0(toupper(substring(x, 1, 1)), substring(x, 2, nchar(x)))
  first <- strsplit(x, "\\_")[[1]][1]
  rest <- sapply(strsplit(x, "\\_")[[1]][-1],
                 function(x) paste(capit(x), sep=""))

  paste(first, rest, collapse="")
}

get_units_radioactivity <- function(x) {

  rad <- gsub("[\\.\\/].*", "", x)
  vol <- gsub("\\w*[\\.\\/]", "", x)
  vol <- gsub("\\-1$", "", vol)

  list(rad=rad, vol=vol)

}


#' Predict Blood Model Values: Generic Function
#'
#' Function to predict values from blood models
#'
#' @param modelout The output object of the model fitting procedure.
#' @param ... Additional optional arguments.
#'
#' @return The predicted values
#'
#' @details This function uses the \code{class(model)} name to call the correct function to plot the model fit.
#'
#' @examples
#'
#' blooddata <- pbr28$blooddata[[1]]
#' blooddata <- bd_blood_dispcor(blooddata)
#' blood <- bd_extract(blooddata, output = "Blood")
#' blood_fit <- blmod_splines(blood$time,
#'                            blood$activity,
#'                            Method = blood$Method)
#'
#' predict_blmod(blood_fit)
#'
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
#'
#' @export
predict_blmod <- function(modelout, ...) {
  modelname <- class(modelout)[1]
  input_list <- as.list(substitute(list(...)))

  out <- do.call(what = paste0("predict_", modelname), args = list(modelout, ...))

  as.numeric(out)

}



#' Extract coefficients
#'
#' Function to extract the coefficients from kinfit objects
#'
#' @param x The output object of the model fitting procedure.
#' @param ... Additional optional arguments.
#'
#' @examples
#' \dontrun{
#' loganout <- Loganplot(t_tac, tac, input, 10, weights)
#' coef(loganout)
#' }
#'
#' @author Granville J Matheson, \email{mathesong@@gmail.com}
coef_kinfit <- function(x, ...) {
  out <- x$par
  return(out)
}

interpends <- function(x, y, xi, method = "linear", yzero = NULL) {
  if (is.null(yzero)) { # start point
    yzero <- head(y, 1)
  }

  if (head(xi, 1) < head(x, 1)) {
    x <- c(head(xi, 1), x)
    y <- c(yzero, y)
  }

  if (tail(xi, 1) > tail(x, 1)) {
    x <- c(x, tail(xi, 1))
    y <- c(y, tail(y, 1))
  }

  pracma::interp1(x, y, xi, method)
}