R/calculations.R

In detzrcr: Compare Detrital Zircon Suites

#' Check concordancy of input ages
#'
#' Check the concordancy of the U-Pb data and return the data within the desired
#' discordancy limit.
#'
#' @param dat data.frame containing at least ages and percentage of discordancy
#' @param disc_lim Discordancy limit
#' @return Concordant data
#' @export
check_conc <- function(dat, disc_lim = 10) {
  if ('disc' %in% names(dat)) {
    llim <- -disc_lim
    ulim <- disc_lim
    conc <- dat[(dat$disc >= llim & dat$disc <= ulim), ]
  } else {
    stop('Requires column with name disc')
  }
}

#' Calculate 1D density of age data
#'
#' Calculate the 1d density of U-Pb age data using KDE or PDD.
#'
#' @param dat data.frame containing at least ages and percentage of discordancy
#' @param bw Bandwidth
#' @param type Type to calculate 'kde': proper KDE; 'pdd': detrital zircon PDD
#' @param age_range Range over which to calculate density
#' @return Density
#' @export
calc_dens <- function(dat, bw=30, type='kde', age_range=c(0, 4560)) {
  n <- age_range[2] - age_range[1]
  if (type == 'kde') {
    dens <- stats::density(dat$age, bw=bw, n=n, from=age_range[1],
                           to=age_range[2])
    x <- dens$x
    y <- dens$y
   }
  if (type == 'pdd') {
    if ('uncert' %in%  names(dat)) {
      x <- seq(age_range[1], age_range[2])
      y <- sapply(x, FUN=function(x) sum(stats::dnorm(x, dat$age, dat$uncert)))
      y <- y / length(x)
    } else {
      stop('Cannot calculate pdd without individual uncertainty each grain')
    }
  }
  dens <- data.frame(x=x, y=y)
  if(!is.null(dat$sample)) {
    dens$sample <- rep(dat$sample[1], nrow(dens))
  }
  return(dens)
}

#' Calculate scaled 1d density
#'
#' Calculates 1d density of age data and scales it so that it can be plotted
#' in the same plot of a histogram of the age data
#'
#' @param dat data.frame
#' @param binwidth Histogram binwidth
#' @param bw Density bandwidth
#' @param type 'kde': KDE; 'pdd': detrital zircon PDD
#' @param age_range Age range to calculated density over
#'
#' @return Returns density
#' @export
calc_dens_hist <- function(dat, binwidth=50, bw=30, type='kde',
                           age_range=c(0, 4560)) {
  dens <- calc_dens(dat, bw, type, age_range)
  max_y <- max(dens$y)
  # binmin <- find_plot_min(dat$age, accuracy=binwidth)
  # binmax <- find_plot_max(dat$age, accuracy=binwidth)
  # breaks <- seq(binmin, binmax, binwidth)
  hist_data <- graphics::hist(dat$age, breaks=seq(0, 4560, binwidth),
                              plot=FALSE, right=TRUE)
  dens$y <- dens$y * (max(hist_data$counts) / max_y)
  return(dens)
}

#' Find minimum value for plotting
#'
#' Find the minimum value for histogram plotting.
#'
#' @param x vector of values
#' @param accuracy round to nearest
#' @export
find_plot_min <- function(x, accuracy=100) {
  minimum <- floor(min(x) / accuracy) * accuracy
}

#' Find maximum value for plotting.
#'
#' Find the maximum value for histogram plotting.
#'
#' @param x vector of values
#' @param accuracy round to nearest
#' @export
find_plot_max <- function(x, accuracy=100) {
  maximum <- ceiling(max(x) / accuracy) * accuracy
}

#' Wrapper function for \code{find_plot_min} and \code{find_plot_max}
#'
#' Find the minimum and maximum values for histogram plotting.
#'
#' @param x Age data
#' @param accuracy Round to nearest
#'
#' @return Returns vector of minimum and maximum plotting values
#' @export
find_plot_min_max <- function(x, accuracy=100) {
  return(c(find_plot_min(x, accuracy), find_plot_max(x, accuracy)))
}

#' Calculate 1d likeness of detrital zircon populations
#'
#' Calculates the likeness of detrital zircon populations in 1 dimension after
#' Satkoski et al. (2013).
#'
#' @param x vector
#' @param y vector
#' @param bw bandwidth
#' @param digits number, round result to significant digits
#'
#' @references Satkoski, A.M., Wilkinson, B.H., Hietpas, J., Samson, S.D., 2013.
#' Likeness among detrital zircon populations - An approach to the comparison of
#' age frequency data in time and space. GSA Bulletin 125, 1783-1799.
#' @export
satkoski_1d <- function(x, y, bw=30, digits=3) {

  a <- calc_dens(x, bw=bw)
  b <- calc_dens(y, bw=bw)

  L <- 1 - (sum(abs(a$y - b$y)) / 2)
  round(L, digits)
}

#' Calculate 2d (age and Lu-Hf) likeness of detrital zircon populations
#'
#' Calculates the likeness of detrital zircon populations in 2 dimensions after
#' Satoski et al. (2013).
#'
#' @param x vector
#' @param y vector
#' @param bw vector of density bandwidths
#' @param digits number, round result to significant digits
#'
#' @references Satkoski, A.M., Wilkinson, B.H., Hietpas, J., Samson, S.D., 2013.
#' Likeness among detrital zircon populations - An approach to the comparison of
#' age frequency data in time and space. GSA Bulletin 125, 1783-1799.
#' @export
satkoski_2d <- function(x, y, bw=c(30, 2.5), digits=3) {
  if (all(is.na(x$ehf_i)) | all(is.na(y$ehf_i))) {
    return(NA)
  } else {
    n <- 100
    bw <- bw * 4
    lims <- c(0, 4560, -30, 30)
    x <- x[!is.na(x$ehf_i), ]
    y <- y[!is.na(y$ehf_i), ]
    a <- MASS::kde2d(x=x$age, y=x$ehf_i, h=bw, n=n, lims=lims)$z
    b <- MASS::kde2d(x=y$age, y=y$ehf_i, h=bw, n=n, lims=lims)$z
    a <- a / sum(a)
    b <- b / sum(b)
    L <- 1 - (sum(abs(a - b)) / 2)
    round(L, digits)
  }
}

#' Pairwise Satkoski likeness
#'
#' Populate a matrix with pairwise Satkoski 1d likeness.
#'
#'@param dat data.frame
#'@param bw density bandwidth
#'@param digits number, round result to significant digits
#' @references Satkoski, A.M., Wilkinson, B.H., Hietpas, J., Samson, S.D., 2013.
#' Likeness among detrital zircon populations - An approach to the comparison of
#' age frequency data in time and space. GSA Bulletin 125, 1783-1799.
#' @export
satkoski_1d_matrix <- function(dat, bw=30, digits=3) {
  populate_matrix(dat, FUN=satkoski_1d, bw=bw, digits=3)
}

#' Pairwise 2d Satkoski likeness
#'
#' Populate a matrix with pairwise Satkoski 12 likeness.
#'
#'@param dat data.frame
#'@param bw vector of density bandwidths
#'@param digits number, round result to significant digits
#'
#' @references Satkoski, A.M., Wilkinson, B.H., Hietpas, J., Samson, S.D., 2013.
#' Likeness among detrital zircon populations - An approach to the comparison of
#' age frequency data in time and space. GSA Bulletin 125, 1783-1799.
#' @export
satkoski_2d_matrix <- function(dat, bw=c(30, 2.5), digits=3) {
  populate_matrix(dat, FUN=satkoski_2d, bw=bw, digits=digits)
}

#' Populate matrix
#'
#' @param dat data.frame
#' @param FUN Function used to populate matrix
#' @param ... Additional parameters passed to function
#'
#' @return Populated matrix
#' @export
populate_matrix <- function(dat, FUN, ...) {
  n <- length(unique(dat$sample))
  if (n < 2) stop('Select more samples')
  name <- as.character(unique(dat$sample))
  len <- length(name)
  mat <- matrix(nrow=n, ncol=n)
  colnames(mat) <- name
  rownames(mat) <- name
  for (i in 1:(len - 1)) {
    for (j in 2:len) {
      if (!(i == j) & ((is.na(mat[i, j])) | is.na(mat[j, i]))) {
        L <- FUN(dat[dat$sample == name[i], ],
                         dat[dat$sample == name[j], ], ...)
        mat[i, j] <- L
        mat[j, i] <- L
      }
    }
  }
  return(mat)
}

#' Calculate hafnium values.
#'
#' Calculates the initial 176Hf/177Hf values, the initial epsilon hafnium
#' values, the model age using the measured 176Lu/177Hf value and
#' the model age assuming the parental magma was produced from an average continental crust
#' (176Lu/177Hf = 0.015) that originally was derived from the depleted mantle
#' (Griffin, 2004).
#'
#' @param dat data.frame, list or matrix of hafnium values
#' @param constants vector of constants which must be in the order
#' decay constant 176Lu, 176/177Hf CHUR, 176Lu/177Hf CHUR, 176/177Hf DM,
#' 176Lu/177Hf DM and 176Lu/177Hf value used for two-stage depleted mantle
#' model age calculations
#' @export
#' @references Bouvier, A., Vervoort, J.D. & Patchett, P.J. 2008.
#' The Lu-Hf and Sm-Nd isotopic composition of CHUR:
#' Constraints from unequilibrated chondrites and implications for the bulk composition of terrestrial planets.
#' Earth And Planetary Science Letters 273(1-2), 48-57.
#' @references Griffin, W., Belousova, E., Shee, S., Pearson, N. and O'Reilly, S. 2004.
#' Archean crustal evolution in the northern Yilgam Craton:
#' U-Pb and Hf-isotope evidence from detrital zircons. Precambrian Research, 231-282.
#' @references Soderlund, U., Patchett, J., Vervoort, J. & Isachsen, C. 2004.
#' The Lu-176 decay constant determined by Lu-Hf and U-Pb isotope systematics of Precambrian mafic intrusions.
#' Earth And Planetary Science Letters 219(3-4), 311-324.
calc_hf <- function(dat, constants) {
  lambda_lu <- constants[1]
  hfhf_chur <- constants[2]
  luhf_chur <- constants[3]
  hfhf_dm <- constants[4]
  luhf_dm <- constants[5]
  luhf_zrc <- constants[6]
  if ('ehf_i' %in% names(dat)) {
    dat <- dat[, !(names(dat) %in% c('hfhf', 'hfhfse', 'luhf'))]
  }
  if (is.vector(dat) || ncol(dat) < 3) {
    stop('Data must contain at least three columns')
  }
  hf_chur <- hfhf_chur - (luhf_chur * (exp(lambda_lu * (dat$age * (10^6))) - 1))
  if ('ehf_i' %in% names(dat)) {
    hf_i <- hf_chur * ((dat$ehf_i / 10^4) + 1)
  } else {
    if ('hfhf' %in% names(dat)) {
      hf_i <- dat$hfhf - (dat$luhf * (exp(lambda_lu * (dat$age * (10^6))) - 1))
      ehf_i <- ((hf_i / hf_chur) - 1) * 10^4
    }
  }
  if ('hfhfse' %in% names(dat)) {
    hf_2se <- 2 * dat$hfhfse
    ehf_2se <- ((((hf_i + hf_2se) / (hfhf_chur - (luhf_chur * (exp(lambda_lu *
                                                 (dat$age * (10^6))) - 1)))) -
                                                  1)*10^4) - ehf_i
  } else {
    hf_2se <- rep(NA, nrow(dat))
    ehf_2se <- hf_2se
  }
  hf_dm <- hf_i + (luhf_zrc * (exp(lambda_lu * (dat$age * (10^6))) -1))
  t_dm2 <- ((1 / (lambda_lu)) *
              log(((hf_dm - hfhf_dm) / (luhf_zrc - luhf_dm)) + 1)) / 10^6
  if ('ehf_i' %in% names(dat)) {
    dat <- cbind(dat, ehf_2se, hf_i, hf_2se, hf_chur, t_dm2)
  } else {
    if ('hfhf' %in% names(dat)) {
      t_dm <- (((1 / (lambda_lu)) *
                  log(((dat$hfhf - hfhf_dm) / (dat$luhf - luhf_dm)) + 1)) / 10^6)
      dat <- cbind(dat, hf_i, hf_2se, ehf_i, ehf_2se, t_dm, t_dm2)
    }
  }
}

#' Dvoretzky-Kiefer-Wolfowitz inequality
#'
#' Calculate confidence bands for ecdfs using the Dvoretzky-Kiefer-Wolfowitz
#' inequality.
#'
#' @param dat data.frame
#' @param column which column to use
#' @param alpha Desired alpha level
#'
#' @return data.frame with ecdf and confidence bands
#' @export
#'
#' @references Dvoretzky, A., Kiefer, J., Wolfowitz, J., 1956. Asymptotic
#' Minimax Character of the Sample Distribution Function and of the Classical
#' Multinomial Estimator. Ann. Math. Stat. 27, 642-669.
#' doi:10.1214/aoms/1177728174
calc_dkw <- function(dat, column='age', alpha=0.05) {
  if ('sample' %in% names(dat)) {
    sample_name <- dat$sample[1]
  } else {
    sample_name <- NA
  }
  x <- dat[, column]
  x <- x[!is.na(x)]
  y <- stats::ecdf(x)(sort(x))
  n <- length(x)
  x <- c(0, sort(x), 4560)
  y <- c(0, y, 1)
  epsilon <- sqrt(log(2 / alpha) / (2 * n))
  low <- pmax(y - epsilon, 0)
  high <- pmin(y + epsilon, 1)
  sample <- rep(sample_name, length(x))
  data.frame(x=x, y=y, low=low, high=high, sample=sample)
}

#' Calculate 1-O
#'
#' @param dat1 data.frame
#' @param dat2 data.frame
#' @param column string of name of column to use ('age' or 't_dm2')
#' @param alpha alpha level
#' @param digits number of digits
#'
#' @return 1-O
#' @export
#' @references Andersen, T., Elburg, M., Cawthorn-Blazeby, A., 2015.
#' U-Pb and Lu-Hf zircon data in young sediments reflect sedimentary recycling
#' in eastern South Africa. J. Geol. Soc. London. 2006-2015.
#' doi:10.1144/jgs2015-006
#'
calc_o_param <- function(dat1, dat2, column, alpha=0.05, digits=2) {
  x <- dat1[, column]
  y <- dat2[, column]
  if (all(is.na(x)) | all(is.na(y))) {
    return(NA)
  }
  else {
    epsilon_one <- sqrt(log(2 / alpha) / (2 * length(x)))
    epsilon_two <- sqrt(log(2 / alpha) / (2 * length(y)))
    x_sort <- sort(x)
    x_y <- stats::ecdf(x_sort)(x_sort)
    y_sort <- sort(y)
    y_y <- stats::ecdf(y_sort)(y_sort)
    x_sort <- c(0, x_sort, 4560)
    x_y <- c(0, x_y, 1)
    y_sort <- c(0, y_sort, 4560)
    y_y <- c(0, y_y, 1)
    x_out <- seq(0, 4560)
    interpolated_one <- stats::approx(x_sort, x_y, xout=x_out, ties='ordered')$y
    interpolated_two <- stats::approx(y_sort, y_y, xout=x_out, ties='ordered')$y
    interpolated_one_low <- pmax(interpolated_one - epsilon_one, 0)
    interpolated_one_up <- pmin(interpolated_one + epsilon_one, 1)
    interpolated_two_low <- pmax(interpolated_two - epsilon_two, 0)
    interpolated_two_up <- pmin(interpolated_two + epsilon_two, 1)
    up <- ifelse((((interpolated_two_up <= interpolated_one_up) &
                     (interpolated_two_up >= interpolated_one_low)) |
                    ((interpolated_one_up <= interpolated_two_up) &
                       (interpolated_one_up >= interpolated_two_low))), 1, 0)
    low <- ifelse((((interpolated_two_low >= interpolated_one_low) &
                      (interpolated_two_low <= interpolated_one_up)) |
                     ((interpolated_one_low >= interpolated_two_low) &
                        (interpolated_one_low <= interpolated_two_up))), 1, 0)
    O <- (length(up[up == 1]) + length(low[low == 1])) / (2 * length(x_out))
    round(1 - O, digits)
  }
  # delta <- epsilon_one + epsilon_two
  # absolute <- pmax((abs(interpolated_one - interpolated_two) - delta), 0)
  # absolute <- absolute[absolute != 0]
  # round(1- (1 - length(absolute) / length(x_out)), digits)
}

#' Populate matrix with age 1-O
#'
#' @param dat data.frame
#' @param alpha alpha level
#' @param digits number of digits
#'
#' @return matrix of 1-O for ages
#' @export
#' @references Andersen, T., Elburg, M., Cawthorn-Blazeby, A., 2015.
#' U-Pb and LuHf zircon data in young sediments reflect sedimentary recycling
#' in eastern South Africa. J. Geol. Soc. London. 2006-2015.
#' doi:10.1144/jgs2015-006
o_param_matrix_age <- function(dat, alpha=0.05, digits=2) {
  populate_matrix(dat, FUN=calc_o_param, column='age', alpha=alpha,
                  digits=digits)
}

#' Populate matrix with model age 1-O
#'
#' @param dat data.frame
#' @param alpha alpha level
#' @param digits number of digits
#'
#' @return matrix of 1-O for model ages
#' @export
#' @references Andersen, T., Elburg, M., Cawthorn-Blazeby, A., 2015.
#' U-Pb and Lu-Hf zircon data in young sediments reflect sedimentary recycling
#' in eastern South Africa. J. Geol. Soc. London. 2006-2015.
#' doi:10.1144/jgs2015-006
o_param_matrix_tdm <- function(dat, alpha=0.05, digits=2) {
  populate_matrix(dat, FUN=calc_o_param, column='t_dm2', alpha=alpha,
                  digits=digits)
}

#' Combine two square matrices
#'
#' @param mat1 Matrix for upper triangle
#' @param mat2 Matrix for lower triangle
#'
#' @export
#'
combine_matrices <- function(mat1, mat2) {
  mat1[lower.tri(mat1)] <- NA
  mat1[is.na(mat1)] <- mat2[is.na(mat1)]
  mat1
}

#' Produce CHUR and DM lines
#'
#' Calculate CHUR and DM lines used for epsilon-Hf vs. age and 176/177Hf vs. age
#' plots.
#'
#' @param range range over which to calculate lines
#' @param plot_type 'ehf' = epsilon-Hf; any thing else gives 176/177Hf
#' @param constants vector of constants which must be in the order
#' decay constant 176Lu, 176/177Hf CHUR, 176Lu/177Hf CHUR, 176/177Hf DM and
#' 176Lu/177Hf DM
#'
#' @export
#'
#' @references Griffin, W., Pearson, N., Belousova, E., Jackson, S., van
#' Achterbergh, E., O'Reilly, S. and Shee, S. 2000. The Hf isotope composition
#' of cratonic mantle:
#' LAM-MC-ICPMS analysis of zircon megacrysts in kimberlites.
#' Geochimica et Cosmochimica Acta 64(1), 133-147.
#' @references Soderlund, U., Jonathan Patchett, P., Vervoort, J.D. and Isachsen
#' C.E. 2004. The 176Lu decay constant determined by Lu-Hf and U-Pb isotope
#' systematics of Precambrian mafic intrusions. Earth and Planetary Science
#' Letters 219, 311-324.
#' @references Bouvier, A., Vervoort, J.D. and Jonathan Patchett P. 2008. The
#' Lu-Hf and Sm-Nd isotopic composition of CHUR: Constraints from unequilibrated
#' chondrites and implications for the bulk composition of terrestrial planets.
#' Earth and Planetary Science Letters 273, 48-57.
hf_lines <- function(range=c(0, 4560), plot_type='ehf', constants) {
  lambda_lu <- constants[1]
  hfhf_chur <- constants[2]
  luhf_chur <- constants[3]
  hfhf_dm <- constants[4]
  luhf_dm <- constants[5]
  x <- seq(range[1], range[2], by=1)
  length_x <- length(x)
  #CHUR (Bouvier,2008)
  chur <- hfhf_chur - luhf_chur * ((exp(lambda_lu * x * 10^6)) - 1)
  #DM (Griffin,2000)
  dm <- hfhf_dm - luhf_dm *((exp(lambda_lu * x * 10^6)) - 1)
  if (plot_type == 'ehf') {
    dm <- ((dm / chur) - 1) * 10^4
    chur <- rep(0, length_x)
  }
  type_chur <- rep('CHUR', length_x)
  type_dm <- rep('DM', length_x)
  type <- c(type_chur, type_dm)
  x <- rep(x, 2)
  line <- c(chur, dm)
  linedata <- data.frame(x=x, line=line, type=type)
}


#' Calculate quantiles
#'
#' Split up data.frame by sample-column and calculate quantiles
#'
#' @param dat data.frame
#' @param column which column in data.frame to use
#' @param alpha alpha-level (not yet used)
#' @param type type of quantile calculation (passed on to stats::quantile)
#'
#' @export
#'
calc_quantiles <- function(dat, column='t_dm2', alpha=0.05, type=8) {
  x <- dat[, column]
  quantiles <- stats::aggregate(x=x, by=list(dat$sample),
                                FUN=stats::quantile, type=type, na.rm=TRUE)
  x <- as.data.frame(quantiles$x)
  x$iqr <- x$`75%` - x$`25%`
  x$sample <- quantiles$Group.1
  names(x) <- c('zero', 'twentyfive', 'fifty', 'seventyfive', 'onehundred',
                'iqr', 'sample')
  return(x)
}

#' Calculate confidence bands for lower and upper quartile
#'
#' @param dat data.frame
#' @param column column to use for calculations
#' @param alpha alpha-level
#'
#' @export
#'
quant_bounds <- function(dat, column='t_dm2', alpha=0.05) {
  column <- dat[, column]
  sample <- dat$sample[1]
  if (all(is.na(column))) {
    data.frame(x=NA, y=NA, ymin=NA, ymax=NA, xmin=NA, xmax=NA, sample=sample)
  } else {
    sort_column <- sort(column)
    y <- stats::ecdf(column)(sort(column))
    epsilon <- sqrt(log(2 / alpha) / (2 * length(column)))
    low <- pmax(y - epsilon, 0)
    high <- pmin(y + epsilon, 1)
    sort_age_low <- data.frame(x=sort(column), y=low)
    sort_age_high <- data.frame(x=sort(column), y=high)
    ll <- stats::approx(x=sort_age_low$y, y=sort_age_low$x, xout=c(0.25),
                        ties='ordered')$y
    lu <- stats::approx(x=sort_age_low$y, y=sort_age_low$x, xout=c(0.75),
                        ties='ordered')$y
    if(is.na(lu)) lu <- 4560
    ul <- stats::approx(x=sort_age_high$y, y=sort_age_high$x, xout=c(0.25),
                        ties='ordered')$y
    if(is.na(ul)) ul <- 0
    uu <- stats::approx(x=sort_age_high$y, y=sort_age_high$x, xout=c(0.75),
                        ties='ordered')$y
    lq_dist <- stats::quantile(column, probs=c(0.25), type=8, na.rm=TRUE)
    uq_dist <- stats::quantile(column, probs=c(0.75), type=8, na.rm=TRUE)
    data.frame(x=lq_dist, y=uq_dist, ymin=uu, ymax=lu, xmin=ul,  xmax=ll, sample)
  }
}

#' Calculate mixing model
#'
#' Gaussian mixing model for detrital zircon data, using lower quantile upper
#' quantile plot
#'
#' @param mu1 first mean
#' @param sig1 first standard deviation
#' @param mu2 second mean
#' @param sig2 second standard deviation
#'
#' @examples
#' dzr_mix(500, 50, 1000, 100)
#' @export
#'
dzr_mix <- function(mu1, sig1, mu2, sig2) {
  X <- rep(NA, 100)
  S <- rep(NA, 100)
  lq <- matrix(nrow=11, ncol=100)
  uq <- lq
  for (j in seq(1, 100, 1)) {
    for (i in seq(1, 11, 1)) {
      X[1: (100 - 10 * (i - 1))] <- mu1
      X[(101 - 10 * (i - 1)): 100] <- mu2
      S[1: (100 - 10 * (i - 1))] <- sig1
      S[(101 - 10 * (i - 1)): 100] <- sig2

      x <- stats::rnorm(100, mean=X, sd=S)
      y <- stats::ecdf(x)(x)
      lq[i, j] <- stats::quantile(x, probs=0.25, type=8)
      uq[i, j] <- stats::quantile(x, probs=0.75, type=8)
    }
  }
  data.frame(lq=rowMeans(lq), uq=rowMeans(uq))
}

#' Ready 1-O matrix for tile plot
#'
#' @param x 1-O parameter vector
#'
#' @export
tile_func <- function(x) {
  if (is.na(x)) {
    return(NA)
  }
  if (x >= 0.05) {
    return(2)
  }
  if (x > 0 & x < 0.05) {
    return(1)
  } else {
    return(x)
  }
}

#' Apply tile_func to vector
#'
#' @param z 1-O parameter vector
#'
#' @export
#'
tiling <- function(z) {
  sapply(z, tile_func)
}

#' Produce data.frame of 1-O matrix suitable for geom_tile
#'
#' @param dat data.frame
#' @param type What to calculate
#'
#' @export
#'
make_tiling <- function(dat, type) {
  if (type == 'age') {
    mat <- o_param_matrix_age(dat$age)
    mat <- mat[, rev(seq_len(ncol(mat)))]
    tile_mat <- data.frame(x=rownames(mat)[row(mat)], y=colnames(mat)[col(mat)],
                           z=as.factor(tiling(c(mat))))
    tile_mat$x <- factor(tile_mat$x, levels=rownames(mat))
    tile_mat$y <- factor(tile_mat$y, levels=colnames(mat))
  }
  if (type == 'tdm') {
    mat <- o_param_matrix_tdm(dat$hf)
    mat <- mat[, rev(seq_len(ncol(mat)))]
    tile_mat <- data.frame(x=rownames(mat)[row(mat)], y=colnames(mat)[col(mat)],
                           z=as.factor(tiling(c(mat))))
    tile_mat$x <- factor(tile_mat$x, levels=rownames(mat))
    tile_mat$y <- factor(tile_mat$y, levels=colnames(mat))
  }
  if (type == 'combine') {
    mat_age <- o_param_matrix_age(dat$age)
    mat_tdm <- o_param_matrix_tdm(dat$hf)
    mat_age[lower.tri(mat_age)] <- NA
    mat <- mat_age
    mat[is.na(mat)] <- mat_tdm[is.na(mat)]
    mat <- t(mat)
    mat <- mat[, rev(seq_len(ncol(mat)))]
    tile_mat <- data.frame(x=rownames(mat)[row(mat)], y=colnames(mat)[col(mat)],
                           z=as.factor(tiling(c(mat))))
    tile_mat$x <- factor(tile_mat$x, levels=rownames(mat))
    tile_mat$y <- factor(tile_mat$y, levels=colnames(mat))
  }
  tile_mat
}

#' Calculate upper and lower concordia intercepts from discordant detrital
#' zircon data
#'
#' @param dat data.frame
#' @param step Chord spacing
#'
#' @export
#'
#' @references Reimink, J.R., Davies, J.H.F.L., Waldron, J.W.F., Rojas, X.
#' (2016). Dealing with discordance: a novel approach for analysing U-Pb
#' detrital zircon datasets. Journal of the Geological Society.
#' doi: 10.1144/jgs2015-114
#'
reimink <- function(dat, step=5) {
  t1 <- seq(0, 4500 - step, step)
  t2 <- seq(step, 4500, step)
  grid <- expand.grid(t1, t2)
  grid <- grid[grid$Var2 > grid$Var1,]
  dat$sigma75 <- stats::median(dat$sigma75)
  dat$sigma68 <- stats::median(dat$sigma68)
  result <- calc_p_apply(dat, grid$Var2, grid$Var1)
  lower <- stats::aggregate(result$p_disc, by=list(result$t1), max)
  upper <- stats::aggregate(result$p_disc, by=list(result$t2), max)
  names(lower) <- c('x', 'y')
  names(upper) <- c('x', 'y')
  lower$x <- as.numeric(as.character(lower$x))
  upper$x <- as.numeric(as.character(upper$x))
  lower$type <- rep('lower', nrow(lower))
  upper$type <- rep('upper', nrow(upper))
  out <- rbind(lower, upper)
}


#' Calculate intercepts and associated p-value
#'
#' @param dat data.frame
#' @param t2 upper intercept age
#' @param t1 lower intercept age
#'
#' @export
#'
calc_p_apply <- function(dat, t2, t1) {
  lambda238 <- 1.55125 * (10^-10)
  X <- dat$r75
  Y <- dat$r68
  sigma75 <- dat$sigma75
  sigma68 <- dat$sigma68
  rho <- dat$rho
  age <- dat$age * 10^6
  disc <- 1 - (Y / ((exp(lambda238 * age)-1)))
  ab <- calc_ab(t2, t1)
  slope <- ab$a
  yintercept <- ab$b
  n <- nrow(dat)
  p_disc <- rep(NA, length(slope))
  for (i in seq_len(length(slope))) {
    S <- atan2((2 * rho * sigma75 * sigma68) / (sigma75^2 - sigma68^2), 2)
    a <- ((slope[i] * cos(S)) - sin(S)) / (cos(S) + (slope[i] * sin(S)))
    b <- (yintercept[i] + (slope[i] * X) - Y) / (cos(S) + (slope[i] * sin(S)))
    B <- (sigma75 / sigma68) * a
    A <- b / sigma68
    p <- (1 / (2 * pi * sigma75 * sigma68)) * exp(-0.5 * (A^2 / (1 + B^2)))
    p <- p / n
    p_disc[i] <- sum(p * disc)
  }
  return(data.frame(t1, t2, p_disc))
}


#' Calculate slope and intercept
#'
#' @param t2 upper intercept
#' @param t1 lower intercept
#'
#' @export
#'
calc_ab <- function (t2, t1) {
  y2 <- concY(t2)
  y1 <- concY(t1)
  x2 <- concX(t2)
  x1 <- concX(t1)
  a = (y2 - y1) / (x2 - x1)
  b = y2 - a * x2
  ab <- data.frame(a=a, b=b)
}

#' Calculate U235 at given age
#'
#' @param age input age
#'
#' @export
#'
concX <- function(age) {
  age <- age * 10^6
  lambda235 <- 9.8485 * (10^-10)
  cx <- exp(lambda235*age) - 1
  return(cx)

}

#' Calculate U238 at given age
#'
#' @param age input age
#'
#' @export
#'
concY <- function(age) {
  age <- age * 10^6
  lambda238 <- 1.55125 * (10^-10)
  cy <- exp(lambda238*age) - 1
  return(cy)
}

#' Find maxima.
#'
#'
#' @param dist distribution.
#' @param xmin minimum value of distribution.
#' @param inc increment.
#'
#' @export
find_maxima <- function(dist, xmin, inc) {
  difference <- c(0, diff(dist))
  difference2 <- c(difference[2:length(difference)], 0)
  maximadata <- ((which(difference > 0 & difference2 <= 0) * inc) + xmin - 1)
}