rhoR coverage - 96.62%

Files
Source

###
#' @title Create Simulated codeSet
#' 
#' @description
#' Creates a simulated \code{\link{codeSet}} with the given parameters
#' 
#' @param length the length of the simulated \code{\link{codeSet}} to be created
#' @param baserate the \code{\link{baserate}} of the simulated \code{\link{codeSet}}
#' @param kappaMin the minimum kappa of the simulated \code{\link{codeSet}}
#' @param kappaMax the maximum kappa of the simulated \code{\link{codeSet}}
#' @param precisionMin the minimum precision of the simulated \code{\link{codeSet}}
#' @param precisionMax the maximum precision of the simulated \code{\link{codeSet}}
#' @param ... Parameters passed to createRandomSet (e.g. type = "set" or type = "ct")
#' @param tries the maximum number of tries to generate a valid set, smaller set lengths may require an increased number of tries
#' 
#' @details
#' \code{\link{codeSet}}s are generated by first picking a random kappa within its range and a random precision within its range.  If the random kappa, random precision, and baserate are not mathematically possible, then the precision is resampled from a range of mathematically possible values within its range.  A unique simulated \code{\link{codeSet}} is then constructed given these parameters.
#' 
#' @export
#' @return A \code{\link{codeSet}} that fulfills the given parameters
###
createSimulatedCodeSet = function(length, baserate, kappaMin, kappaMax, precisionMin, precisionMax, ..., tries = 50){
  if(length <= 0){stop("length must be positive")}
  if(baserate < 0){stop("baserate must be positive")}
  if(kappaMin < 0 | kappaMin > 1) stop("kappaMin must be between 0 and 1.")
  if(kappaMax < 0 | kappaMax > 1) stop("kappaMax must be between 0 and 1.")
  if(kappaMin > kappaMax){stop("kappaMin must be less than kappaMax")}
  if(precisionMin < 0 | precisionMin > 1) stop("precisionMin must be between 0 and 1.")
  if(precisionMax < 0 | precisionMax > 1) stop("precisionMax must be between 0 and 1.")
  if(precisionMin > precisionMax){stop("precisionMin must be less than precisionMax")}
  if(tries < 1){stop("tries must be greater than 1")}
  
  setFound = NULL
  validSet = F
  triesLeft = tries
  Ks = c()
  while(!validSet && triesLeft > 0) {
    set = createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax)
    k = kappa(set)
    
    if ((k > kappaMin && k < kappaMax) || any_equal(k, c(kappaMin, kappaMax))) {
      setFound = set
      validSet = T
    }
    
    Ks = c(Ks, k)
    triesLeft = triesLeft - 1
  }
  
  if (is.null(setFound)) {
    Ks_under_min <- Ks[which(Ks < kappaMin)]
    Kdist.to.min <- kappaMin - Ks_under_min
    Kdist.to.min.low <- ifelse(length(Ks_under_min) > 0, min(Kdist.to.min), NA)
    
    Ks_over_min  <- Ks[which(Ks > kappaMax)]
    Kdist.to.max <- Ks_over_min - kappaMax
    Kdist.to.max.low <- ifelse(length(Ks_over_min) > 0, min(Kdist.to.max), NA)
    
    Kdist.all_mins <- c(Kdist.to.min, Kdist.to.max)
    
    both_mins <- c(Kdist.to.min.low, Kdist.to.max.low)
    both_mins <- both_mins[!is.na(both_mins)]
    closest <- c(Ks_under_min, Ks_over_min) [which(Kdist.all_mins == min(both_mins))][1]
    
    stop(paste0(
      "Unable to create a set of the provided length (", length, ") ",
        "in the provided kappa range (", kappaMin, "\U2014", kappaMax, "). ",
        "The closest kappa obtained after ", tries, " tries was: ", 
        round(x = closest, digits = 5), "\n",
        "\tTry increasing the set length or the range of valid kappas"
    ))
  }
  
  return(setFound)
  
  # return(createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax, ...))
}

generateKPs = function(
  numNeeded, baserate,
  kappaMin, kappaMax,
  precisionMin, precisionMax,
  distributionType = 'FLAT',
  distributionLength = 100000
) {
  kappaDistribution = seq(kappaMin, kappaMax, length = distributionLength)

  if (distributionType == 'FLAT') {
    #probability for a flat distribution
    kappaProbability = NULL
  }
  else if (distributionType == 'BELL') {
    #probability for a bell curve
    kappaProbability = stats::dnorm(kappaDistribution, mean = 0.9, sd = 0.1)
  }

  precisionDistribution = seq(precisionMin, precisionMax, length = 10000)
  precisionProbability = NULL

  KPs = tryCatch({
    KPs = list()
    for (i in 1:numNeeded) {
      KPs[[length(KPs) + 1]] = genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate)
    }
    KPs
  },error=function(e){
    NULL
  })
  if(is.null(KPs)){stop("Could not generate a valid set given ranges of kappa and precision")}

  return (KPs);
}

###
#' @title Get Test Set
#' 
#' @description
#' This function gets a \emph{testSet} from a larger \code{\link{codeSet}} given certain sampling parameters.
#' 
#' @details
#' A \emph{testSet} is a \code{\link{codeSet}} that is a subset of a larger \code{\link{codeSet}} with a given set of properties.  A \emph{testSet} is constructed by sampling (without replacement) P rows from rows in the larger \code{\link{codeSet}} where the first rater's code was 1, and then appending an additional sample (without replacement) of R rows taken at random from the larger \code{\link{codeSet}} excluding rows included in the first P rows sampled. P is computed as the minbaserate * length of the \emph{testset}. R is computed as testSetLength - P. The result of this sampling procedure is to create a sample with a minimum baserate regardless of the baserate of the larger \code{\link{codeSet}}.If \emph{testSetBaserateInflation} is set to zero, the function selects rows at random.
#' 
#' @param set The \code{\link{codeSet}} from which the \emph{testSet} is taken
#' @param testSetLength The length of the \emph{testSet} to be taken
#' @param testSetBaserateInflation The minimum guaranteed \code{\link{baserate}} of the \emph{testSet}. Default to 0
#' 
#' 
#' @export
#' @return A \code{\link{codeSet}} with the properties specified
###

getTestSet = function(set, testSetLength, testSetBaserateInflation = 0){
  if(length(unlist(set)) < testSetLength){stop("testSetLength must be less than the length of set")}
  if(testSetLength %% 1 | testSetLength < 1) {stop("testSetLength value must be a positive integer.")}
  if(testSetBaserateInflation < 0){stop("testSetBaserateInfation must be positive")}
  if(testSetBaserateInflation >= 1){stop("testSetBaserateInflation must be below 1")}
  if(length(which(set[,1] == 1 & set[,2] == 1)) == nrow(set)){
    return(set[1:testSetLength,])
  }else if(length(which(set[,1] == 0 & set[,2] == 0)) == nrow(set)){
    if(testSetBaserateInflation == 0){
      return(set[1:testSetLength,])
    }else{
      stop("Not enough positives in first rater to inflate to this level")
    }
  }
  return(getHandSet(set = set, handSetLength = testSetLength, handSetBaserate = testSetBaserateInflation, returnSet = T))
}

###
# Precision Recall Set
#
# This creates a set with a given recall and precision
# @param precision This is the precision of the given set
# @param recall This is the recall of the given set
# @param length This is the length of the set
# @param baserate This is the baserate of the set
# @return returns a set with the given parameters
###
prset = function(precision, recall, length, baserate) {
  #generates a set given the different 4 quadrants
  gold1s = max(round(baserate * length , 0), 1); # TP + FN
  gold0s = length - gold1s; # FP + FN
  TP = max(round(gold1s * recall , 0) , 1); # TP
  FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
  gold = c(rep(1, gold1s), rep(0, gold0s));
  silver = c(rep(1, TP),rep(0, gold1s - TP), rep(1, FP), rep(0, gold0s - FP));

  return(cbind(gold, silver));
}

pr.ct = function(precision, recall, length, baserate) {
  #generates a set given the different 4 quadrants
  gold1s = max(round(baserate * length , 0), 1); # TP + FN
  gold0s = length - gold1s; # FP + FN
  TP = max(round(gold1s * recall , 0) , 1); # TP
  FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
  
  ct = matrix(0, ncol = 2, nrow = 2)
  ct[1,1] = TP
  ct[1,2] = gold1s - TP
  ct[2,1] = FP
  ct[2,2] = gold0s - FP
  return(ct);
}

any_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
  any(sapply(target, function(targ) abs(current - targ) < tolerance ))
}
all_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
  all(sapply(target, function(targ) abs(current - targ) < tolerance ))
}


#' Convert a codeset to a contingency table
#'
#' @param x codeset 
#'
#' @return contingency table as a 2x2 matrix
#' 
#' @export
as.contingency.table <- function(x) {
  y <- NULL
  
  if ( nrow(x) > 2 && ncol(x) == 2 ) {
    tp = length(which(x[,1] == 1 & x[,2] == 1))
    tn = length(which(x[,1] == 0 & x[,2] == 0))
    fp = length(which(x[,1] == 0 & x[,2] == 1))
    fn = length(which(x[,1] == 1 & x[,2] == 0))
    
    y <- matrix(
      c(tp, fp, fn, tn), 
      ncol = 2, nrow = 2
    )
  }
  else {
    y <- x
  }
  
  dimnames(y) = list(c("first_positive", "first_negative"), c("second_positive", "second_negative"))
  class(y) <- c("contingency.table", "rating.set", class(y))
  attr(y, "agreement") <- list(
    "true_positive" = y[1, 1],
    "false_negative" = y[1, 2],
    "false_positive" = y[2, 1],
    "true_negative" = y[2, 2]
  )
  attr(y, "baserate") <- baserateCT(y)
  attr(y, "kappa") <- kappaCT(y)
  
  return(y)
}

#' Convert codeset to contingency table
#'
#' @param x matrix contingency table (2x2)
#'
#' @return 2-column matrix representation of the contingency table
#' 
#' @export
as.code.set <- function(x) {
  y <- NULL
  if (all(dim(x) == c(2,2))) {
    y <- matrix(
      c(rep(1, x[1,1]*2), rep(1:0, x[1,2]), rep(0:1, x[2,1]), rep(0, x[2,2]*2)),
      byrow = TRUE,
      ncol = 2
    )
  } 
  else {
    y <- x
  }
  
  class(y) <- c("code.set", "rating.set", class(y))
  attr(y, "baserate") <- baserateSet(y)
  attr(y, "kappa") <- kappaSet(y)
  colnames(y) <- c("first_rater", "second_rater")
  
  return(y)
}

#' Helper function to return special values on a rating set
#' 
#' @param x Set or Contingency.Table
#'
#' @param i Value to search for
#'
#' @export
"$.rating.set" <- function (x, i) {
  if (i %in% names(attributes(x))) {
    return(attr(x, i))
  }
  else {
    if (is(x, "code.set")) {
      if (i %in% colnames(x)) {
        return(x[, which(colnames(x) == i)])
      }
    }
    else if (is(x, "contingency.table")) {
      if (i %in% colnames(x)) {
        return(sum(x[, which(colnames(x) == i)]))
      }
      else if (i %in% rownames(x)) {
        return(sum(x[which(rownames(x) == i), ]))
      }
    }
  }
  
  return(NULL)
}

#' @export
.DollarNames.rating.set <- function(x, pattern) {
  vals <- .codegen_dollar(x, pattern)
  
  cols <- colnames(x)
  if (!is.null(cols)) {
    vals <- c(cols, vals)
  }

  rows <- rownames(x)
  if (!is.null(rows)) {
    vals <- c(rows, vals)
  }
  
  vals
}

.codegen_dollar <- function(x, pattern = "") {
  attrs <- names(attributes(x))
  vals <- attrs[attrs %in% c("baserate", "kappa", "agreement")]
  return(vals)
}

#' @export
print.rating.set <- function(x, ...) {
  y <- unclass(x)

  for(a in additional_attributes) {
    attr(y, a) <- NULL
  }
  
  print(y, ...)
}

###
#' @title Rho (kappa)
#' 
#' @description 
#' This function calculates rho for an observed kappa value with associated set parameters
#' (testSetLength and OcSBaserate). Called by \code{\link{rho}}. A p-value is returned and if 
#' this value is less than 0.05, it is said that the handset does generalize to the entire set
#' 
#' @template rho-params
#' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
#' @param method set to "c" to calculate using the C++ implmentation. Defaults to "standard"
#' 
#' @export
#' @return rho for the given parameters
###
rhoK = function(
  x, 
  OcSBaserate,
  testSetLength,
  testSetBaserateInflation = 0, 
  OcSLength = 10000, 
  replicates = 800, 
  ScSKappaThreshold = 0.9, 
  ScSKappaMin = 0.40, 
  ScSPrecisionMin = 0.6, 
  ScSPrecisionMax = 1,
  method = "standard"
) {
  if (
    (x > 1 | x < 0) || 
    (OcSBaserate > 1 | OcSBaserate < 0)
  ) stop("x and OcSBaserate must be between 0 and 1.")

  if (OcSBaserate < 0 || ScSKappaThreshold < 0 || ScSKappaMin < 0) 
    stop("OcSBaserate, ScSKappaThreshold, and ScSKappaMin must be positive.")

  if (
    testSetLength %% 1 != 0 ||
    OcSLength %% 1 != 0 ||
    replicates %% 1 != 0
  ) stop("testSetLength, OcSLength, and replicates values must be a positive integer.")

  if (ScSKappaThreshold < ScSKappaMin) 
    stop("ScSKappaThreshold must be greater than ScSKappaMin.")

  if (
    (ScSPrecisionMin < 0 | ScSPrecisionMin > 1) ||
    (ScSPrecisionMax < 0 | ScSPrecisionMax > 1)
  ) stop("ScSPrecisionMin and ScSPrecisionMax must be between 0 and 1.")

  if (ScSPrecisionMax < ScSPrecisionMin) 
    stop("ScSPrecisionMax must be greater than ScSPrecisionMin.")

  result <- NULL
  if (method == "standard") {
    result <- calcRho(  x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
  }
  else {
    result <- calcRho_c(x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
  }
  return(result)
}

###
#' @title Get Handset
#' @description This function is to get a handset of a set and calculate the kappa
#' @param set This is the set to take a handset of
#' @param handSetLength This is the length of the handset to take
#' @param handSetBaserate This is the minimum baserate to inflate the handset to
#' @param returnSet If TRUE, then return the handSet if FALSE, return the kappa of the handSet
#' 
#' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
#' @export
###
getHandSet = function(set, handSetLength, handSetBaserate, returnSet=FALSE) {
  #positives is the minimum number of positive pairs in the first rater
  positives = ceiling(handSetLength * handSetBaserate);
  posInd = which(set[,1] == 1);

  if(positives > length(posInd)){stop("Not enough positives in first rater to inflate to this level")}
  #if there is an inflated baserate positives will be > 0, so if handSetBaserate > 0, positives > 0
  if (positives > 0) {
    #goes through and picks a set that fits the number of positives first then randomly samples from the rest
    #if the set created is not a valid set, then it reruns the process until a valid set is found
    positiveIndices = posInd[sample.int(length(posInd),size=positives,replace=FALSE)];
    others = set[!(1:nrow(set) %in% positiveIndices),];
    otherIndices = sample.int(nrow(others),size=(handSetLength - positives),replace=FALSE);
    this.set = rbind(set[positiveIndices,], others[otherIndices,]);
  } else if (positives == 0){
    #if there is no restriction on positives, then the set is generated by randomly sampling from the entire set
    #this.set=matrix(0,2,2);
    #while ((sum(this.set[,2]) == 0 | sum(this.set[,2]) == nrow(this.set))) {
      theseIndices = sample.int(nrow(set),size=handSetLength,replace=FALSE);
      this.set = set[theseIndices,];
    #}
  }

  #if return set is true return set rather than kappa
  if (returnSet) {
    this.set <- as.code.set(this.set);
    return(this.set);
  }

  #otherwise return the kappa of the set
  handSetKappa = calcKappa(this.set);
  
  return(handSetKappa);
}

###
#' @title Get Handset
#' @description This function is to get a handset of a set and calculate the kappa
#' 
#' @param full.ct This is the set to take a handset of
#' @param handSetLength This is the length of the handset to take
#' @param handSetBaserate This is the minimum baserate to inflate the handset to
#' @param as_kappa If FALSE then return the handSet, if TRUE (default) return the kappa of the handSet
#' 
#' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
#' @export
###
getHandCT <- function(full.ct, handSetLength, handSetBaserate, as_kappa = TRUE) {
  positives <- ceiling(handSetLength * handSetBaserate)
  if (positives > sum(full.ct[1,])) {
    stop("Not enough positives in first rater to inflate to this level")
  }
  
  # positiveIndices <- sample_ct(full.ct[1,], positives)
  # 
  # other_ct <- full.ct
  # other_ct[1, ] <- other_ct[1, ] - positiveIndices
  # 
  # other_ct_vec <- as.vector(other_ct)
  # this_ct <- sample_ct(other_ct_vec,  handSetLength - positives)
  # this_ct[1, ] <- this_ct[1, ] + positiveIndices
  
  this_ct <- getHand_ct(ct = full.ct, handSetLength = handSetLength, handSetBaserate = handSetBaserate)

  if (!as_kappa) {
    this_ct <- as.contingency.table(this_ct);
    return(this_ct)
  } 

  return(kappa_ct(this_ct))
}

# sample_ct <- function(x, size) {
#   x_vec <- as.vector(x)
#   x_sampled <- matrix(rep(0, length(x)), ncol = length(x_vec) / 2)
#   for(i in seq(size)) {
#     s <- sample(seq(x_vec), 1, replace = F, prob = x_vec / sum(x_vec))
#     x_sampled[[s]] = x_sampled[[s]] + 1
#     x_vec[s] = x_vec[s] - 1
#   }
# 
#   return(x_sampled)
# }

###
#' @title Rho (contingency Table)
#' 
#' @description 
#' This function calculates rho and kappa for a given \code{\link{contingencyTable}}, and returns a list containing both values. Called by \code{\link{rho}}.
#' 
#' 
#' @template rho-params
#' 
#' @export
#' @return A list of the format:\describe{
#' \item{rho}{The rho of the \code{\link{contingencyTable}}}
#' \item{kappa}{The Cohen's Kappa of the \code{\link{contingencyTable}}}
#' }
###
rhoCT = function(
  x,
  OcSBaserate = NULL,
  testSetBaserateInflation = 0,
  OcSLength = 10000,
  replicates = 800,
  ScSKappaThreshold = 0.9,
  ScSKappaMin = .40,
  ScSPrecisionMin = 0.6,
  ScSPrecisionMax = 1,
  method = "standard"
){
  if (any(x < 0)) {
    stop("Values in Contingency Table must be positive")
  }
  kappa = kappa_ct(x);
  
  if (is.null(OcSBaserate)) {
    OcSBaserate = (sum(x[1,])) / (sum(x))
  }
  
  rho_result <- rhoK(
    kappa, 
    OcSBaserate,
    sum(x),
    testSetBaserateInflation,
    OcSLength,
    replicates,
    ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax
  )
  
  return(list(
    rho = rho_result, 
    kappa = kappa,
    recall = rating_set_recall(x),
    precision = rating_set_precision(x)
  ))
}

###
#' @title Rho using a file
#' 
#' @description 
#' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}} as defined by the file 
#' and columns (col1, col2), and returns a list containing both values. Called by \code{\link{rho}}.
#' 
#' @export
#' 
#' @param col1 The first column from file
#' @param col2 The second column from file
#' @template rho-params
#' 
#' @return A list of the format:\describe{
#' \item{rho}{The rho of the \code{\link{codeSet}}}
#' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
#' }
###
rho.file <- function(
  x, col1, col2, 
  OcSBaserate = NULL,
  testSetBaserateInflation = 0,
  OcSLength = 10000,
  replicates = 800,
  ScSKappaThreshold = 0.9,
  ScSKappaMin = .40,
  ScSPrecisionMin = 0.6,
  ScSPrecisionMax = 1,
  method = "standard"
) {
  setFile = utils::read.csv(x)
  set = as.matrix(setFile[,c(col1, col2)])
  
  if (!any(is.na(set))) {
    rhoSet(
      set,
      OcSBaserate,
      testSetBaserateInflation, 
      OcSLength, replicates, 
      ScSKappaThreshold, ScSKappaMin, 
      ScSPrecisionMin, ScSPrecisionMax
    )
  }
  else {
    stop("The columns provided for col1 and col2 must not contain NA")
  }
}

### 
#' @title Rho
#' 
#' @description 
#' This function calculates rho for a \code{\link[=getTestSet]{testSet}}, \code{\link{contingencyTable}}, or an observed kappa value with associated set parameters (testSetLength and OcSBaserate). 
#' 
#' @details 
#' Rho is a Monte Carlo rejective method of interrater reliability statistics, implemented here for Cohen's Kappa. Rho constructs a collection of data sets in which kappa is below a specified threshold, and computes the empirical distribution on kappa based on the specified sampling procedure. Rho returns the percent of the empirical distribution greater than or equal to an observed kappa. As a result, Rho quantifies the type 1 error in generalizing from an observed test set to a true value of agreement between two raters.
#' 
#' Rho starts with an observed kappa value, calculated on a subset of a \code{\link{codeSet}}, known as an observed \code{\link[=getTestSet]{testSet}}, and a \emph{kappa threshold} which indicates what is considered significant agreement between raters.
#' 
#' It then generates a collection of fully-coded, simulated 
#' \code{\link[=getTestSet]{codeSets}} (ScS), further described in 
#' \code{\link{createSimulatedCodeSet}}, all of which have a kappa value below the kappa 
#' threshold and similar properties as the original \code{\link{codeSet}}.
#' 
#' Then, kappa is calculated on a \code{\link[=getTestSet]{testSet}} sampled from each of the ScSs in the
#' collection to create a null hypothesis distribution. These \code{\link[=getTestSet]{testSets}} mirror the observed \code{\link[=getTestSet]{testSets}} in their size and sampling method.  How these \code{\link[=getTestSet]{testSets}} are sampled is futher described in \code{\link{getTestSet}}.
#' 
#' The null hypothesis is that the observed \code{\link[=getTestSet]{testSet}}, was sampled from a data set, which, if both raters were to code in its entirety, would result in a level of agreement below the kappa threshold.
#' 
#' For example, using an alpha level of 0.05, if the observed kappa is greater than 95 percent of the kappas in the null hypothesis distribution, the null hypothesis is rejected. Then one can conclude that the two raters would have acceptable agreement had they coded the entire data set.
#'
#' @template rho-params
#' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
#' 
#' @examples
#' # Given an observed kappa value
#' rho(x = 0.88, OcSBaserate = 0.2, testSetLength = 80)
#' 
#' # Given a test Set
#' rho(x = codeSet)
#' 
#' # Given a contingency Table
#' rho(x = contingencyTable)
#'
#' @export
#' @return rho and kappa for the given data and parameters (unless kappa is given)
###
rho = function(
  x,
  OcSBaserate = NULL,
  testSetLength = NULL,
  testSetBaserateInflation = 0,
  OcSLength = 10000,
  replicates = 800,
  ScSKappaThreshold = 0.9,
  ScSKappaMin = 0.40,
  ScSPrecisionMin = 0.6,
  ScSPrecisionMax = 1,
  method = "standard"
) {

  ### Calculate rho given a kappa value
  if (is(x, "numeric")) {
    if (is.null(OcSBaserate)) { stop("Must give a baserate when a kappa value is given") }
    if (is.null(testSetLength)) { stop("Must give a testSetLength when a kappa value is given") }

    result <- rhoK(x, OcSBaserate,
               testSetLength,
               testSetBaserateInflation,
               OcSLength,
               replicates,
               ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax) 
  }

  ### Calculate rho given a code set
  else if (nrow(x) != 2) {
    result <- rhoSet(x, OcSBaserate,
                  testSetBaserateInflation,
                  OcSLength,
                  replicates,
                  ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
  }

  ### Calculate rho given a contingency table
  else {
    result <- rhoCT(x, OcSBaserate,
                 testSetBaserateInflation,
                 OcSLength,
                 replicates,
                 ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
  }

  return(result)
}

genPKcombo = function(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate){
  # currKappa = sample(kappaDistribution, 1, replace = FALSE, prob = kappaProbability)
  # currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
  currKappa = kappaDistribution[sample.int(length(kappaDistribution), 1, replace = FALSE, prob = kappaProbability)]
  currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]

  if(!checkBRPKcombo(baserate, currPrecision, currKappa)){
    precisionMin = (2*baserate*currKappa - 2*baserate - currKappa)/(currKappa - 2)
    indicies = which(precisionDistribution > precisionMin)
    
    if (length(indicies) == 0) {
      return(genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate))
    }
    
    precisionDistribution = precisionDistribution[indicies]
    precisionProbability = precisionProbability[indicies]
    # currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
    currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]
    return(list(precision = currPrecision, kappa = currKappa))
  } else {
    return(list(precision = currPrecision, kappa = currKappa))
  }
}

###
# Generate Precision Combo
#
# This function generates a precision that fits within the range and will create a valid set
# @param kappa This is the kappa of the desired pair
# @param Pmin This is the minimum precision desired
# @param Pmax This is the maximum precision desired
# @param BR This is the baserate of the desired pair
# @return This returns a precision value that will generate a valid set with the other parameters
####

genPcombo = function(kappa, Pmin, Pmax, BR){
  #gets a random precision in the precision range
  currPrecision = stats::runif(1, Pmin, Pmax);
  #checks is the combination of precision, kappa, and baserate is valid.  If it is vaid, return the precision, if not, recall the function until a valid precision is found
  if(!checkBRPKcombo(BR, currPrecision, kappa)){
    genPcombo(kappa, Pmin, Pmax, BR);
  }else{
    return(currPrecision);
  }
}

###
#' @title Calculate Baserate (CT)
#' 
#' @description
#' This function calculates the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
#' 
#' @export
#' 
#' @param CT The \code{\link{contingencyTable}} for which the baserate is calculated
#' 
#' @seealso \code{\link{baserate}} and \code{\link{baserateSet}}
#' 
#' @return A list of the format:\describe{
#' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
#' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
#' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
#' }
#' 
###

baserateCT = function(CT){
  if(any(CT < 0)){stop("Values in Contingency Table must be positive")}
  firstBaserate = sum(CT[1,])/sum(CT)
  secondBaserate = sum(CT[,1])/sum(CT)
  averageBaserate = mean(c(firstBaserate,secondBaserate))
  return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
}

###
# Get Boot P Value
#
# This function gets the p value of the result in the distribution
# @param distribution The distribution that the result falls into
# @param result The value that you are evaluating where in the distribution it falls
# @return This returns the p value from the result into the distribution
###
getBootPvalue = function(distribution, result) {
  #returns the percentage of the time that the distribution was greater or equal to the observed kappa
  N = length(distribution);
  if(result < mean(distribution, na.rm = T)) {
    #if the result is less than the mean of the distribution, than the p value is 1
    return(1);
  } else {
    #return the number of times that the distribution is greater than the result as a percentage of the total number of items in the distribution
    return(length(which(distribution >= result)) / N);
  }
}

#' Generate a Handset
#' 
#' @description Generate a vector representing indices of set, using the handSetBaserate
#' to determine the minimum number of indices that are positive
#' 
#' @param set matrix of two columns
#' @param handSetLength number of indices to find
#' @param handSetBaserate number between 0 and 1 to use as a minimum number of positive indices
#'
#' @return vector of indices from set
#' @export
getHandSetIndices = function(set, handSetLength = 20, handSetBaserate = 0.2) {
  positives = ceiling(handSetLength * handSetBaserate);
  posInd = which(set[, 1] == 1);
  
  if (positives > length(posInd)) {
    stop("Not enough positives in first rater to inflate to this level")
  }
  
  this.set = NULL
  if (positives > 0) {
    positiveIndices = posInd[sample.int(length(posInd), size=positives, replace=FALSE)]
    others = which(!(1:nrow(set) %in% positiveIndices))
    otherIndices = sample.int(length(others), size = (handSetLength - positives), replace = FALSE)
    this.set = c(positiveIndices, others[otherIndices])
  }
  else if (positives == 0) {
    theseIndices = sample.int(length(set),size=handSetLength,replace=FALSE)
    this.set = theseIndices
  }
  
  return(sample(this.set))
}

###
#' @title Calculate Baserate (Set)
#' 
#' @description
#' This function will calculate the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
#' 
#' @param set The \code{\link{codeSet}} for which the baserate is calculated
#' 
#' @seealso \code{\link{baserate}} and \code{\link{baserateCT}}
#' 
#' @export
#' @return A list of the format:\describe{
#' \item{firstBaserate}{The percentage that a positive code, or a 1, appears in the first rater}
#' \item{secondBaserate}{The percentage that a positive code, or a 1, appears in the second rater}
#' \item{averageBaserate}{The average percentage that a positive code, or a 1, appears in either of the two raters}
#' }
#' 
###
baserateSet = function(set){
  firstBaserate = length(which(set[,1] == 1))/nrow(set)
  secondBaserate = length(which(set[,2] == 1))/nrow(set)
  averageBaserate = mean(c(firstBaserate,secondBaserate))
  return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
}

###
# Calculate Kappa
#
# This function calculates kappa given a set
# @param set This is the set to calculate kappa from
# @param isSet TRUE if set FALSE if contigency table
# @param kappaThreshold if null then return kappa otherwise return list of kappa and above or not
# @return This returns the kappa of the set given
###
calcKappa <- function(set, isSet = TRUE, kappaThreshold = NULL) {
  if (!isSet) {
    set <- contingencyToSet(set[1, 1], set[2, 1], set[1, 2], set[2, 2])
  }

  if (
    length(which(set[, 1] == set[, 2] & set[, 1] == 1)) == nrow(set) |
    length(which(set[, 1] == set[, 2] & set[, 1] == 0)) == nrow(set)
  ) {
    return(1)
  }

  #calculates kappa by calculating the agreement and the baserates and then creating the adjacenty matrix
  agreement <- length(which(set[, 1] == set[, 2])) / nrow(set);
  baserate2 <- length(which(1 == set[, 2])) / nrow(set);
  baserate1 <- length(which(1 == set[, 1])) / nrow(set);
  randomAgreement <- baserate1 * baserate2 + (1 - baserate1) * (1 - baserate2);
  kappa <- (agreement - randomAgreement) / (1 - randomAgreement);

  if (is.null(kappaThreshold)) {
    return(kappa);
  } else {
    return(list(kappa = kappa, above = kappa > kappaThreshold))
  }
}

###
# Check Baserate, Precision, Kappa Combo
#
# This function checks to make sure the combination of BR, P, and K will create a valid set
# @param BR This is the supplied baserate
# @param P This is the supplied precision
# @param K This is the supplied kappa
# @return returns TRUE if the combination is valid, FALSE otherwise
###
checkBRPKcombo = function(BR, P, K) {
  #if right is less than P, then it is a valid combination. If not, then the pair of kappa, precision, and baserate does not correspond to a valid set
  right = (2*BR*K - 2*BR - K)/(K - 2);
  if(P > right){
    return(TRUE);
  }else{
    return(FALSE);
  }
}

###
# This function calculates rho given the parameters of the initial set and the kappa of the observed handset and is
# called by rhoK()
#
# The computation is libqe-stats' qestats::rho_standard (src/rhor.cpp): the
# same draws and arithmetic as this function's former R implementation
# (generateKPs, getR, prset, getHandSet, getBootPvalue), so seeded results
# are unchanged. KPs, when given, is a list of genPKcombo() results.
###
calcRho = function(
  x, 
  OcSBaserate, 
  testSetLength, 
  testSetBaserateInflation = 0,
  OcSLength = 10000, 
  replicates = 800, 
  ScSKappaThreshold = 0.9, 
  ScSKappaMin = 0.40, 
  ScSPrecisionMin = 0.6, 
  ScSPrecisionMax = 1, 
  KPs = NULL
) {
  if (!is.null(KPs)) {
    KPs = cbind(
      vapply(KPs, function(KP) KP$precision, numeric(1)),
      vapply(KPs, function(KP) KP$kappa, numeric(1))
    )
  }

  return(calcRho_standard_c(
    x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates,
    ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax, KPs
  ));
}

rating_set_recall <- function (x) {
  if (!is(x, "contingency.table")) {
    x <- as.contingency.table(x)
  }
  
  x$agreement$true_positive / (x$agreement$true_positive + x$agreement$false_negative)
}

rating_set_precision <- function (x) {
  if (!is(x, "contingency.table")) {
    x <- as.contingency.table(x)
  }
  
  x$agreement$true_positive / (x$agreement$true_positive + x$agreement$false_positive)
}

###
# Create A Random Set
#
# This creates a set with given parameters
# @param setLength The length of the set to be created
# @param baserate The baserate of the set to be created
# @param kappa The kappa of the set to be created
# @param minPrecision The minimum precision accepted in this set. Defaults to 0
# @param maxPrecision The maximum precision accepted in this set. Defaults to 1
# @return This returns a set with the given parameters
###
createRandomSet = function(setLength, baserate, kappaMin, kappaMax, minPrecision = 0, maxPrecision = 1, type = "set"){
  #generates a valid precision with the given kappa and baserate
  KP = generateKPs(1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision)[[1]]
  kappa = KP$kappa
  precision = KP$precision
  recall = getR(kappa, baserate, precision)
  
  if(type == "ct") {
    fullSet = pr.ct(precision = precision, recall = recall, length = setLength, baserate = baserate);
  } 
  else {
    fullSet = prset(precision = precision, recall = recall, length = setLength, baserate = baserate);
  }
  
  return(fullSet);
}

###
#' @title Rho Min
#' 
#' @description 
#' This function calculates the minimum testSetLength where it is possible to get a rho less than alpha for the given parameters of rho.
#' 
#' @param baserate A \code{\link{baserate}}
#' @param alpha The threshold of significance for rho (similar to an alpha level for a p value), defaulted to 0.05
#' @param inc An integer indicating by how much the testSetLength should increase each iteration
#' @param printInc A boolean indicating whether to print out each increment value with it's corresponding significance for rho
#' @param ... Any additional parameters passed into \code{\link{rho}}
#' 
#' @examples
#' #Add testSetBaserateInflation as an additional parameter
#' rhoMin(0.2, testSetBaserateInflation = 0.33)
#' 
#' #Add testSetBaserateInflation as well as changing inc and selecting printInc
#' rhoMin(0.2, inc = 5, printInc = TRUE, testSetBaserateInflation = 0.33)
#' 
#' @export
#' @return  The minimum length of testSet, to the nearest multiple of inc, greater than the minimum length, that would give a value where rho less than alpha becomes mathematically possible.
###
rhoMin = function(baserate, alpha = 0.05, inc = 10, printInc = FALSE, ...){
  rhoVal = 1
  testSetLength = 0
  if(inc%%1 != 0 | inc <= 0) stop("Inc value must be a positive integer.")
  if(alpha > 1 | alpha < 0) stop("Alpha must be between 0 and 1.")
  
  
  while(rhoVal > alpha){
    testSetLength = testSetLength + inc
    rhoVal = rho(1, baserate, testSetLength, ...)
    if(printInc) {
      cat(testSetLength, rhoVal,"\n")
    }
  }
  return(testSetLength)
  
}

###
#' @title Calculate Baserate
#' 
#' @description
#' This function calculates the baserate of the first rater, second rater, and the average of both the raters.
#' 
#' @details
#' A baserate is the percentage, as a decimal, that a positive code appears in data (either a \code{\link{codeSet}} or \code{\link{contingencyTable}}) for a given rater.  It is assumed that the first rater is more experienced and thus provides a better estimation of the actual baserate for a given code, so the first rater's baserate is often used as if it is the actual baserate.  If the raters are assumed to have the same experience level, the average baserate may give a better estimation.  If the second rater is more experienced, the second rater's baserate may give a better estimation.  Functions assume that the first rater is the more experienced rater and thus uses the first rater's baserate as the overall baserate estimation.
#' 
#' 
#' @param data The \code{\link[=getTestSet]{testSet}} or \code{\link{contingencyTable}} for which the baserate is calculatede
#' 
#' @seealso \code{\link{baserateSet}} and \code{\link{baserateCT}}
#'
#' @examples
#' #Given a code set
#' baserate(data = codeSet)
#' 
#' #Given a contingency Table
#' baserate(data = contingencyTable)
#'
#' @export
#' @return A list of the format:\describe{
#' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
#' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
#' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
#' }
#' 
###
baserate = function(data){
  if (nrow(data) != 2) {
    return(baserateSet(data))
  } else {
    return(baserateCT(data))
  }
}

###
#' @title Calculate kappa (contingency Table)
#' 
#' @description 
#' This function calculates Cohen's kappa on a \code{\link{contingencyTable}}. Called by \code{\link{kappa}}. 
#' 
#' @param ct A \code{\link{contingencyTable}}
#' 
#' @seealso \code{\link{kappa}} and \code{\link{kappaSet}}
#' 
#' @export
#' @return  The kappa of the \code{\link{contingencyTable}}
###
kappaCT = function(ct){
  if(any(ct < 0)){stop("Values in Contingency Table must be positive")}
  if(ncol(ct) != 2 || nrow(ct) != 2) stop("Incorrect number of dimensions: Contingency table must be 2x2.")
  for(i in 1:nrow(ct)) {
    for(j in 1:ncol(ct)) {
      if (ct[i,j]%%1 != 0) stop("Contingency table values must be positive integers.")
    }  
  }
  
  return(calcKappa(ct,isSet = FALSE,kappaThreshold = NULL));
}

###
# Get Recall
#
# This gets the recall from a pair
# @param kappa This is the kappa of the pair
# @param BR This is the baserate of the pair
# @param P this is the precision of the pair
# @return returns the recall given the parameters
###
getR = function(kappa, BR, P){
  #gets recall that corresponds with the kappa, baserate, and precision
  top = kappa*P;
  R = top / (2*P-2*BR-kappa+2*BR*kappa);
  return (R);
}

###
#' @title Calculate kappa (Set)
#' 
#' @description 
#' This function calculates Cohen's kappa for a given \code{\link{codeSet}}. Called by \code{\link{kappa}}.
#' 
#' @export
#' @param set A \code{\link{codeSet}}
#' 
#' @seealso \code{\link{kappa}} and \code{\link{kappaCT}}
#' 
#' @return 
#' The kappa of the \code{\link{codeSet}}
#' 
###
kappaSet = function(set){
  if(any(set != 1 & set != 0)) {
    stop("Each set must only consist of zeros and ones.")
  }
  
  return(calcKappa(set,isSet = TRUE,kappaThreshold = NULL));
}


###
#' @include rhoK.R
#' @title Rho (set)
#' 
#' @description 
#' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}}, and returns a list containing both values. Called by \code{\link{rho}}.
#' 
#' @export
#' 
#' @template rho-params
#' 
#' @return A list of the format:\describe{
#' \item{rho}{The rho of the \code{\link{codeSet}}}
#' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
#' }
###
rhoSet = function(
  x,
  OcSBaserate = NULL,
  testSetBaserateInflation = 0,
  OcSLength = 10000,
  replicates = 800,
  ScSKappaThreshold = 0.9,
  ScSKappaMin = 0.40,
  ScSPrecisionMin = 0.6,
  ScSPrecisionMax = 1,
  method = "standard"
) {
  kappa = kappaSet(x);
  
  row.count = nrow(x);
  if(is.null(OcSBaserate)) {
    OcSBaserate = length(which(x[,1] == 1)) / row.count;
  }

  return(list(
    rho = rhoK(
      kappa, 
      OcSBaserate,
      row.count,
      testSetBaserateInflation,
      OcSLength,
      replicates,
      ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax,
      method = method
    ), 
    kappa = kappa,
    recall = rating_set_recall(x),
    precision = rating_set_precision(x)
  ))
}

###
#' @title Calculate kappa
#'
#' @description
#' This function calculates Cohen's kappa on a \code{\link{contingencyTable}} or a \code{\link{codeSet}}
#'
#' @param data A \code{\link{contingencyTable}} or a \code{\link{codeSet}}
#'
#'
#' @seealso \code{\link{kappaSet}} and \code{\link{kappaCT}}
#'
#' @examples
#' #Given a code set
#' kappa(data = codeSet)
#'
#' #Given a contingency Table
#' kappa(data = contingencyTable)
#'
#' @export
#' @return The kappa of the \code{\link{contingencyTable}} or \code{\link{codeSet}}
###
kappa <- function(data) {
  if (all(dim(data) == 2)) {
    return(kappaCT(data))
  } else {
    return(kappaSet(data))
  }
}

#pragma once
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>

// Functions that must reproduce R's arithmetic (rhoR parity) start their body
// with QESTATS_R_ARITHMETIC and carry QESTATS_R_ARITHMETIC_FN: R never fuses
// a * b + c, so these keep the compiler from contracting it into an FMA.
#if defined(__clang__)
#  define QESTATS_R_ARITHMETIC _Pragma("STDC FP_CONTRACT OFF")
#  define QESTATS_R_ARITHMETIC_FN
#elif defined(__GNUC__)
#  define QESTATS_R_ARITHMETIC
#  define QESTATS_R_ARITHMETIC_FN __attribute__((optimize("fp-contract=off")))
#else
#  define QESTATS_R_ARITHMETIC
#  define QESTATS_R_ARITHMETIC_FN
#endif

namespace qestats {

// Confusion matrix proportions (must sum to 1 after normalization).
struct ConfusionMatrix {
    double tp = 0, fp = 0, fn = 0, tn = 0;
};

// All derived classification metrics.
struct Metrics {
    // confusion matrix proportions
    double tp = 0, fp = 0, fn = 0, tn = 0;
    // derived
    double b1 = 0;        // classifier positive rate  (tp + fp)
    double b2 = 0;        // true positive rate        (tp + fn)
    double precision = 0; // tp / (tp + fp)
    double recall    = 0; // tp / (tp + fn)
    double accuracy  = 0; // tp + tn
    double fdr       = 0; // fp / (tp + fp)   — False Discovery Rate
    double for_      = 0; // fn / (fn + tn)   — False Omission Rate
    double kappa     = 0;
};

// Range result — min/max of a metric over the key variable interval.
struct MetricRange {
    double min = 0;
    double max = 0;
};

// Confidence interval from a hypergeometric estimate.
struct Interval {
    double lower = 0;
    double mean  = 0;
    double upper = 0;
};

// One row in the metrics_range output table.
struct MetricRangeRow {
    std::string metric;
    double lower = 0, mean = 0, upper = 0;
};

// ── Core arithmetic ───────────────────────────────────────────────────────────

// Cohen's kappa from confusion-matrix counts (or proportions), computed in
// the same order as rhoR's calcKappa() so the two agree to the bit: perfect
// agreement on one class (all TP or all TN) is 1.
QESTATS_R_ARITHMETIC_FN
inline double kappa_from_counts(double tp, double fp, double fn, double tn) {
    QESTATS_R_ARITHMETIC
    const double s = tp + fp + fn + tn;
    if (tp == s || tn == s) return 1.0;
    const double agreement = (tp + tn) / s;
    const double rate1 = (tp + fn) / s;   // first rater's positive rate
    const double rate2 = (tp + fp) / s;   // second rater's positive rate
    const double pe = rate1 * rate2 + (1 - rate1) * (1 - rate2);
    return (pe < 1.0) ? (agreement - pe) / (1 - pe) : 0.0;
}

// Rates are computed from the counts, not the normalized cells, so precision
// and recall match rhoR's TP / (TP + FP) and TP / (TP + FN) exactly.
inline Metrics metrics_from_cm(const ConfusionMatrix& cm) {
    const double s = cm.tp + cm.fp + cm.fn + cm.tn;
    if (s <= 0) throw std::domain_error("metrics_from_cm: confusion matrix is empty");
    Metrics m;
    m.tp = cm.tp / s; m.fp = cm.fp / s;
    m.fn = cm.fn / s; m.tn = cm.tn / s;
    const double pos1 = cm.tp + cm.fp, pos2 = cm.tp + cm.fn, neg1 = cm.fn + cm.tn;
    m.b1        = pos1 / s;
    m.b2        = pos2 / s;
    m.precision = (pos1 > 0) ? cm.tp / pos1 : 0.0;
    m.recall    = (pos2 > 0) ? cm.tp / pos2 : 0.0;
    m.accuracy  = (cm.tp + cm.tn) / s;
    m.fdr       = (pos1 > 0) ? cm.fp / pos1 : 0.0;
    m.for_      = (neg1 > 0) ? cm.fn / neg1 : 0.0;
    m.kappa     = kappa_from_counts(cm.tp, cm.fp, cm.fn, cm.tn);
    return m;
}

// Convenience: from raw counts (negative values clamped to 0).
inline Metrics calculate_metrics(double tp, double fp, double fn, double tn) {
    if (tp < 0) tp = 0; if (fp < 0) fp = 0;
    if (fn < 0) fn = 0; if (tn < 0) tn = 0;
    return metrics_from_cm({tp, fp, fn, tn});
}

// Retrieve a single metric value by name (supports alias names).
inline double get_metric(const Metrics& m, const std::string& name) {
    if (name == "b1" || name == "baserate")  return m.b1;
    if (name == "b2" || name == "baserate2") return m.b2;
    if (name == "precision")  return m.precision;
    if (name == "recall")     return m.recall;
    if (name == "accuracy")   return m.accuracy;
    if (name == "FDR")        return m.fdr;
    if (name == "FOR")        return m.for_;
    if (name == "kappa")      return m.kappa;
    if (name == "tp")         return m.tp;
    if (name == "fp")         return m.fp;
    if (name == "fn")         return m.fn;
    if (name == "tn")         return m.tn;
    throw std::invalid_argument("get_metric: unknown metric '" + name + "'");
}

} // namespace qestats

#pragma once
#include "metrics.hpp"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

// Rho (Shaffer's rho): how likely a test set's kappa this high would be if the
// coding process's population kappa were below a threshold. A port of rhoR:
// given an adapter over R's RNG, a seeded rho() returns exactly what rhoR does.
//
// rhoR computes rho two ways, and both are kept, quirks included:
//
//   RhoMethod::Standard   rhoR's default path, written in R (calcRho): grids
//                         from seq(), draws from sample.int(), the hand set
//                         sampled as rows of the full set. qestats::rho_standard.
//   RhoMethod::Fast       rhoR's C++ path (calcRho_c, any other `method`):
//                         Armadillo regspace() grids, RcppArmadillo sample() /
//                         randi() draws, the hand set sampled from the
//                         contingency table. qestats::rho_fast.
//
// Random numbers come from an Rng: any type with
//
//   double unif_rand();          // uniform on (0, 1)          — R's unif_rand()
//   double unif_index(double n); // uniform integer in [0, n)  — R's R_unif_index()
//
// The Standard method draws only through unif_index, the Fast method only
// through unif_rand. StdRng is a standalone default.

namespace qestats {

// mt19937_64 with the draws derived by hand (the std distributions differ
// between standard libraries), so a seed gives the same rho everywhere.
struct StdRng {
    std::mt19937_64 gen;
    explicit StdRng(std::uint64_t seed = 5489u) : gen(seed) {}
    double unif_rand() {
        double u;
        do { u = static_cast<double>(gen() >> 11) * 0x1.0p-53; } while (u <= 0.0);
        return u;
    }
    // Rejection on the top bits, like R_unif_index.
    double unif_index(double n) {
        const auto m = static_cast<std::uint64_t>(n);
        if (m <= 1) return 0;
        int bits = 0;
        while (bits < 64 && (std::uint64_t(1) << bits) < m) ++bits;
        std::uint64_t v;
        do { v = gen() >> (64 - bits); } while (v >= m);
        return static_cast<double>(v);
    }
};

enum class RhoMethod { Standard, Fast };

// rhoR's rhoK() arguments after the observed kappa, base rate and test set
// length, with the same defaults.
struct RhoOptions {
    double test_set_baserate_inflation = 0;
    int    ocs_length          = 10000;   // OcSLength
    int    replicates          = 800;
    double scs_kappa_threshold = 0.9;     // ScSKappaThreshold
    double scs_kappa_min       = 0.40;    // ScSKappaMin
    double scs_precision_min   = 0.6;     // ScSPrecisionMin
    double scs_precision_max   = 1;       // ScSPrecisionMax
    RhoMethod method = RhoMethod::Standard;
};

struct KappaPrecision {
    double kappa = 0, precision = 0;
};

// Both rhoR paths recurse until a valid (kappa, precision) pair turns up and
// give up with this message; here the retries are capped instead.
inline constexpr int RHO_MAX_ATTEMPTS = 10000;
inline const char* const RHO_NO_VALID_SET =
    "Could not generate a valid set given ranges of kappa and precision";

// ── Standard method: rhoR's R implementation ─────────────────────────────────

namespace rho_standard {

// getR(): the recall that gives `kappa` at this base rate and precision.
QESTATS_R_ARITHMETIC_FN
inline double recall(double kappa, double baserate, double precision) {
    QESTATS_R_ARITHMETIC
    const double top = kappa * precision;
    return top / (2 * precision - 2 * baserate - kappa + 2 * baserate * kappa);
}

// checkBRPKcombo(): whether (base rate, precision, kappa) can form a set.
QESTATS_R_ARITHMETIC_FN
inline bool valid_combo(double baserate, double precision, double kappa) {
    QESTATS_R_ARITHMETIC
    const double right = (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
    return precision > right;
}

// seq(from, to, length.out = n).
QESTATS_R_ARITHMETIC_FN
inline std::vector<double> seq(double from, double to, int n) {
    QESTATS_R_ARITHMETIC
    if (n <= 0) return {};
    if (n <= 2) {
        std::vector<double> v = {from, to};
        v.resize(n);
        return v;
    }
    if (from == to) return std::vector<double>(n, from);
    const double by = (to - from) / (n - 1);
    std::vector<double> v(n);
    v[0] = from;
    for (int i = 1; i < n - 1; ++i) v[i] = from + static_cast<double>(i) * by;
    v[n - 1] = to;
    return v;
}

// sample.int(n, k) without replacement, 0-based.
template <class Rng>
std::vector<std::size_t> sample_int(std::size_t n, std::size_t k, Rng& rng) {
    if (k > n)
        throw std::invalid_argument(
            "cannot take a sample larger than the population when 'replace = FALSE'");
    std::vector<std::size_t> out(k);
    if (k < 2) {
        for (std::size_t i = 0; i < k; ++i)
            out[i] = static_cast<std::size_t>(rng.unif_index(static_cast<double>(n)));
        return out;
    }
    std::vector<std::size_t> x(n);
    for (std::size_t i = 0; i < n; ++i) x[i] = i;
    for (std::size_t i = 0; i < k; ++i) {
        const auto j = static_cast<std::size_t>(rng.unif_index(static_cast<double>(n)));
        out[i] = x[j];
        x[j] = x[--n];
    }
    return out;
}

// genPKcombo(): one (kappa, precision) draw from the two grids. When the pair
// is invalid, precision is redrawn from the grid values above the minimum
// valid precision; when none are, kappa is redrawn too.
template <class Rng>
QESTATS_R_ARITHMETIC_FN
KappaPrecision kappa_precision(const std::vector<double>& kappas,
                               const std::vector<double>& precisions,
                               double baserate, Rng& rng) {
    QESTATS_R_ARITHMETIC
    for (int attempt = 0; attempt < RHO_MAX_ATTEMPTS; ++attempt) {
        const double kappa = kappas[sample_int(kappas.size(), 1, rng)[0]];
        double precision = precisions[sample_int(precisions.size(), 1, rng)[0]];
        if (valid_combo(baserate, precision, kappa)) return {kappa, precision};
        const double precision_min =
            (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
        std::vector<double> above;
        for (double p : precisions) if (p > precision_min) above.push_back(p);
        if (above.empty()) continue;
        precision = above[sample_int(above.size(), 1, rng)[0]];
        return {kappa, precision};
    }
    throw std::runtime_error(RHO_NO_VALID_SET);
}

// generateKPs() with the flat distribution: kappa from a 100,000-point grid,
// precision from a 10,000-point grid.
template <class Rng>
std::vector<KappaPrecision>
generate_kps(int n, double baserate, double kappa_min, double kappa_max,
             double precision_min, double precision_max, Rng& rng) {
    const auto kappas = seq(kappa_min, kappa_max, 100000);
    const auto precisions = seq(precision_min, precision_max, 10000);
    std::vector<KappaPrecision> kps;
    kps.reserve(n > 0 ? n : 0);
    for (int i = 0; i < n; ++i) kps.push_back(kappa_precision(kappas, precisions, baserate, rng));
    return kps;
}

// prset() / pr.ct(): the full set's contingency table (counts). R's round()
// rounds halves to even.
QESTATS_R_ARITHMETIC_FN
inline ConfusionMatrix contingency_table(double precision, double recall,
                                         int length, double baserate) {
    QESTATS_R_ARITHMETIC
    const double gold1s = std::max(std::nearbyint(baserate * length), 1.0);
    const double gold0s = length - gold1s;
    const double tp = std::max(std::nearbyint(gold1s * recall), 1.0);
    const double fp = std::min(gold0s, std::max(std::nearbyint(tp * (1 - precision) / precision), 1.0));
    return {tp, fp, gold1s - tp, gold0s - fp};
}

// getHandSet(): kappa of a test set of `length` rows sampled from the full
// set laid out as prset() does (TP rows, FN rows, FP rows, TN rows). With an
// inflated base rate, ceiling(length * inflation) rows are first sampled from
// the first rater's positives.
template <class Rng>
QESTATS_R_ARITHMETIC_FN
double hand_kappa(const ConfusionMatrix& full, int length, double inflation, Rng& rng) {
    QESTATS_R_ARITHMETIC
    const auto tp = static_cast<std::size_t>(full.tp), fn = static_cast<std::size_t>(full.fn);
    const auto fp = static_cast<std::size_t>(full.fp), tn = static_cast<std::size_t>(full.tn);
    const std::size_t gold1s = tp + fn, rows = gold1s + fp + tn;
    double counts[4] = {0, 0, 0, 0};   // tp, fp, fn, tn
    auto add = [&](std::size_t row) {
        if (row < tp)                counts[0] += 1;
        else if (row < gold1s)       counts[2] += 1;
        else if (row < gold1s + fp)  counts[1] += 1;
        else                         counts[3] += 1;
    };
    const double positives = std::ceil(length * inflation);
    if (positives > static_cast<double>(gold1s))
        throw std::invalid_argument("Not enough positives in first rater to inflate to this level");
    if (positives > 0) {
        const auto k = static_cast<std::size_t>(positives);
        std::vector<bool> taken(rows, false);
        for (std::size_t row : sample_int(gold1s, k, rng)) { taken[row] = true; add(row); }
        std::vector<std::size_t> others;
        others.reserve(rows - k);
        for (std::size_t row = 0; row < rows; ++row) if (!taken[row]) others.push_back(row);
        for (std::size_t i : sample_int(others.size(), static_cast<std::size_t>(length) - k, rng))
            add(others[i]);
    } else {
        for (std::size_t row : sample_int(rows, static_cast<std::size_t>(length), rng)) add(row);
    }
    return kappa_from_counts(counts[0], counts[1], counts[2], counts[3]);
}

// mean(x, na.rm = TRUE): R sums in long double and adds a correction pass.
inline double r_mean(const std::vector<double>& x) {
    long double s = 0;
    std::size_t n = 0;
    for (double v : x) if (!std::isnan(v)) { s += v; ++n; }
    s /= static_cast<long double>(n);
    if (std::isfinite(static_cast<double>(s))) {
        long double t = 0;
        for (double v : x) if (!std::isnan(v)) t += (v - s);
        s += t / static_cast<long double>(n);
    }
    return static_cast<double>(s);
}

// getBootPvalue(): 1 below the distribution's mean, else the share of the
// distribution at or above `observed`.
inline double boot_pvalue(const std::vector<double>& distribution, double observed) {
    if (observed < r_mean(distribution)) return 1;
    std::size_t at_or_above = 0;
    for (double v : distribution) if (v >= observed) ++at_or_above;
    return static_cast<double>(at_or_above) / static_cast<double>(distribution.size());
}

// calcRho() with given (kappa, precision) pairs.
template <class Rng>
double rho(double observed_kappa, double baserate, int test_set_length,
           const std::vector<KappaPrecision>& kps, const RhoOptions& opts, Rng& rng) {
    const std::size_t replicates =
        std::min(kps.size(), static_cast<std::size_t>(std::max(opts.replicates, 0)));
    std::vector<double> kappas(replicates, 0.0);
    for (std::size_t i = 0; i < replicates; ++i) {
        const double r = recall(kps[i].kappa, baserate, kps[i].precision);
        const auto full = contingency_table(kps[i].precision, r, opts.ocs_length, baserate);
        kappas[i] = hand_kappa(full, test_set_length, opts.test_set_baserate_inflation, rng);
    }
    return boot_pvalue(kappas, observed_kappa);
}

} // namespace rho_standard

// ── Fast method: rhoR's C++ implementation ───────────────────────────────────
//
// Compiled as plain C++, like rhoR's src/rhor.cpp, so the compiler treats its
// arithmetic the way it treated rhoR's.

namespace rho_fast {

// recall(): the same formula as rho_standard::recall.
inline double recall(double kappa, double baserate, double precision) {
    const double top = kappa * precision;
    return top / (2 * precision - 2 * baserate - kappa + 2 * baserate * kappa);
}

// check_BRK_combo().
inline bool valid_combo(double baserate, double precision, double kappa) {
    const double right = ((2 * baserate * kappa) - (2 * baserate) - kappa) / (kappa - 2);
    return precision > right;
}

// Armadillo's regspace(start, delta, end).
inline std::vector<double> regspace(double start, double delta, double end) {
    if (((start < end) && (delta < 0)) || ((start > end) && (delta > 0)) || (delta == 0)) return {};
    const bool ascend = (start <= end);
    const double inc = (delta < 0) ? -delta : delta;
    const double m = (ascend ? (end - start) : (start - end)) / inc;
    const auto n = static_cast<std::size_t>(1 + static_cast<std::size_t>(std::floor(m)));
    // Two loops, as Armadillo writes them: a compiler that fuses
    // start + i * inc into an FMA then does it the same way here.
    std::vector<double> v(n);
    double* mem = v.data();
    if (ascend) {
        for (std::size_t i = 0; i < n; ++i) { mem[i] = double(start + double(i * inc)); }
    } else {
        for (std::size_t i = 0; i < n; ++i) { mem[i] = double(start - double(i * inc)); }
    }
    return v;
}

// One index from RcppArmadillo's sample(x, 1, FALSE): int(n * unif_rand()).
template <class Rng>
std::size_t sample_one(std::size_t n, Rng& rng) {
    const int n_orig = static_cast<int>(n);
    const int j = static_cast<int>(n_orig * rng.unif_rand());
    return static_cast<std::size_t>(j);
}

// Armadillo's randi<int>(distr_param(a, b)) on R's RNG (RcppArmadillo's
// Alt_R_RNG.h): runif(0, RAND_MAX) scaled to b - a + 1 values.
template <class Rng>
int randi(int a, int b, Rng& rng) {
    if (!(a <= b))
        throw std::logic_error("randi(): incorrect distribution parameters; a must be less than b");
    const double rand_max = static_cast<double>(RAND_MAX);
    const double draw = 0.0 + (rand_max - 0.0) * rng.unif_rand();
    if (a == 0 && b == RAND_MAX) return static_cast<int>(draw);
    const auto length = static_cast<unsigned long long>(b - a + 1);
    const double scale = double(length) / rand_max;
    return (std::min)(b, int(draw * scale) + a);
}

// find_valid_pk().
template <class Rng>
KappaPrecision kappa_precision(const std::vector<double>& kappas,
                               const std::vector<double>& precisions,
                               double baserate, Rng& rng) {
    for (int attempt = 0; attempt < RHO_MAX_ATTEMPTS; ++attempt) {
        const double kappa = kappas[sample_one(kappas.size(), rng)];
        double precision = precisions[sample_one(precisions.size(), rng)];
        if (!valid_combo(baserate, precision, kappa)) {
            const double precision_min =
                (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
            std::vector<double> above;
            for (double p : precisions) if (p > precision_min) above.push_back(p);
            if (above.empty()) continue;
            precision = above[sample_one(above.size(), rng)];
        }
        return {kappa, precision};
    }
    throw std::runtime_error(RHO_NO_VALID_SET);
}

// generate_kp_list() with the flat distribution (distributionType 0); rhoR's
// bell-shaped option (distributionType 1) is not ported.
template <class Rng>
std::vector<KappaPrecision>
generate_kps(int n, double baserate, double kappa_min, double kappa_max,
             double precision_min, double precision_max, Rng& rng,
             long distribution_length = 10000) {
    const double kappa_step = (kappa_max - kappa_min) / (distribution_length - 1);
    const auto kappas = regspace(kappa_min, kappa_step, kappa_max);
    const double precision_step = (precision_max - precision_min) / (10000 - 1);
    const auto precisions = regspace(precision_min, precision_step, precision_max);
    if (kappas.empty() || precisions.empty())
        throw std::invalid_argument("generate_kps: empty kappa or precision range");
    std::vector<KappaPrecision> kps;
    kps.reserve(n > 0 ? n : 0);
    for (int i = 0; i < n; ++i) kps.push_back(kappa_precision(kappas, precisions, baserate, rng));
    return kps;
}

// contingency_table(): the full set's table (counts). C's round() rounds
// halves away from zero.
inline ConfusionMatrix contingency_table(double precision, double recall,
                                         int length, double baserate) {
    const int gold1s = static_cast<int>(std::max(std::round(baserate * length), 1.0));
    const int gold0s = length - gold1s;
    const int tp = static_cast<int>(std::max(std::round(gold1s * recall), 1.0));
    const int fp = static_cast<int>(std::min(static_cast<double>(gold0s),
                                             std::max(std::round(tp * (1 - precision) / precision), 1.0)));
    return {double(tp), double(fp), double(gold1s - tp), double(gold0s - fp)};
}

// sample_contingency_table(): draw n items without replacement from cells
// (a 1x2 row or a 2x2 table in column order: TP, FP, FN, TN), returning each
// draw's 0-based cell. Kept as rhoR has it: the draw spans sum + 1 values and
// the first cell takes the extra one, so an empty first cell can be drawn.
// The draw spans every cell given; draws past the third cell land in the
// fourth.
template <class Rng>
std::vector<int> sample_contingency_table(std::vector<int> cells, int n, Rng& rng) {
    if (cells.size() < 4) cells.resize(4, 0);
    std::vector<int> out(n > 0 ? n : 0);
    for (int w = 0; w < n; ++w) {
        int max_int = 0;
        for (int c : cells) max_int += c;
        const int s = randi(0, max_int, rng);
        const int cell = (s <= cells[0]) ? 0
                       : (s <= cells[0] + cells[1]) ? 1
                       : (s <= cells[0] + cells[1] + cells[2]) ? 2 : 3;
        out[w] = cell;
        cells[cell]--;
    }
    return out;
}

// getHand_ct(): the test set's table, drawn from the full set's table.
template <class Rng>
ConfusionMatrix hand_contingency_table(const ConfusionMatrix& full, int length,
                                       double inflation, Rng& rng) {
    const int positives = static_cast<int>(std::ceil(length * inflation));
    int tp = static_cast<int>(full.tp), fp = static_cast<int>(full.fp);
    int fn = static_cast<int>(full.fn), tn = static_cast<int>(full.tn);
    int counts[4] = {0, 0, 0, 0};   // column order: tp, fp, fn, tn
    if (positives > 0) {
        int pos_tp = 0, pos_fn = 0;
        for (int cell : sample_contingency_table({tp, fn}, positives, rng))
            (cell == 0 ? pos_tp : pos_fn)++;
        tp -= pos_tp;
        fn -= pos_fn;
        counts[0] += pos_tp;
        counts[2] += pos_fn;
    }
    for (int cell : sample_contingency_table({tp, fp, fn, tn}, length - positives, rng))
        counts[cell]++;
    return {double(counts[0]), double(counts[1]), double(counts[2]), double(counts[3])};
}

// kappa_ct(): kappa as rhoR's C++ computes it. Within an ulp of
// kappa_from_counts, but NaN (not 1) when one cell holds every item.
inline double kappa_ct(const ConfusionMatrix& ct) {
    const double a = ct.tp, c = ct.fn, b = ct.fp, d = ct.tn;
    const double size = a + b + c + d;
    const double p_zero = (a + d) / size;
    const double p_yes = ((a + b) / size) * ((a + c) / size);
    const double p_no = ((c + d) / size) * ((b + d) / size);
    const double p_e = (p_yes + p_no);
    return (p_zero - p_e) / (1 - p_e);
}

// Armadillo's mean(): two interleaved accumulators.
inline double arma_mean(const std::vector<double>& x) {
    const std::size_t n = x.size();
    if (n == 0) throw std::logic_error("mean(): object has no elements");
    double acc1 = 0, acc2 = 0;
    std::size_t j;
    for (j = 1; j < n; j += 2) { acc1 += x[j - 1]; acc2 += x[j]; }
    if ((j - 1) < n) acc1 += x[j - 1];
    const double mean = (acc1 + acc2) / double(n);
    if (std::isfinite(mean)) return mean;
    for (double v : x) if (!std::isfinite(v)) return mean;
    double r_mean = 0;
    for (std::size_t i = 0; i < n; ++i) r_mean = r_mean + (x[i] - r_mean) / double(i + 1);
    return r_mean;
}

// getBootPvalue_c().
inline double boot_pvalue(const std::vector<double>& distribution, double observed) {
    if (observed < arma_mean(distribution)) return 1.0;
    std::size_t at_or_above = 0;
    for (double v : distribution) if (v >= observed) ++at_or_above;
    return (at_or_above * 1.0) / distribution.size();
}

// calcRho_c() with given (kappa, precision) pairs.
template <class Rng>
double rho(double observed_kappa, double baserate, int test_set_length,
           const std::vector<KappaPrecision>& kps, const RhoOptions& opts, Rng& rng) {
    const std::size_t replicates =
        std::min(kps.size(), static_cast<std::size_t>(std::max(opts.replicates, 0)));
    std::vector<double> kappas(replicates);
    for (std::size_t i = 0; i < replicates; ++i) {
        const double r = recall(kps[i].kappa, baserate, kps[i].precision);
        const auto full = contingency_table(kps[i].precision, r, opts.ocs_length, baserate);
        kappas[i] = kappa_ct(hand_contingency_table(full, test_set_length,
                                                    opts.test_set_baserate_inflation, rng));
    }
    return boot_pvalue(kappas, observed_kappa);
}

} // namespace rho_fast

// ── rho ──────────────────────────────────────────────────────────────────────

// rhoK()'s argument checks, with its messages.
inline void check_rho_args(double observed_kappa, double baserate, const RhoOptions& o) {
    if ((observed_kappa > 1 || observed_kappa < 0) || (baserate > 1 || baserate < 0))
        throw std::invalid_argument("x and OcSBaserate must be between 0 and 1.");
    if (o.scs_kappa_threshold < 0 || o.scs_kappa_min < 0)
        throw std::invalid_argument("OcSBaserate, ScSKappaThreshold, and ScSKappaMin must be positive.");
    if (o.scs_kappa_threshold < o.scs_kappa_min)
        throw std::invalid_argument("ScSKappaThreshold must be greater than ScSKappaMin.");
    if ((o.scs_precision_min < 0 || o.scs_precision_min > 1) ||
        (o.scs_precision_max < 0 || o.scs_precision_max > 1))
        throw std::invalid_argument("ScSPrecisionMin and ScSPrecisionMax must be between 0 and 1.");
    if (o.scs_precision_max < o.scs_precision_min)
        throw std::invalid_argument("ScSPrecisionMax must be greater than ScSPrecisionMin.");
}

// The (kappa, precision) pairs rho() simulates: opts.replicates draws with
// kappa in [scs_kappa_min, scs_kappa_threshold] and precision in
// [scs_precision_min, scs_precision_max].
template <class Rng>
std::vector<KappaPrecision> rho_kps(double baserate, const RhoOptions& opts, Rng& rng) {
    if (opts.method == RhoMethod::Fast)
        return rho_fast::generate_kps(opts.replicates, baserate, opts.scs_kappa_min,
                                      opts.scs_kappa_threshold, opts.scs_precision_min,
                                      opts.scs_precision_max, rng);
    return rho_standard::generate_kps(opts.replicates, baserate, opts.scs_kappa_min,
                                      opts.scs_kappa_threshold, opts.scs_precision_min,
                                      opts.scs_precision_max, rng);
}

// rhoK(): rho for an observed kappa on a test set of `test_set_length` items
// from a population with this base rate.
template <class Rng>
double rho(double observed_kappa, double baserate, int test_set_length,
           const RhoOptions& opts, Rng& rng) {
    check_rho_args(observed_kappa, baserate, opts);
    const auto kps = rho_kps(baserate, opts, rng);
    if (opts.method == RhoMethod::Fast)
        return rho_fast::rho(observed_kappa, baserate, test_set_length, kps, opts, rng);
    return rho_standard::rho(observed_kappa, baserate, test_set_length, kps, opts, rng);
}

inline double rho(double observed_kappa, double baserate, int test_set_length,
                  const RhoOptions& opts = {}, std::uint64_t seed = 5489u) {
    StdRng rng(seed);
    return rho(observed_kappa, baserate, test_set_length, opts, rng);
}

} // namespace qestats

// rhoR's C++ functions, on libqe-stats: rho and its building blocks are
// qestats::rho_fast (this file's functions, as they were) and
// qestats::rho_standard (rhoK's default method), vendored in
// inst/include/qestats/ by tools/sync-qestats.sh. Given R's RNG they return
// exactly what rhoR returned before the move (tests/testthat/test-regression-pins.R).
//
// Only the bell-shaped kappa distribution (generate_kp_list(distributionType = 1))
// is still computed here, with RcppArmadillo's weighted sample().

#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
#include <Rcpp.h>
#include "qestats/rho.hpp"
using namespace Rcpp;

namespace {

// R's RNG for qestats: unif_rand() for the fast method, R_unif_index() (what
// sample.int() draws with) for the standard one.
struct RRng {
  double unif_rand() { return ::unif_rand(); }
  double unif_index(double n) { return R_unif_index(n); }
};

// rhoR's tables are gold (rows) by silver (columns): [1,1] TP, [1,2] FN,
// [2,1] FP, [2,2] TN.
qestats::ConfusionMatrix to_cm(const IntegerMatrix& ct) {
  return {double(ct(0, 0)), double(ct(1, 0)), double(ct(0, 1)), double(ct(1, 1))};
}

IntegerMatrix to_ct(const qestats::ConfusionMatrix& cm) {
  IntegerMatrix ct(2, 2);
  ct(0, 0) = cm.tp;
  ct(0, 1) = cm.fn;
  ct(1, 0) = cm.fp;
  ct(1, 1) = cm.tn;
  return ct;
}

std::vector<qestats::KappaPrecision> to_kps(const NumericMatrix& KPs) {
  std::vector<qestats::KappaPrecision> kps(KPs.nrow());
  for (int i = 0; i < KPs.nrow(); i++) kps[i] = {KPs(i, 1), KPs(i, 0)};
  return kps;
}

NumericMatrix from_kps(const std::vector<qestats::KappaPrecision>& kps) {
  NumericMatrix KPs(kps.size(), 2);
  for (std::size_t i = 0; i < kps.size(); i++) {
    KPs(i, 0) = kps[i].precision;
    KPs(i, 1) = kps[i].kappa;
  }
  return KPs;
}

// find_valid_pk() with kappa probabilities (the bell-shaped distribution).
// Precision is drawn unweighted, as before.
NumericVector find_valid_pk_weighted(
  arma::colvec kappaDistribution, arma::vec kappaProbability,
  arma::colvec precisionDistribution, double baserate
) {
  for (int attempt = 0; attempt < qestats::RHO_MAX_ATTEMPTS; attempt++) {
    arma::ivec kappaRng = arma::regspace<arma::ivec>(0, kappaDistribution.size() - 1);
    arma::ivec whKappa = Rcpp::RcppArmadillo::sample(kappaRng, 1, false, kappaProbability);
    double currKappa = kappaDistribution[whKappa.at(0)];

    arma::ivec precRng = arma::regspace<arma::ivec>(0, precisionDistribution.size() - 1);
    arma::ivec whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
    double currPrec = precisionDistribution[whPrec.at(0)];

    if (!qestats::rho_fast::valid_combo(baserate, currPrec, currKappa)) {
      double precisionMin = (2 * baserate * currKappa - 2 * baserate - currKappa) / (currKappa - 2);
      arma::uvec indicies = arma::find(precisionDistribution > precisionMin);
      if (indicies.size() == 0) continue;
      arma::colvec above = precisionDistribution.elem(indicies);
      precRng = arma::regspace<arma::ivec>(0, above.size() - 1);
      whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
      currPrec = above[whPrec.at(0)];
    }
    return NumericVector::create(currPrec, currKappa);
  }
  stop(qestats::RHO_NO_VALID_SET);
}

} // namespace

//' @title sample_contingency_table
//' @name sample_contingency_table
//'
//' @param xx contingency table matrix
//' @param n int size of the contingency table
//' @param forR bool if true, add 1 to the results accounting for R indices starting at 1
//'
//' @export
// [[Rcpp::export]]
IntegerMatrix sample_contingency_table(IntegerMatrix xx, int n, bool forR = true) {
  RRng rng;
  std::vector<int> cells = qestats::rho_fast::sample_contingency_table(
    std::vector<int>(xx.begin(), xx.end()), n, rng);
  IntegerMatrix ret(cells.size(), 1);
  for (std::size_t w = 0; w < cells.size(); w++) ret[w] = cells[w] + (forR ? 1 : 0);
  return ret;
}

//' @name getBootPvalue_c
//' @title getBootPvalue_c
//'
//' @param distribution vector of calculated kappas
//' @param result double calculated kappa to compare against
//'
//' @description
//' returns the percentage of the time that the distribution was greater or equal to the observed kappa
//' if the result is less than the mean of the distribution, than the p value is 1
//' else return the number of times that the distribution is greater than the result as a percentage of the total number
//' of items in the distribution
//'
//' @return double calculated p-value
//'
//' @export
// [[Rcpp::export]]
double getBootPvalue_c(NumericVector distribution, double result) {
  return qestats::rho_fast::boot_pvalue(
    std::vector<double>(distribution.begin(), distribution.end()), result);
}

// Validate a Baserate/Kappa combo against the supplied Precision.
// Return true if Precision is greater than the calculated right value
// [[Rcpp::export]]
bool check_BRK_combo(double BR, double P, double K) {
  return qestats::rho_fast::valid_combo(BR, P, K);
}

//' @name recall
//' @title recall
//'
//' @param kappa double
//' @param BR double
//' @param P double
//'
//' @return Recall calculated from provided kappa, BR, and P
//'
//' @export
// [[Rcpp::export]]
double recall(double kappa, double BR, double P) {
  return qestats::rho_fast::recall(kappa, BR, P);
}

// Find a valid precision and kappa combo from a distribution of Kappas and Precisions given a baserate
//
// [[Rcpp::export]]
NumericVector find_valid_pk(
  NumericVector kappaDistribution, NumericVector kappaProbability,
  NumericVector precisionDistribution, NumericVector precisionProbability,
  double baserate
) {
  if (kappaDistribution.size() == 0 || precisionDistribution.size() == 0) {
    stop("find_valid_pk: empty kappa or precision distribution");
  }
  if (kappaProbability.size() > 0) {
    return find_valid_pk_weighted(as<arma::colvec>(kappaDistribution), as<arma::vec>(kappaProbability),
                                  as<arma::colvec>(precisionDistribution), baserate);
  }
  RRng rng;
  qestats::KappaPrecision kp = qestats::rho_fast::kappa_precision(
    std::vector<double>(kappaDistribution.begin(), kappaDistribution.end()),
    std::vector<double>(precisionDistribution.begin(), precisionDistribution.end()),
    baserate, rng);
  return NumericVector::create(kp.precision, kp.kappa);
}


//' @name generateKPs_c
//' @title generate_kp_list
//'
//' @param numNeeded int
//' @param baserate double
//' @param kappaMin double
//' @param kappaMax double
//' @param precisionMin double
//' @param precisionMax double
//' @param distributionType int 0 - normal (default), 1 - bell
//' @param distributionLength long
//'
//' @return matrix of kappa and precision values (column 1 as precision)
//'
//' @export
// [[Rcpp::export]]
Rcpp::NumericMatrix generate_kp_list(
    int numNeeded, double baserate,
    double kappaMin, double kappaMax,
    double precisionMin, double precisionMax,
    int distributionType = 0, long distributionLength = 10000
) {
  if (distributionType != 1) {
    RRng rng;
    return from_kps(qestats::rho_fast::generate_kps(
      numNeeded, baserate, kappaMin, kappaMax, precisionMin, precisionMax, rng, distributionLength));
  }

  double kappaStep = ((kappaMax - kappaMin) / (distributionLength - 1));
  arma::colvec kappaDistribution = arma::regspace(kappaMin, kappaStep, kappaMax);
  arma::vec kappaProbability = arma::normpdf(kappaDistribution, 0.9, 0.1);
  double precStep = (precisionMax - precisionMin) / (10000 - 1);
  arma::colvec precisionDistribution = arma::regspace(precisionMin, precStep, precisionMax);

  Rcpp::NumericMatrix KPs(numNeeded, 2);
  for (int i = 0; i < numNeeded; i++) {
    KPs(i, _) = find_valid_pk_weighted(kappaDistribution, kappaProbability, precisionDistribution, baserate);
  }
  return KPs;
}

//' @name contingency_table
//' @title contingency_table
//' @description Create a contingency table using the provied precision, recall, baserate, and length.
//'
//' @param precision double
//' @param rec double
//' @param length int
//' @param baserate double
//'
//' @export
// [[Rcpp::export]]
IntegerMatrix contingency_table(double precision, double rec, int length, double baserate) {
  return to_ct(qestats::rho_fast::contingency_table(precision, rec, length, baserate));
}

//' @name random_contingency_table
//' @title random_contingency_table
//' @param setLength [TBD]
//' @param baserate [TBD]
//' @param kappaMin [TBD]
//' @param kappaMax [TBD]
//' @param minPrecision [TBD]
//' @param maxPrecision [TBD]
//'
//' @export
// [[Rcpp::export]]
IntegerMatrix random_contingency_table(
    int setLength,double baserate,
    double kappaMin, double kappaMax,
    double minPrecision = 0, double maxPrecision = 1
) {
  RRng rng;
  qestats::KappaPrecision kp = qestats::rho_fast::generate_kps(
    1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision, rng)[0];
  double rec = qestats::rho_fast::recall(kp.kappa, baserate, kp.precision);
  return to_ct(qestats::rho_fast::contingency_table(kp.precision, rec, setLength, baserate));
}

//' @title kappa_ct
//' @description Calculate kappa from a contingency table
//' @param ct [TBD]
//'
//' @export
// [[Rcpp::export]]
double kappa_ct(IntegerMatrix ct) {
  return qestats::rho_fast::kappa_ct(to_cm(ct));
}

// [[Rcpp::export]]
IntegerMatrix getHand_ct(IntegerMatrix ct, int handSetLength, double handSetBaserate) {
  RRng rng;
  return to_ct(qestats::rho_fast::hand_contingency_table(to_cm(ct), handSetLength, handSetBaserate, rng));
}

//' @name getHand_kappa
//' @title getHand_kappa
//'
//' @description This function returns kappa calculated from a Handset taken from  a larger Contingency Table
//'
//' @param ct KPs matrix of kappa (column 1) and precision (column 2) values
//' @param handSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
//' @param handSetBaserate baserate to inflate the sampled contingency table to
//'
//' @return Kappa as double
//'
//' @export
// [[Rcpp::export]]
double getHand_kappa(IntegerMatrix ct, int handSetLength, double handSetBaserate) {
  RRng rng;
  return qestats::rho_fast::kappa_ct(
    qestats::rho_fast::hand_contingency_table(to_cm(ct), handSetLength, handSetBaserate, rng));
}

// This function calculates rho given the parameters of the initial set and
// the kappa of the observed handset and is called by rhoK()
//
// [[Rcpp::export]]
double calcRho_c(
  double x,
  double OcSBaserate,
  int testSetLength,
  double testSetBaserateInflation = 0,
  int OcSLength = 10000,
  int replicates = 800,
  double ScSKappaThreshold = 0.9,
  double ScSKappaMin = 0.40,
  double ScSPrecisionMin = 0.6,
  double ScSPrecisionMax = 1.0,
  NumericMatrix KPs = NumericMatrix(0)
) {
  RRng rng;
  qestats::RhoOptions opts;
  opts.test_set_baserate_inflation = testSetBaserateInflation;
  opts.ocs_length = OcSLength;
  opts.replicates = replicates;
  std::vector<qestats::KappaPrecision> kps;
  if (KPs.size() == 1 && KPs[0] == 0) {
    kps = qestats::rho_fast::generate_kps(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold,
                                          ScSPrecisionMin, ScSPrecisionMax, rng);
  } else {
    kps = to_kps(KPs);
  }
  return qestats::rho_fast::rho(x, OcSBaserate, testSetLength, kps, opts, rng);
}

// calcRho()'s computation (rhoK's default method): qestats::rho_standard.
// KPs, when given, is a matrix of precision (column 1) and kappa (column 2).
// [[Rcpp::export]]
double calcRho_standard_c(
  double x,
  double OcSBaserate,
  int testSetLength,
  double testSetBaserateInflation,
  int OcSLength,
  int replicates,
  double ScSKappaThreshold,
  double ScSKappaMin,
  double ScSPrecisionMin,
  double ScSPrecisionMax,
  Nullable<NumericMatrix> KPs = R_NilValue
) {
  RRng rng;
  qestats::RhoOptions opts;
  opts.test_set_baserate_inflation = testSetBaserateInflation;
  opts.ocs_length = OcSLength;
  opts.replicates = replicates;
  std::vector<qestats::KappaPrecision> kps;
  if (KPs.isNull()) {
    kps = qestats::rho_standard::generate_kps(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold,
                                              ScSPrecisionMin, ScSPrecisionMax, rng);
  } else {
    kps = to_kps(NumericMatrix(KPs.get()));
  }
  return qestats::rho_standard::rho(x, OcSBaserate, testSetLength, kps, opts, rng);
}

1		###
2		# Precision Recall Set
3		#
4		# This creates a set with a given recall and precision
5		# @param precision This is the precision of the given set
6		# @param recall This is the recall of the given set
7		# @param length This is the length of the set
8		# @param baserate This is the baserate of the set
9		# @return returns a set with the given parameters
10		###
11		prset = function(precision, recall, length, baserate) {
12		#generates a set given the different 4 quadrants
13	59x	gold1s = max(round(baserate * length , 0), 1); # TP + FN
14	59x	gold0s = length - gold1s; # FP + FN
15	59x	TP = max(round(gold1s * recall , 0) , 1); # TP
16	59x	FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
17	59x	gold = c(rep(1, gold1s), rep(0, gold0s));
18	59x	silver = c(rep(1, TP),rep(0, gold1s - TP), rep(1, FP), rep(0, gold0s - FP));
19
20	59x	return(cbind(gold, silver));
21		}
22
23		pr.ct = function(precision, recall, length, baserate) {
24		#generates a set given the different 4 quadrants
25	4x	gold1s = max(round(baserate * length , 0), 1); # TP + FN
26	4x	gold0s = length - gold1s; # FP + FN
27	4x	TP = max(round(gold1s * recall , 0) , 1); # TP
28	4x	FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
29
30	4x	ct = matrix(0, ncol = 2, nrow = 2)
31	4x	ct[1,1] = TP
32	4x	ct[1,2] = gold1s - TP
33	4x	ct[2,1] = FP
34	4x	ct[2,2] = gold0s - FP
35	4x	return(ct);
36		}

1		any_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
2	54x	any(sapply(target, function(targ) abs(current - targ) < tolerance ))
3		}
4		all_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
5	5x	all(sapply(target, function(targ) abs(current - targ) < tolerance ))
6		}
7
8
9		#' Convert a codeset to a contingency table
10		#'
11		#' @param x codeset
12		#'
13		#' @return contingency table as a 2x2 matrix
14		#'
15		#' @export
16		as.contingency.table <- function(x) {
17	31x	y <- NULL
18
19	31x	if ( nrow(x) > 2 && ncol(x) == 2 ) {
20	18x	tp = length(which(x[,1] == 1 & x[,2] == 1))
21	18x	tn = length(which(x[,1] == 0 & x[,2] == 0))
22	18x	fp = length(which(x[,1] == 0 & x[,2] == 1))
23	18x	fn = length(which(x[,1] == 1 & x[,2] == 0))
24
25	18x	y <- matrix(
26	18x	c(tp, fp, fn, tn),
27	18x	ncol = 2, nrow = 2
28		)
29		}
30		else {
31	13x	y <- x
32		}
33
34	31x	dimnames(y) = list(c("first_positive", "first_negative"), c("second_positive", "second_negative"))
35	31x	class(y) <- c("contingency.table", "rating.set", class(y))
36	31x	attr(y, "agreement") <- list(
37	31x	"true_positive" = y[1, 1],
38	31x	"false_negative" = y[1, 2],
39	31x	"false_positive" = y[2, 1],
40	31x	"true_negative" = y[2, 2]
41		)
42	31x	attr(y, "baserate") <- baserateCT(y)
43	31x	attr(y, "kappa") <- kappaCT(y)
44
45	31x	return(y)
46		}
47
48		#' Convert codeset to contingency table
49		#'
50		#' @param x matrix contingency table (2x2)
51		#'
52		#' @return 2-column matrix representation of the contingency table
53		#'
54		#' @export
55		as.code.set <- function(x) {
56	8x	y <- NULL
57	8x	if (all(dim(x) == c(2,2))) {
58	3x	y <- matrix(
59	3x	c(rep(1, x[1,1]2), rep(1:0, x[1,2]), rep(0:1, x[2,1]), rep(0, x[2,2]2)),
60	3x	byrow = TRUE,
61	3x	ncol = 2
62		)
63		}
64		else {
65	5x	y <- x
66		}
67
68	8x	class(y) <- c("code.set", "rating.set", class(y))
69	8x	attr(y, "baserate") <- baserateSet(y)
70	8x	attr(y, "kappa") <- kappaSet(y)
71	8x	colnames(y) <- c("first_rater", "second_rater")
72
73	8x	return(y)
74		}
75
76		#' Helper function to return special values on a rating set
77		#'
78		#' @param x Set or Contingency.Table
79		#'
80		#' @param i Value to search for
81		#'
82		#' @export
83		"$.rating.set" <- function (x, i) {
84	91x	if (i %in% names(attributes(x))) {
85	86x	return(attr(x, i))
86		}
87		else {
88	5x	if (is(x, "code.set")) {
89	2x	if (i %in% colnames(x)) {
90	1x	return(x[, which(colnames(x) == i)])
91		}
92		}
93	3x	else if (is(x, "contingency.table")) {
94	3x	if (i %in% colnames(x)) {
95	1x	return(sum(x[, which(colnames(x) == i)]))
96		}
97	2x	else if (i %in% rownames(x)) {
98	1x	return(sum(x[which(rownames(x) == i), ]))
99		}
100		}
101		}
102
103	2x	return(NULL)
104		}
105
106		#' @export
107		.DollarNames.rating.set <- function(x, pattern) {
108	2x	vals <- .codegen_dollar(x, pattern)
109
110	2x	cols <- colnames(x)
111	2x	if (!is.null(cols)) {
112	2x	vals <- c(cols, vals)
113		}
114
115	2x	rows <- rownames(x)
116	2x	if (!is.null(rows)) {
117	1x	vals <- c(rows, vals)
118		}
119
120	2x	vals
121		}
122
123		.codegen_dollar <- function(x, pattern = "") {
124	3x	attrs <- names(attributes(x))
125	3x	vals <- attrs[attrs %in% c("baserate", "kappa", "agreement")]
126	3x	return(vals)
127		}
128
129		#' @export
130		print.rating.set <- function(x, ...) {
131	1x	y <- unclass(x)
132
133	1x	for(a in additional_attributes) {
134	3x	attr(y, a) <- NULL
135		}
136
137	1x	print(y, ...)
138		}

1		contingencyToSet = function(TP, FP, FN, TN){
2	48x	setLength = TP + FP + FN + TN;
3
4	48x	gold1s = TP + FN;
5	48x	gold0s = setLength - gold1s;
6
7	48x	gold = c(rep(1, gold1s), rep(0, gold0s));
8	48x	silver = c(rep(1, TP),rep(0, gold1s - TP), rep(1, FP), rep(0, gold0s - FP));
9
10	48x	return(cbind(gold, silver));
11		};

1		###
2		#' @title Create Simulated codeSet
3		#'
4		#' @description
5		#' Creates a simulated \code{\link{codeSet}} with the given parameters
6		#'
7		#' @param length the length of the simulated \code{\link{codeSet}} to be created
8		#' @param baserate the \code{\link{baserate}} of the simulated \code{\link{codeSet}}
9		#' @param kappaMin the minimum kappa of the simulated \code{\link{codeSet}}
10		#' @param kappaMax the maximum kappa of the simulated \code{\link{codeSet}}
11		#' @param precisionMin the minimum precision of the simulated \code{\link{codeSet}}
12		#' @param precisionMax the maximum precision of the simulated \code{\link{codeSet}}
13		#' @param ... Parameters passed to createRandomSet (e.g. type = "set" or type = "ct")
14		#' @param tries the maximum number of tries to generate a valid set, smaller set lengths may require an increased number of tries
15		#'
16		#' @details
17		#' \code{\link{codeSet}}s are generated by first picking a random kappa within its range and a random precision within its range. If the random kappa, random precision, and baserate are not mathematically possible, then the precision is resampled from a range of mathematically possible values within its range. A unique simulated \code{\link{codeSet}} is then constructed given these parameters.
18		#'
19		#' @export
20		#' @return A \code{\link{codeSet}} that fulfills the given parameters
21		###
22		createSimulatedCodeSet = function(length, baserate, kappaMin, kappaMax, precisionMin, precisionMax, ..., tries = 50){
23	1x	if(length <= 0){stop("length must be positive")}
24	1x	if(baserate < 0){stop("baserate must be positive")}
25	2x	if(kappaMin < 0 \| kappaMin > 1) stop("kappaMin must be between 0 and 1.")
26	2x	if(kappaMax < 0 \| kappaMax > 1) stop("kappaMax must be between 0 and 1.")
27	1x	if(kappaMin > kappaMax){stop("kappaMin must be less than kappaMax")}
28	2x	if(precisionMin < 0 \| precisionMin > 1) stop("precisionMin must be between 0 and 1.")
29	2x	if(precisionMax < 0 \| precisionMax > 1) stop("precisionMax must be between 0 and 1.")
30	1x	if(precisionMin > precisionMax){stop("precisionMin must be less than precisionMax")}
31	1x	if(tries < 1){stop("tries must be greater than 1")}
32
33	3x	setFound = NULL
34	3x	validSet = F
35	3x	triesLeft = tries
36	3x	Ks = c()
37	3x	while(!validSet && triesLeft > 0) {
38	52x	set = createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax)
39	52x	k = kappa(set)
40
41	52x	if ((k > kappaMin && k < kappaMax) \|\| any_equal(k, c(kappaMin, kappaMax))) {
42	2x	setFound = set
43	2x	validSet = T
44		}
45
46	52x	Ks = c(Ks, k)
47	52x	triesLeft = triesLeft - 1
48		}
49
50	3x	if (is.null(setFound)) {
51	1x	Ks_under_min <- Ks[which(Ks < kappaMin)]
52	1x	Kdist.to.min <- kappaMin - Ks_under_min
53	1x	Kdist.to.min.low <- ifelse(length(Ks_under_min) > 0, min(Kdist.to.min), NA)
54
55	1x	Ks_over_min <- Ks[which(Ks > kappaMax)]
56	1x	Kdist.to.max <- Ks_over_min - kappaMax
57	1x	Kdist.to.max.low <- ifelse(length(Ks_over_min) > 0, min(Kdist.to.max), NA)
58
59	1x	Kdist.all_mins <- c(Kdist.to.min, Kdist.to.max)
60
61	1x	both_mins <- c(Kdist.to.min.low, Kdist.to.max.low)
62	1x	both_mins <- both_mins[!is.na(both_mins)]
63	1x	closest <- c(Ks_under_min, Ks_over_min) [which(Kdist.all_mins == min(both_mins))][1]
64
65	1x	stop(paste0(
66	1x	"Unable to create a set of the provided length (", length, ") ",
67	1x	"in the provided kappa range (", kappaMin, "\U2014", kappaMax, "). ",
68	1x	"The closest kappa obtained after ", tries, " tries was: ",
69	1x	round(x = closest, digits = 5), "\n",
70	1x	"\tTry increasing the set length or the range of valid kappas"
71		))
72		}
73
74	2x	return(setFound)
75
76		# return(createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax, ...))
77		}

1		generateKPs = function(
2		numNeeded, baserate,
3		kappaMin, kappaMax,
4		precisionMin, precisionMax,
5		distributionType = 'FLAT',
6		distributionLength = 100000
7		) {
8	62x	kappaDistribution = seq(kappaMin, kappaMax, length = distributionLength)
9
10	62x	if (distributionType == 'FLAT') {
11		#probability for a flat distribution
12	60x	kappaProbability = NULL
13		}
14	2x	else if (distributionType == 'BELL') {
15		#probability for a bell curve
16	2x	kappaProbability = stats::dnorm(kappaDistribution, mean = 0.9, sd = 0.1)
17		}
18
19	62x	precisionDistribution = seq(precisionMin, precisionMax, length = 10000)
20	62x	precisionProbability = NULL
21
22	62x	KPs = tryCatch({
23	62x	KPs = list()
24	62x	for (i in 1:numNeeded) {
25	10082x	KPs[[length(KPs) + 1]] = genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate)
26		}
27	61x	KPs
28	62x	},error=function(e){
29	1x	NULL
30		})
31	1x	if(is.null(KPs)){stop("Could not generate a valid set given ranges of kappa and precision")}
32
33	61x	return (KPs);
34		}

1		###
2		#' @title Get Test Set
3		#'
4		#' @description
5		#' This function gets a \emph{testSet} from a larger \code{\link{codeSet}} given certain sampling parameters.
6		#'
7		#' @details
8		#' A \emph{testSet} is a \code{\link{codeSet}} that is a subset of a larger \code{\link{codeSet}} with a given set of properties. A \emph{testSet} is constructed by sampling (without replacement) P rows from rows in the larger \code{\link{codeSet}} where the first rater's code was 1, and then appending an additional sample (without replacement) of R rows taken at random from the larger \code{\link{codeSet}} excluding rows included in the first P rows sampled. P is computed as the minbaserate * length of the \emph{testset}. R is computed as testSetLength - P. The result of this sampling procedure is to create a sample with a minimum baserate regardless of the baserate of the larger \code{\link{codeSet}}.If \emph{testSetBaserateInflation} is set to zero, the function selects rows at random.
9		#'
10		#' @param set The \code{\link{codeSet}} from which the \emph{testSet} is taken
11		#' @param testSetLength The length of the \emph{testSet} to be taken
12		#' @param testSetBaserateInflation The minimum guaranteed \code{\link{baserate}} of the \emph{testSet}. Default to 0
13		#'
14		#'
15		#' @export
16		#' @return A \code{\link{codeSet}} with the properties specified
17		###
18
19		getTestSet = function(set, testSetLength, testSetBaserateInflation = 0){
20	1x	if(length(unlist(set)) < testSetLength){stop("testSetLength must be less than the length of set")}
21	1x	if(testSetLength %% 1 \| testSetLength < 1) {stop("testSetLength value must be a positive integer.")}
22	1x	if(testSetBaserateInflation < 0){stop("testSetBaserateInfation must be positive")}
23	1x	if(testSetBaserateInflation >= 1){stop("testSetBaserateInflation must be below 1")}
24	5x	if(length(which(set[,1] == 1 & set[,2] == 1)) == nrow(set)){
25	1x	return(set[1:testSetLength,])
26	4x	}else if(length(which(set[,1] == 0 & set[,2] == 0)) == nrow(set)){
27	2x	if(testSetBaserateInflation == 0){
28	1x	return(set[1:testSetLength,])
29		}else{
30	1x	stop("Not enough positives in first rater to inflate to this level")
31		}
32		}
33	2x	return(getHandSet(set = set, handSetLength = testSetLength, handSetBaserate = testSetBaserateInflation, returnSet = T))
34		}

1		###
2		#' @title Rho (kappa)
3		#'
4		#' @description
5		#' This function calculates rho for an observed kappa value with associated set parameters
6		#' (testSetLength and OcSBaserate). Called by \code{\link{rho}}. A p-value is returned and if
7		#' this value is less than 0.05, it is said that the handset does generalize to the entire set
8		#'
9		#' @template rho-params
10		#' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
11		#' @param method set to "c" to calculate using the C++ implmentation. Defaults to "standard"
12		#'
13		#' @export
14		#' @return rho for the given parameters
15		###
16		rhoK = function(
17		x,
18		OcSBaserate,
19		testSetLength,
20		testSetBaserateInflation = 0,
21		OcSLength = 10000,
22		replicates = 800,
23		ScSKappaThreshold = 0.9,
24		ScSKappaMin = 0.40,
25		ScSPrecisionMin = 0.6,
26		ScSPrecisionMax = 1,
27		method = "standard"
28		) {
29		if (
30	85x	(x > 1 \| x < 0) \|\|
31	85x	(OcSBaserate > 1 \| OcSBaserate < 0)
32	3x	) stop("x and OcSBaserate must be between 0 and 1.")
33
34	82x	if (OcSBaserate < 0 \|\| ScSKappaThreshold < 0 \|\| ScSKappaMin < 0)
35	2x	stop("OcSBaserate, ScSKappaThreshold, and ScSKappaMin must be positive.")
36
37		if (
38	80x	testSetLength %% 1 != 0 \|\|
39	80x	OcSLength %% 1 != 0 \|\|
40	80x	replicates %% 1 != 0
41	3x	) stop("testSetLength, OcSLength, and replicates values must be a positive integer.")
42
43	77x	if (ScSKappaThreshold < ScSKappaMin)
44	1x	stop("ScSKappaThreshold must be greater than ScSKappaMin.")
45
46		if (
47	76x	(ScSPrecisionMin < 0 \| ScSPrecisionMin > 1) \|\|
48	76x	(ScSPrecisionMax < 0 \| ScSPrecisionMax > 1)
49	4x	) stop("ScSPrecisionMin and ScSPrecisionMax must be between 0 and 1.")
50
51	72x	if (ScSPrecisionMax < ScSPrecisionMin)
52	1x	stop("ScSPrecisionMax must be greater than ScSPrecisionMin.")
53
54	71x	result <- NULL
55	71x	if (method == "standard") {
56	48x	result <- calcRho( x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
57		}
58		else {
59	23x	result <- calcRho_c(x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
60		}
61	71x	return(result)
62		}

1		###
2		#' @title Get Handset
3		#' @description This function is to get a handset of a set and calculate the kappa
4		#' @param set This is the set to take a handset of
5		#' @param handSetLength This is the length of the handset to take
6		#' @param handSetBaserate This is the minimum baserate to inflate the handset to
7		#' @param returnSet If TRUE, then return the handSet if FALSE, return the kappa of the handSet
8		#'
9		#' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
10		#' @export
11		###
12		getHandSet = function(set, handSetLength, handSetBaserate, returnSet=FALSE) {
13		#positives is the minimum number of positive pairs in the first rater
14	9x	positives = ceiling(handSetLength * handSetBaserate);
15	9x	posInd = which(set[,1] == 1);
16
17	1x	if(positives > length(posInd)){stop("Not enough positives in first rater to inflate to this level")}
18		#if there is an inflated baserate positives will be > 0, so if handSetBaserate > 0, positives > 0
19	8x	if (positives > 0) {
20		#goes through and picks a set that fits the number of positives first then randomly samples from the rest
21		#if the set created is not a valid set, then it reruns the process until a valid set is found
22	5x	positiveIndices = posInd[sample.int(length(posInd),size=positives,replace=FALSE)];
23	5x	others = set[!(1:nrow(set) %in% positiveIndices),];
24	5x	otherIndices = sample.int(nrow(others),size=(handSetLength - positives),replace=FALSE);
25	5x	this.set = rbind(set[positiveIndices,], others[otherIndices,]);
26	3x	} else if (positives == 0){
27		#if there is no restriction on positives, then the set is generated by randomly sampling from the entire set
28		#this.set=matrix(0,2,2);
29		#while ((sum(this.set[,2]) == 0 \| sum(this.set[,2]) == nrow(this.set))) {
30	3x	theseIndices = sample.int(nrow(set),size=handSetLength,replace=FALSE);
31	3x	this.set = set[theseIndices,];
32		#}
33		}
34
35		#if return set is true return set rather than kappa
36	8x	if (returnSet) {
37	5x	this.set <- as.code.set(this.set);
38	5x	return(this.set);
39		}
40
41		#otherwise return the kappa of the set
42	3x	handSetKappa = calcKappa(this.set);
43
44	3x	return(handSetKappa);
45		}
46
47		###
48		#' @title Get Handset
49		#' @description This function is to get a handset of a set and calculate the kappa
50		#'
51		#' @param full.ct This is the set to take a handset of
52		#' @param handSetLength This is the length of the handset to take
53		#' @param handSetBaserate This is the minimum baserate to inflate the handset to
54		#' @param as_kappa If FALSE then return the handSet, if TRUE (default) return the kappa of the handSet
55		#'
56		#' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
57		#' @export
58		###
59		getHandCT <- function(full.ct, handSetLength, handSetBaserate, as_kappa = TRUE) {
60	4x	positives <- ceiling(handSetLength * handSetBaserate)
61	4x	if (positives > sum(full.ct[1,])) {
62	1x	stop("Not enough positives in first rater to inflate to this level")
63		}
64
65		# positiveIndices <- sample_ct(full.ct[1,], positives)
66		#
67		# other_ct <- full.ct
68		# other_ct[1, ] <- other_ct[1, ] - positiveIndices
69		#
70		# other_ct_vec <- as.vector(other_ct)
71		# this_ct <- sample_ct(other_ct_vec, handSetLength - positives)
72		# this_ct[1, ] <- this_ct[1, ] + positiveIndices
73
74	3x	this_ct <- getHand_ct(ct = full.ct, handSetLength = handSetLength, handSetBaserate = handSetBaserate)
75
76	3x	if (!as_kappa) {
77	1x	this_ct <- as.contingency.table(this_ct);
78	1x	return(this_ct)
79		}
80
81	2x	return(kappa_ct(this_ct))
82		}
83
84		# sample_ct <- function(x, size) {
85		# x_vec <- as.vector(x)
86		# x_sampled <- matrix(rep(0, length(x)), ncol = length(x_vec) / 2)
87		# for(i in seq(size)) {
88		# s <- sample(seq(x_vec), 1, replace = F, prob = x_vec / sum(x_vec))
89		# x_sampled[[s]] = x_sampled[[s]] + 1
90		# x_vec[s] = x_vec[s] - 1
91		# }
92		#
93		# return(x_sampled)
94		# }

1		###
2		#' @title Rho (contingency Table)
3		#'
4		#' @description
5		#' This function calculates rho and kappa for a given \code{\link{contingencyTable}}, and returns a list containing both values. Called by \code{\link{rho}}.
6		#'
7		#'
8		#' @template rho-params
9		#'
10		#' @export
11		#' @return A list of the format:\describe{
12		#' \item{rho}{The rho of the \code{\link{contingencyTable}}}
13		#' \item{kappa}{The Cohen's Kappa of the \code{\link{contingencyTable}}}
14		#' }
15		###
16		rhoCT = function(
17		x,
18		OcSBaserate = NULL,
19		testSetBaserateInflation = 0,
20		OcSLength = 10000,
21		replicates = 800,
22		ScSKappaThreshold = 0.9,
23		ScSKappaMin = .40,
24		ScSPrecisionMin = 0.6,
25		ScSPrecisionMax = 1,
26		method = "standard"
27		){
28	6x	if (any(x < 0)) {
29	1x	stop("Values in Contingency Table must be positive")
30		}
31	5x	kappa = kappa_ct(x);
32
33	5x	if (is.null(OcSBaserate)) {
34	5x	OcSBaserate = (sum(x[1,])) / (sum(x))
35		}
36
37	5x	rho_result <- rhoK(
38	5x	kappa,
39	5x	OcSBaserate,
40	5x	sum(x),
41	5x	testSetBaserateInflation,
42	5x	OcSLength,
43	5x	replicates,
44	5x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax
45		)
46
47	5x	return(list(
48	5x	rho = rho_result,
49	5x	kappa = kappa,
50	5x	recall = rating_set_recall(x),
51	5x	precision = rating_set_precision(x)
52		))
53		}

1		###
2		#' @title Rho using a file
3		#'
4		#' @description
5		#' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}} as defined by the file
6		#' and columns (col1, col2), and returns a list containing both values. Called by \code{\link{rho}}.
7		#'
8		#' @export
9		#'
10		#' @param col1 The first column from file
11		#' @param col2 The second column from file
12		#' @template rho-params
13		#'
14		#' @return A list of the format:\describe{
15		#' \item{rho}{The rho of the \code{\link{codeSet}}}
16		#' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
17		#' }
18		###
19		rho.file <- function(
20		x, col1, col2,
21		OcSBaserate = NULL,
22		testSetBaserateInflation = 0,
23		OcSLength = 10000,
24		replicates = 800,
25		ScSKappaThreshold = 0.9,
26		ScSKappaMin = .40,
27		ScSPrecisionMin = 0.6,
28		ScSPrecisionMax = 1,
29		method = "standard"
30		) {
31	3x	setFile = utils::read.csv(x)
32	3x	set = as.matrix(setFile[,c(col1, col2)])
33
34	3x	if (!any(is.na(set))) {
35	2x	rhoSet(
36	2x	set,
37	2x	OcSBaserate,
38	2x	testSetBaserateInflation,
39	2x	OcSLength, replicates,
40	2x	ScSKappaThreshold, ScSKappaMin,
41	2x	ScSPrecisionMin, ScSPrecisionMax
42		)
43		}
44		else {
45	1x	stop("The columns provided for col1 and col2 must not contain NA")
46		}
47		}

1		###
2		#' @title Rho
3		#'
4		#' @description
5		#' This function calculates rho for a \code{\link[=getTestSet]{testSet}}, \code{\link{contingencyTable}}, or an observed kappa value with associated set parameters (testSetLength and OcSBaserate).
6		#'
7		#' @details
8		#' Rho is a Monte Carlo rejective method of interrater reliability statistics, implemented here for Cohen's Kappa. Rho constructs a collection of data sets in which kappa is below a specified threshold, and computes the empirical distribution on kappa based on the specified sampling procedure. Rho returns the percent of the empirical distribution greater than or equal to an observed kappa. As a result, Rho quantifies the type 1 error in generalizing from an observed test set to a true value of agreement between two raters.
9		#'
10		#' Rho starts with an observed kappa value, calculated on a subset of a \code{\link{codeSet}}, known as an observed \code{\link[=getTestSet]{testSet}}, and a \emph{kappa threshold} which indicates what is considered significant agreement between raters.
11		#'
12		#' It then generates a collection of fully-coded, simulated
13		#' \code{\link[=getTestSet]{codeSets}} (ScS), further described in
14		#' \code{\link{createSimulatedCodeSet}}, all of which have a kappa value below the kappa
15		#' threshold and similar properties as the original \code{\link{codeSet}}.
16		#'
17		#' Then, kappa is calculated on a \code{\link[=getTestSet]{testSet}} sampled from each of the ScSs in the
18		#' collection to create a null hypothesis distribution. These \code{\link[=getTestSet]{testSets}} mirror the observed \code{\link[=getTestSet]{testSets}} in their size and sampling method. How these \code{\link[=getTestSet]{testSets}} are sampled is futher described in \code{\link{getTestSet}}.
19		#'
20		#' The null hypothesis is that the observed \code{\link[=getTestSet]{testSet}}, was sampled from a data set, which, if both raters were to code in its entirety, would result in a level of agreement below the kappa threshold.
21		#'
22		#' For example, using an alpha level of 0.05, if the observed kappa is greater than 95 percent of the kappas in the null hypothesis distribution, the null hypothesis is rejected. Then one can conclude that the two raters would have acceptable agreement had they coded the entire data set.
23		#'
24		#' @template rho-params
25		#' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
26		#'
27		#' @examples
28		#' # Given an observed kappa value
29		#' rho(x = 0.88, OcSBaserate = 0.2, testSetLength = 80)
30		#'
31		#' # Given a test Set
32		#' rho(x = codeSet)
33		#'
34		#' # Given a contingency Table
35		#' rho(x = contingencyTable)
36		#'
37		#' @export
38		#' @return rho and kappa for the given data and parameters (unless kappa is given)
39		###
40		rho = function(
41		x,
42		OcSBaserate = NULL,
43		testSetLength = NULL,
44		testSetBaserateInflation = 0,
45		OcSLength = 10000,
46		replicates = 800,
47		ScSKappaThreshold = 0.9,
48		ScSKappaMin = 0.40,
49		ScSPrecisionMin = 0.6,
50		ScSPrecisionMax = 1,
51		method = "standard"
52		) {
53
54		### Calculate rho given a kappa value
55	22x	if (is(x, "numeric")) {
56	1x	if (is.null(OcSBaserate)) { stop("Must give a baserate when a kappa value is given") }
57	1x	if (is.null(testSetLength)) { stop("Must give a testSetLength when a kappa value is given") }
58
59	13x	result <- rhoK(x, OcSBaserate,
60	13x	testSetLength,
61	13x	testSetBaserateInflation,
62	13x	OcSLength,
63	13x	replicates,
64	13x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
65		}
66
67		### Calculate rho given a code set
68	7x	else if (nrow(x) != 2) {
69	3x	result <- rhoSet(x, OcSBaserate,
70	3x	testSetBaserateInflation,
71	3x	OcSLength,
72	3x	replicates,
73	3x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
74		}
75
76		### Calculate rho given a contingency table
77		else {
78	4x	result <- rhoCT(x, OcSBaserate,
79	4x	testSetBaserateInflation,
80	4x	OcSLength,
81	4x	replicates,
82	4x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
83		}
84
85	20x	return(result)
86		}

1		genPKcombo = function(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate){
2		# currKappa = sample(kappaDistribution, 1, replace = FALSE, prob = kappaProbability)
3		# currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
4	12132x	currKappa = kappaDistribution[sample.int(length(kappaDistribution), 1, replace = FALSE, prob = kappaProbability)]
5	12131x	currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]
6
7	12131x	if(!checkBRPKcombo(baserate, currPrecision, currKappa)){
8	2057x	precisionMin = (2baseratecurrKappa - 2*baserate - currKappa)/(currKappa - 2)
9	2057x	indicies = which(precisionDistribution > precisionMin)
10
11	2057x	if (length(indicies) == 0) {
12	2048x	return(genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate))
13		}
14
15	9x	precisionDistribution = precisionDistribution[indicies]
16	9x	precisionProbability = precisionProbability[indicies]
17		# currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
18	9x	currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]
19	9x	return(list(precision = currPrecision, kappa = currKappa))
20		} else {
21	10074x	return(list(precision = currPrecision, kappa = currKappa))
22		}
23		}

1		###
2		# Generate Precision Combo
3		#
4		# This function generates a precision that fits within the range and will create a valid set
5		# @param kappa This is the kappa of the desired pair
6		# @param Pmin This is the minimum precision desired
7		# @param Pmax This is the maximum precision desired
8		# @param BR This is the baserate of the desired pair
9		# @return This returns a precision value that will generate a valid set with the other parameters
10		####
11
12		genPcombo = function(kappa, Pmin, Pmax, BR){
13		#gets a random precision in the precision range
14	!	currPrecision = stats::runif(1, Pmin, Pmax);
15		#checks is the combination of precision, kappa, and baserate is valid. If it is vaid, return the precision, if not, recall the function until a valid precision is found
16	!	if(!checkBRPKcombo(BR, currPrecision, kappa)){
17	!	genPcombo(kappa, Pmin, Pmax, BR);
18		}else{
19	!	return(currPrecision);
20		}
21		}

1		###
2		#' @title Calculate Baserate (CT)
3		#'
4		#' @description
5		#' This function calculates the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
6		#'
7		#' @export
8		#'
9		#' @param CT The \code{\link{contingencyTable}} for which the baserate is calculated
10		#'
11		#' @seealso \code{\link{baserate}} and \code{\link{baserateSet}}
12		#'
13		#' @return A list of the format:\describe{
14		#' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
15		#' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
16		#' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
17		#' }
18		#'
19		###
20
21		baserateCT = function(CT){
22	1x	if(any(CT < 0)){stop("Values in Contingency Table must be positive")}
23	35x	firstBaserate = sum(CT[1,])/sum(CT)
24	35x	secondBaserate = sum(CT[,1])/sum(CT)
25	35x	averageBaserate = mean(c(firstBaserate,secondBaserate))
26	35x	return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
27		}

1		###
2		# Get Boot P Value
3		#
4		# This function gets the p value of the result in the distribution
5		# @param distribution The distribution that the result falls into
6		# @param result The value that you are evaluating where in the distribution it falls
7		# @return This returns the p value from the result into the distribution
8		###
9		getBootPvalue = function(distribution, result) {
10		#returns the percentage of the time that the distribution was greater or equal to the observed kappa
11	6x	N = length(distribution);
12	6x	if(result < mean(distribution, na.rm = T)) {
13		#if the result is less than the mean of the distribution, than the p value is 1
14	2x	return(1);
15		} else {
16		#return the number of times that the distribution is greater than the result as a percentage of the total number of items in the distribution
17	4x	return(length(which(distribution >= result)) / N);
18		}
19		}

1		#' Generate a Handset
2		#'
3		#' @description Generate a vector representing indices of set, using the handSetBaserate
4		#' to determine the minimum number of indices that are positive
5		#'
6		#' @param set matrix of two columns
7		#' @param handSetLength number of indices to find
8		#' @param handSetBaserate number between 0 and 1 to use as a minimum number of positive indices
9		#'
10		#' @return vector of indices from set
11		#' @export
12		getHandSetIndices = function(set, handSetLength = 20, handSetBaserate = 0.2) {
13	5x	positives = ceiling(handSetLength * handSetBaserate);
14	5x	posInd = which(set[, 1] == 1);
15
16	5x	if (positives > length(posInd)) {
17	1x	stop("Not enough positives in first rater to inflate to this level")
18		}
19
20	4x	this.set = NULL
21	4x	if (positives > 0) {
22	3x	positiveIndices = posInd[sample.int(length(posInd), size=positives, replace=FALSE)]
23	3x	others = which(!(1:nrow(set) %in% positiveIndices))
24	3x	otherIndices = sample.int(length(others), size = (handSetLength - positives), replace = FALSE)
25	3x	this.set = c(positiveIndices, others[otherIndices])
26		}
27	1x	else if (positives == 0) {
28	1x	theseIndices = sample.int(length(set),size=handSetLength,replace=FALSE)
29	1x	this.set = theseIndices
30		}
31
32	4x	return(sample(this.set))
33		}

1		###
2		#' @title Calculate Baserate (Set)
3		#'
4		#' @description
5		#' This function will calculate the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
6		#'
7		#' @param set The \code{\link{codeSet}} for which the baserate is calculated
8		#'
9		#' @seealso \code{\link{baserate}} and \code{\link{baserateCT}}
10		#'
11		#' @export
12		#' @return A list of the format:\describe{
13		#' \item{firstBaserate}{The percentage that a positive code, or a 1, appears in the first rater}
14		#' \item{secondBaserate}{The percentage that a positive code, or a 1, appears in the second rater}
15		#' \item{averageBaserate}{The average percentage that a positive code, or a 1, appears in either of the two raters}
16		#' }
17		#'
18		###
19		baserateSet = function(set){
20	12x	firstBaserate = length(which(set[,1] == 1))/nrow(set)
21	12x	secondBaserate = length(which(set[,2] == 1))/nrow(set)
22	12x	averageBaserate = mean(c(firstBaserate,secondBaserate))
23	12x	return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
24		}

1		###
2		# Calculate Kappa
3		#
4		# This function calculates kappa given a set
5		# @param set This is the set to calculate kappa from
6		# @param isSet TRUE if set FALSE if contigency table
7		# @param kappaThreshold if null then return kappa otherwise return list of kappa and above or not
8		# @return This returns the kappa of the set given
9		###
10		calcKappa <- function(set, isSet = TRUE, kappaThreshold = NULL) {
11	126x	if (!isSet) {
12	47x	set <- contingencyToSet(set[1, 1], set[2, 1], set[1, 2], set[2, 2])
13		}
14
15		if (
16	126x	length(which(set[, 1] == set[, 2] & set[, 1] == 1)) == nrow(set) \|
17	126x	length(which(set[, 1] == set[, 2] & set[, 1] == 0)) == nrow(set)
18		) {
19	4x	return(1)
20		}
21
22		#calculates kappa by calculating the agreement and the baserates and then creating the adjacenty matrix
23	122x	agreement <- length(which(set[, 1] == set[, 2])) / nrow(set);
24	122x	baserate2 <- length(which(1 == set[, 2])) / nrow(set);
25	122x	baserate1 <- length(which(1 == set[, 1])) / nrow(set);
26	122x	randomAgreement <- baserate1 * baserate2 + (1 - baserate1) * (1 - baserate2);
27	122x	kappa <- (agreement - randomAgreement) / (1 - randomAgreement);
28
29	122x	if (is.null(kappaThreshold)) {
30	120x	return(kappa);
31		} else {
32	2x	return(list(kappa = kappa, above = kappa > kappaThreshold))
33		}
34		}

1		###
2		# Check Baserate, Precision, Kappa Combo
3		#
4		# This function checks to make sure the combination of BR, P, and K will create a valid set
5		# @param BR This is the supplied baserate
6		# @param P This is the supplied precision
7		# @param K This is the supplied kappa
8		# @return returns TRUE if the combination is valid, FALSE otherwise
9		###
10		checkBRPKcombo = function(BR, P, K) {
11		#if right is less than P, then it is a valid combination. If not, then the pair of kappa, precision, and baserate does not correspond to a valid set
12	12135x	right = (2BRK - 2*BR - K)/(K - 2);
13	12135x	if(P > right){
14	10076x	return(TRUE);
15		}else{
16	2059x	return(FALSE);
17		}
18		}

1		###
2		# This function calculates rho given the parameters of the initial set and the kappa of the observed handset and is
3		# called by rhoK()
4		#
5		# The computation is libqe-stats' qestats::rho_standard (src/rhor.cpp): the
6		# same draws and arithmetic as this function's former R implementation
7		# (generateKPs, getR, prset, getHandSet, getBootPvalue), so seeded results
8		# are unchanged. KPs, when given, is a list of genPKcombo() results.
9		###
10		calcRho = function(
11		x,
12		OcSBaserate,
13		testSetLength,
14		testSetBaserateInflation = 0,
15		OcSLength = 10000,
16		replicates = 800,
17		ScSKappaThreshold = 0.9,
18		ScSKappaMin = 0.40,
19		ScSPrecisionMin = 0.6,
20		ScSPrecisionMax = 1,
21		KPs = NULL
22		) {
23	159x	if (!is.null(KPs)) {
24	100x	KPs = cbind(
25	100x	vapply(KPs, function(KP) KP$precision, numeric(1)),
26	100x	vapply(KPs, function(KP) KP$kappa, numeric(1))
27		)
28		}
29
30	159x	return(calcRho_standard_c(
31	159x	x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates,
32	159x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax, KPs
33		));
34		}

1		rating_set_recall <- function (x) {
2	14x	if (!is(x, "contingency.table")) {
3	14x	x <- as.contingency.table(x)
4		}
5
6	14x	x$agreement$true_positive / (x$agreement$true_positive + x$agreement$false_negative)
7		}
8
9		rating_set_precision <- function (x) {
10	14x	if (!is(x, "contingency.table")) {
11	14x	x <- as.contingency.table(x)
12		}
13
14	14x	x$agreement$true_positive / (x$agreement$true_positive + x$agreement$false_positive)
15		}

1		###
2		# Create A Random Set
3		#
4		# This creates a set with given parameters
5		# @param setLength The length of the set to be created
6		# @param baserate The baserate of the set to be created
7		# @param kappa The kappa of the set to be created
8		# @param minPrecision The minimum precision accepted in this set. Defaults to 0
9		# @param maxPrecision The maximum precision accepted in this set. Defaults to 1
10		# @return This returns a set with the given parameters
11		###
12		createRandomSet = function(setLength, baserate, kappaMin, kappaMax, minPrecision = 0, maxPrecision = 1, type = "set"){
13		#generates a valid precision with the given kappa and baserate
14	56x	KP = generateKPs(1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision)[[1]]
15	56x	kappa = KP$kappa
16	56x	precision = KP$precision
17	56x	recall = getR(kappa, baserate, precision)
18
19	56x	if(type == "ct") {
20	2x	fullSet = pr.ct(precision = precision, recall = recall, length = setLength, baserate = baserate);
21		}
22		else {
23	54x	fullSet = prset(precision = precision, recall = recall, length = setLength, baserate = baserate);
24		}
25
26	56x	return(fullSet);
27		}

1		###
2		#' @title Rho Min
3		#'
4		#' @description
5		#' This function calculates the minimum testSetLength where it is possible to get a rho less than alpha for the given parameters of rho.
6		#'
7		#' @param baserate A \code{\link{baserate}}
8		#' @param alpha The threshold of significance for rho (similar to an alpha level for a p value), defaulted to 0.05
9		#' @param inc An integer indicating by how much the testSetLength should increase each iteration
10		#' @param printInc A boolean indicating whether to print out each increment value with it's corresponding significance for rho
11		#' @param ... Any additional parameters passed into \code{\link{rho}}
12		#'
13		#' @examples
14		#' #Add testSetBaserateInflation as an additional parameter
15		#' rhoMin(0.2, testSetBaserateInflation = 0.33)
16		#'
17		#' #Add testSetBaserateInflation as well as changing inc and selecting printInc
18		#' rhoMin(0.2, inc = 5, printInc = TRUE, testSetBaserateInflation = 0.33)
19		#'
20		#' @export
21		#' @return The minimum length of testSet, to the nearest multiple of inc, greater than the minimum length, that would give a value where rho less than alpha becomes mathematically possible.
22		###
23		rhoMin = function(baserate, alpha = 0.05, inc = 10, printInc = FALSE, ...){
24	5x	rhoVal = 1
25	5x	testSetLength = 0
26	1x	if(inc%%1 != 0 \| inc <= 0) stop("Inc value must be a positive integer.")
27	1x	if(alpha > 1 \| alpha < 0) stop("Alpha must be between 0 and 1.")
28
29
30	3x	while(rhoVal > alpha){
31	12x	testSetLength = testSetLength + inc
32	12x	rhoVal = rho(1, baserate, testSetLength, ...)
33	12x	if(printInc) {
34	4x	cat(testSetLength, rhoVal,"\n")
35		}
36		}
37	3x	return(testSetLength)
38
39		}

1		###
2		#' @title Calculate Baserate
3		#'
4		#' @description
5		#' This function calculates the baserate of the first rater, second rater, and the average of both the raters.
6		#'
7		#' @details
8		#' A baserate is the percentage, as a decimal, that a positive code appears in data (either a \code{\link{codeSet}} or \code{\link{contingencyTable}}) for a given rater. It is assumed that the first rater is more experienced and thus provides a better estimation of the actual baserate for a given code, so the first rater's baserate is often used as if it is the actual baserate. If the raters are assumed to have the same experience level, the average baserate may give a better estimation. If the second rater is more experienced, the second rater's baserate may give a better estimation. Functions assume that the first rater is the more experienced rater and thus uses the first rater's baserate as the overall baserate estimation.
9		#'
10		#'
11		#' @param data The \code{\link[=getTestSet]{testSet}} or \code{\link{contingencyTable}} for which the baserate is calculatede
12		#'
13		#' @seealso \code{\link{baserateSet}} and \code{\link{baserateCT}}
14		#'
15		#' @examples
16		#' #Given a code set
17		#' baserate(data = codeSet)
18		#'
19		#' #Given a contingency Table
20		#' baserate(data = contingencyTable)
21		#'
22		#' @export
23		#' @return A list of the format:\describe{
24		#' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
25		#' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
26		#' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
27		#' }
28		#'
29		###
30		baserate = function(data){
31	7x	if (nrow(data) != 2) {
32	4x	return(baserateSet(data))
33		} else {
34	3x	return(baserateCT(data))
35		}
36		}

1		###
2		#' @title Calculate kappa (contingency Table)
3		#'
4		#' @description
5		#' This function calculates Cohen's kappa on a \code{\link{contingencyTable}}. Called by \code{\link{kappa}}.
6		#'
7		#' @param ct A \code{\link{contingencyTable}}
8		#'
9		#' @seealso \code{\link{kappa}} and \code{\link{kappaSet}}
10		#'
11		#' @export
12		#' @return The kappa of the \code{\link{contingencyTable}}
13		###
14		kappaCT = function(ct){
15	1x	if(any(ct < 0)){stop("Values in Contingency Table must be positive")}
16	1x	if(ncol(ct) != 2 \|\| nrow(ct) != 2) stop("Incorrect number of dimensions: Contingency table must be 2x2.")
17	42x	for(i in 1:nrow(ct)) {
18	83x	for(j in 1:ncol(ct)) {
19	1x	if (ct[i,j]%%1 != 0) stop("Contingency table values must be positive integers.")
20		}
21		}
22
23	41x	return(calcKappa(ct,isSet = FALSE,kappaThreshold = NULL));
24		}

1		###
2		# Get Recall
3		#
4		# This gets the recall from a pair
5		# @param kappa This is the kappa of the pair
6		# @param BR This is the baserate of the pair
7		# @param P this is the precision of the pair
8		# @return returns the recall given the parameters
9		###
10		getR = function(kappa, BR, P){
11		#gets recall that corresponds with the kappa, baserate, and precision
12	63x	top = kappa*P;
13	63x	R = top / (2P-2BR-kappa+2BRkappa);
14	63x	return (R);
15		}

1		###
2		#' @title Calculate kappa (Set)
3		#'
4		#' @description
5		#' This function calculates Cohen's kappa for a given \code{\link{codeSet}}. Called by \code{\link{kappa}}.
6		#'
7		#' @export
8		#' @param set A \code{\link{codeSet}}
9		#'
10		#' @seealso \code{\link{kappa}} and \code{\link{kappaCT}}
11		#'
12		#' @return
13		#' The kappa of the \code{\link{codeSet}}
14		#'
15		###
16		kappaSet = function(set){
17	76x	if(any(set != 1 & set != 0)) {
18	1x	stop("Each set must only consist of zeros and ones.")
19		}
20
21	75x	return(calcKappa(set,isSet = TRUE,kappaThreshold = NULL));
22		}
23

1		###
2		#' @include rhoK.R
3		#' @title Rho (set)
4		#'
5		#' @description
6		#' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}}, and returns a list containing both values. Called by \code{\link{rho}}.
7		#'
8		#' @export
9		#'
10		#' @template rho-params
11		#'
12		#' @return A list of the format:\describe{
13		#' \item{rho}{The rho of the \code{\link{codeSet}}}
14		#' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
15		#' }
16		###
17		rhoSet = function(
18		x,
19		OcSBaserate = NULL,
20		testSetBaserateInflation = 0,
21		OcSLength = 10000,
22		replicates = 800,
23		ScSKappaThreshold = 0.9,
24		ScSKappaMin = 0.40,
25		ScSPrecisionMin = 0.6,
26		ScSPrecisionMax = 1,
27		method = "standard"
28		) {
29	8x	kappa = kappaSet(x);
30
31	8x	row.count = nrow(x);
32	8x	if(is.null(OcSBaserate)) {
33	7x	OcSBaserate = length(which(x[,1] == 1)) / row.count;
34		}
35
36	8x	return(list(
37	8x	rho = rhoK(
38	8x	kappa,
39	8x	OcSBaserate,
40	8x	row.count,
41	8x	testSetBaserateInflation,
42	8x	OcSLength,
43	8x	replicates,
44	8x	ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax,
45	8x	method = method
46		),
47	8x	kappa = kappa,
48	8x	recall = rating_set_recall(x),
49	8x	precision = rating_set_precision(x)
50		))
51		}

1		###
2		#' @title Calculate kappa
3		#'
4		#' @description
5		#' This function calculates Cohen's kappa on a \code{\link{contingencyTable}} or a \code{\link{codeSet}}
6		#'
7		#' @param data A \code{\link{contingencyTable}} or a \code{\link{codeSet}}
8		#'
9		#'
10		#' @seealso \code{\link{kappaSet}} and \code{\link{kappaCT}}
11		#'
12		#' @examples
13		#' #Given a code set
14		#' kappa(data = codeSet)
15		#'
16		#' #Given a contingency Table
17		#' kappa(data = contingencyTable)
18		#'
19		#' @export
20		#' @return The kappa of the \code{\link{contingencyTable}} or \code{\link{codeSet}}
21		###
22		kappa <- function(data) {
23	62x	if (all(dim(data) == 2)) {
24	5x	return(kappaCT(data))
25		} else {
26	57x	return(kappaSet(data))
27		}
28		}

1		#pragma once
2		#include <algorithm>
3		#include <cmath>
4		#include <stdexcept>
5		#include <string>
6
7		// Functions that must reproduce R's arithmetic (rhoR parity) start their body
8		// with QESTATS_R_ARITHMETIC and carry QESTATS_R_ARITHMETIC_FN: R never fuses
9		// a * b + c, so these keep the compiler from contracting it into an FMA.
10		#if defined(__clang__)
11		# define QESTATS_R_ARITHMETIC _Pragma("STDC FP_CONTRACT OFF")
12		# define QESTATS_R_ARITHMETIC_FN
13		#elif defined(__GNUC__)
14		# define QESTATS_R_ARITHMETIC
15		# define QESTATS_R_ARITHMETIC_FN __attribute__((optimize("fp-contract=off")))
16		#else
17		# define QESTATS_R_ARITHMETIC
18		# define QESTATS_R_ARITHMETIC_FN
19		#endif
20
21		namespace qestats {
22
23		// Confusion matrix proportions (must sum to 1 after normalization).
24		struct ConfusionMatrix {
25		double tp = 0, fp = 0, fn = 0, tn = 0;
26		};
27
28		// All derived classification metrics.
29		struct Metrics {
30		// confusion matrix proportions
31		double tp = 0, fp = 0, fn = 0, tn = 0;
32		// derived
33		double b1 = 0; // classifier positive rate (tp + fp)
34		double b2 = 0; // true positive rate (tp + fn)
35		double precision = 0; // tp / (tp + fp)
36		double recall = 0; // tp / (tp + fn)
37		double accuracy = 0; // tp + tn
38		double fdr = 0; // fp / (tp + fp) — False Discovery Rate
39		double for_ = 0; // fn / (fn + tn) — False Omission Rate
40		double kappa = 0;
41		};
42
43		// Range result — min/max of a metric over the key variable interval.
44		struct MetricRange {
45		double min = 0;
46		double max = 0;
47		};
48
49		// Confidence interval from a hypergeometric estimate.
50		struct Interval {
51		double lower = 0;
52		double mean = 0;
53		double upper = 0;
54		};
55
56		// One row in the metrics_range output table.
57		struct MetricRangeRow {
58		std::string metric;
59		double lower = 0, mean = 0, upper = 0;
60		};
61
62		// ── Core arithmetic ───────────────────────────────────────────────────────────
63
64		// Cohen's kappa from confusion-matrix counts (or proportions), computed in
65		// the same order as rhoR's calcKappa() so the two agree to the bit: perfect
66		// agreement on one class (all TP or all TN) is 1.
67		QESTATS_R_ARITHMETIC_FN
68	51800x	inline double kappa_from_counts(double tp, double fp, double fn, double tn) {
69		QESTATS_R_ARITHMETIC
70	51800x	const double s = tp + fp + fn + tn;
71	51800x	if (tp == s \|\| tn == s) return 1.0;
72	50850x	const double agreement = (tp + tn) / s;
73	50850x	const double rate1 = (tp + fn) / s; // first rater's positive rate
74	50850x	const double rate2 = (tp + fp) / s; // second rater's positive rate
75	50850x	const double pe = rate1 * rate2 + (1 - rate1) * (1 - rate2);
76		return (pe < 1.0) ? (agreement - pe) / (1 - pe) : 0.0;
77		}
78
79		// Rates are computed from the counts, not the normalized cells, so precision
80		// and recall match rhoR's TP / (TP + FP) and TP / (TP + FN) exactly.
81		inline Metrics metrics_from_cm(const ConfusionMatrix& cm) {
82		const double s = cm.tp + cm.fp + cm.fn + cm.tn;
83		if (s <= 0) throw std::domain_error("metrics_from_cm: confusion matrix is empty");
84		Metrics m;
85		m.tp = cm.tp / s; m.fp = cm.fp / s;
86		m.fn = cm.fn / s; m.tn = cm.tn / s;
87		const double pos1 = cm.tp + cm.fp, pos2 = cm.tp + cm.fn, neg1 = cm.fn + cm.tn;
88		m.b1 = pos1 / s;
89		m.b2 = pos2 / s;
90		m.precision = (pos1 > 0) ? cm.tp / pos1 : 0.0;
91		m.recall = (pos2 > 0) ? cm.tp / pos2 : 0.0;
92		m.accuracy = (cm.tp + cm.tn) / s;
93		m.fdr = (pos1 > 0) ? cm.fp / pos1 : 0.0;
94		m.for_ = (neg1 > 0) ? cm.fn / neg1 : 0.0;
95		m.kappa = kappa_from_counts(cm.tp, cm.fp, cm.fn, cm.tn);
96		return m;
97		}
98
99		// Convenience: from raw counts (negative values clamped to 0).
100		inline Metrics calculate_metrics(double tp, double fp, double fn, double tn) {
101		if (tp < 0) tp = 0; if (fp < 0) fp = 0;
102		if (fn < 0) fn = 0; if (tn < 0) tn = 0;
103		return metrics_from_cm({tp, fp, fn, tn});
104		}
105
106		// Retrieve a single metric value by name (supports alias names).
107		inline double get_metric(const Metrics& m, const std::string& name) {
108		if (name == "b1" \|\| name == "baserate") return m.b1;
109		if (name == "b2" \|\| name == "baserate2") return m.b2;
110		if (name == "precision") return m.precision;
111		if (name == "recall") return m.recall;
112		if (name == "accuracy") return m.accuracy;
113		if (name == "FDR") return m.fdr;
114		if (name == "FOR") return m.for_;
115		if (name == "kappa") return m.kappa;
116		if (name == "tp") return m.tp;
117		if (name == "fp") return m.fp;
118		if (name == "fn") return m.fn;
119		if (name == "tn") return m.tn;
120		throw std::invalid_argument("get_metric: unknown metric '" + name + "'");
121		}
122
123		} // namespace qestats

1		// rhoR's C++ functions, on libqe-stats: rho and its building blocks are
2		// qestats::rho_fast (this file's functions, as they were) and
3		// qestats::rho_standard (rhoK's default method), vendored in
4		// inst/include/qestats/ by tools/sync-qestats.sh. Given R's RNG they return
5		// exactly what rhoR returned before the move (tests/testthat/test-regression-pins.R).
6		//
7		// Only the bell-shaped kappa distribution (generate_kp_list(distributionType = 1))
8		// is still computed here, with RcppArmadillo's weighted sample().
9
10		#include <RcppArmadilloExtensions/sample.h>
11		// [[Rcpp::depends(RcppArmadillo)]]
12
13		#include <RcppArmadillo.h>
14		#include <Rcpp.h>
15		#include "qestats/rho.hpp"
16		using namespace Rcpp;
17
18		namespace {
19
20		// R's RNG for qestats: unif_rand() for the fast method, R_unif_index() (what
21		// sample.int() draws with) for the standard one.
22		struct RRng {
23	3180904x	double unif_rand() { return ::unif_rand(); }
24	4194943x	double unif_index(double n) { return R_unif_index(n); }
25		};
26
27		// rhoR's tables are gold (rows) by silver (columns): [1,1] TP, [1,2] FN,
28		// [2,1] FP, [2,2] TN.
29	18x	qestats::ConfusionMatrix to_cm(const IntegerMatrix& ct) {
30	18x	return {double(ct(0, 0)), double(ct(1, 0)), double(ct(0, 1)), double(ct(1, 1))};
31		}
32
33	11x	IntegerMatrix to_ct(const qestats::ConfusionMatrix& cm) {
34	11x	IntegerMatrix ct(2, 2);
35	11x	ct(0, 0) = cm.tp;
36	11x	ct(0, 1) = cm.fn;
37	11x	ct(1, 0) = cm.fp;
38	11x	ct(1, 1) = cm.tn;
39	11x	return ct;
40		}
41
42	200x	std::vector<qestats::KappaPrecision> to_kps(const NumericMatrix& KPs) {
43	200x	std::vector<qestats::KappaPrecision> kps(KPs.nrow());
44	20200x	for (int i = 0; i < KPs.nrow(); i++) kps[i] = {KPs(i, 1), KPs(i, 0)};
45	200x	return kps;
46		}
47
48	2x	NumericMatrix from_kps(const std::vector<qestats::KappaPrecision>& kps) {
49	2x	NumericMatrix KPs(kps.size(), 2);
50	10007x	for (std::size_t i = 0; i < kps.size(); i++) {
51	10005x	KPs(i, 0) = kps[i].precision;
52	10005x	KPs(i, 1) = kps[i].kappa;
53		}
54	2x	return KPs;
55		}
56
57		// find_valid_pk() with kappa probabilities (the bell-shaped distribution).
58		// Precision is drawn unweighted, as before.
59	50007x	NumericVector find_valid_pk_weighted(
60		arma::colvec kappaDistribution, arma::vec kappaProbability,
61		arma::colvec precisionDistribution, double baserate
62		) {
63	50065x	for (int attempt = 0; attempt < qestats::RHO_MAX_ATTEMPTS; attempt++) {
64	50065x	arma::ivec kappaRng = arma::regspace<arma::ivec>(0, kappaDistribution.size() - 1);
65	50065x	arma::ivec whKappa = Rcpp::RcppArmadillo::sample(kappaRng, 1, false, kappaProbability);
66	100130x	double currKappa = kappaDistribution[whKappa.at(0)];
67
68	50065x	arma::ivec precRng = arma::regspace<arma::ivec>(0, precisionDistribution.size() - 1);
69	50065x	arma::ivec whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
70	100130x	double currPrec = precisionDistribution[whPrec.at(0)];
71
72	50065x	if (!qestats::rho_fast::valid_combo(baserate, currPrec, currKappa)) {
73	41218x	double precisionMin = (2 * baserate * currKappa - 2 * baserate - currKappa) / (currKappa - 2);
74	41218x	arma::uvec indicies = arma::find(precisionDistribution > precisionMin);
75	41218x	if (indicies.size() == 0) continue;
76	82320x	arma::colvec above = precisionDistribution.elem(indicies);
77	41160x	precRng = arma::regspace<arma::ivec>(0, above.size() - 1);
78	41160x	whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
79	41160x	currPrec = above[whPrec.at(0)];
80		}
81	100014x	return NumericVector::create(currPrec, currKappa);
82		}
83	!	stop(qestats::RHO_NO_VALID_SET);
84		}
85
86		} // namespace
87
88		//' @title sample_contingency_table
89		//' @name sample_contingency_table
90		//'
91		//' @param xx contingency table matrix
92		//' @param n int size of the contingency table
93		//' @param forR bool if true, add 1 to the results accounting for R indices starting at 1
94		//'
95		//' @export
96		// [[Rcpp::export]]
97	1x	IntegerMatrix sample_contingency_table(IntegerMatrix xx, int n, bool forR = true) {
98		RRng rng;
99		std::vector<int> cells = qestats::rho_fast::sample_contingency_table(
100	1x	std::vector<int>(xx.begin(), xx.end()), n, rng);
101	1x	IntegerMatrix ret(cells.size(), 1);
102		for (std::size_t w = 0; w < cells.size(); w++) ret[w] = cells[w] + (forR ? 1 : 0);
103	1x	return ret;
104		}
105
106		//' @name getBootPvalue_c
107		//' @title getBootPvalue_c
108		//'
109		//' @param distribution vector of calculated kappas
110		//' @param result double calculated kappa to compare against
111		//'
112		//' @description
113		//' returns the percentage of the time that the distribution was greater or equal to the observed kappa
114		//' if the result is less than the mean of the distribution, than the p value is 1
115		//' else return the number of times that the distribution is greater than the result as a percentage of the total number
116		//' of items in the distribution
117		//'
118		//' @return double calculated p-value
119		//'
120		//' @export
121		// [[Rcpp::export]]
122	5x	double getBootPvalue_c(NumericVector distribution, double result) {
123	5x	return qestats::rho_fast::boot_pvalue(
124	15x	std::vector<double>(distribution.begin(), distribution.end()), result);
125		}
126
127		// Validate a Baserate/Kappa combo against the supplied Precision.
128		// Return true if Precision is greater than the calculated right value
129		// [[Rcpp::export]]
130	1x	bool check_BRK_combo(double BR, double P, double K) {
131	1x	return qestats::rho_fast::valid_combo(BR, P, K);
132		}
133
134		//' @name recall
135		//' @title recall
136		//'
137		//' @param kappa double
138		//' @param BR double
139		//' @param P double
140		//'
141		//' @return Recall calculated from provided kappa, BR, and P
142		//'
143		//' @export
144		// [[Rcpp::export]]
145	3x	double recall(double kappa, double BR, double P) {
146	3x	return qestats::rho_fast::recall(kappa, BR, P);
147		}
148
149		// Find a valid precision and kappa combo from a distribution of Kappas and Precisions given a baserate
150		//
151		// [[Rcpp::export]]
152	2x	NumericVector find_valid_pk(
153		NumericVector kappaDistribution, NumericVector kappaProbability,
154		NumericVector precisionDistribution, NumericVector precisionProbability,
155		double baserate
156		) {
157		if (kappaDistribution.size() == 0 \|\| precisionDistribution.size() == 0) {
158	!	stop("find_valid_pk: empty kappa or precision distribution");
159		}
160	2x	if (kappaProbability.size() > 0) {
161	4x	return find_valid_pk_weighted(as<arma::colvec>(kappaDistribution), as<arma::vec>(kappaProbability),
162	6x	as<arma::colvec>(precisionDistribution), baserate);
163		}
164		RRng rng;
165	!	qestats::KappaPrecision kp = qestats::rho_fast::kappa_precision(
166	!	std::vector<double>(kappaDistribution.begin(), kappaDistribution.end()),
167	!	std::vector<double>(precisionDistribution.begin(), precisionDistribution.end()),
168		baserate, rng);
169	!	return NumericVector::create(kp.precision, kp.kappa);
170		}
171
172
173		//' @name generateKPs_c
174		//' @title generate_kp_list
175		//'
176		//' @param numNeeded int
177		//' @param baserate double
178		//' @param kappaMin double
179		//' @param kappaMax double
180		//' @param precisionMin double
181		//' @param precisionMax double
182		//' @param distributionType int 0 - normal (default), 1 - bell
183		//' @param distributionLength long
184		//'
185		//' @return matrix of kappa and precision values (column 1 as precision)
186		//'
187		//' @export
188		// [[Rcpp::export]]
189	4x	Rcpp::NumericMatrix generate_kp_list(
190		int numNeeded, double baserate,
191		double kappaMin, double kappaMax,
192		double precisionMin, double precisionMax,
193		int distributionType = 0, long distributionLength = 10000
194		) {
195	4x	if (distributionType != 1) {
196		RRng rng;
197	4x	return from_kps(qestats::rho_fast::generate_kps(
198	2x	numNeeded, baserate, kappaMin, kappaMax, precisionMin, precisionMax, rng, distributionLength));
199		}
200
201	2x	double kappaStep = ((kappaMax - kappaMin) / (distributionLength - 1));
202	2x	arma::colvec kappaDistribution = arma::regspace(kappaMin, kappaStep, kappaMax);
203	2x	arma::vec kappaProbability = arma::normpdf(kappaDistribution, 0.9, 0.1);
204	2x	double precStep = (precisionMax - precisionMin) / (10000 - 1);
205	2x	arma::colvec precisionDistribution = arma::regspace(precisionMin, precStep, precisionMax);
206
207	2x	Rcpp::NumericMatrix KPs(numNeeded, 2);
208	50007x	for (int i = 0; i < numNeeded; i++) {
209	50005x	KPs(i, _) = find_valid_pk_weighted(kappaDistribution, kappaProbability, precisionDistribution, baserate);
210		}
211	2x	return KPs;
212		}
213
214		//' @name contingency_table
215		//' @title contingency_table
216		//' @description Create a contingency table using the provied precision, recall, baserate, and length.
217		//'
218		//' @param precision double
219		//' @param rec double
220		//' @param length int
221		//' @param baserate double
222		//'
223		//' @export
224		// [[Rcpp::export]]
225	5x	IntegerMatrix contingency_table(double precision, double rec, int length, double baserate) {
226	5x	return to_ct(qestats::rho_fast::contingency_table(precision, rec, length, baserate));
227		}
228
229		//' @name random_contingency_table
230		//' @title random_contingency_table
231		//' @param setLength [TBD]
232		//' @param baserate [TBD]
233		//' @param kappaMin [TBD]
234		//' @param kappaMax [TBD]
235		//' @param minPrecision [TBD]
236		//' @param maxPrecision [TBD]
237		//'
238		//' @export
239		// [[Rcpp::export]]
240	2x	IntegerMatrix random_contingency_table(
241		int setLength,double baserate,
242		double kappaMin, double kappaMax,
243		double minPrecision = 0, double maxPrecision = 1
244		) {
245		RRng rng;
246	2x	qestats::KappaPrecision kp = qestats::rho_fast::generate_kps(
247	2x	1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision, rng)[0];
248	2x	double rec = qestats::rho_fast::recall(kp.kappa, baserate, kp.precision);
249	2x	return to_ct(qestats::rho_fast::contingency_table(kp.precision, rec, setLength, baserate));
250		}
251
252		//' @title kappa_ct
253		//' @description Calculate kappa from a contingency table
254		//' @param ct [TBD]
255		//'
256		//' @export
257		// [[Rcpp::export]]
258	13x	double kappa_ct(IntegerMatrix ct) {
259	13x	return qestats::rho_fast::kappa_ct(to_cm(ct));
260		}
261
262		// [[Rcpp::export]]
263	4x	IntegerMatrix getHand_ct(IntegerMatrix ct, int handSetLength, double handSetBaserate) {
264		RRng rng;
265	4x	return to_ct(qestats::rho_fast::hand_contingency_table(to_cm(ct), handSetLength, handSetBaserate, rng));
266		}
267
268		//' @name getHand_kappa
269		//' @title getHand_kappa
270		//'
271		//' @description This function returns kappa calculated from a Handset taken from a larger Contingency Table
272		//'
273		//' @param ct KPs matrix of kappa (column 1) and precision (column 2) values
274		//' @param handSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
275		//' @param handSetBaserate baserate to inflate the sampled contingency table to
276		//'
277		//' @return Kappa as double
278		//'
279		//' @export
280		// [[Rcpp::export]]
281	1x	double getHand_kappa(IntegerMatrix ct, int handSetLength, double handSetBaserate) {
282		RRng rng;
283	1x	return qestats::rho_fast::kappa_ct(
284	2x	qestats::rho_fast::hand_contingency_table(to_cm(ct), handSetLength, handSetBaserate, rng));
285		}
286
287		// This function calculates rho given the parameters of the initial set and
288		// the kappa of the observed handset and is called by rhoK()
289		//
290		// [[Rcpp::export]]
291	134x	double calcRho_c(
292		double x,
293		double OcSBaserate,
294		int testSetLength,
295		double testSetBaserateInflation = 0,
296		int OcSLength = 10000,
297		int replicates = 800,
298		double ScSKappaThreshold = 0.9,
299		double ScSKappaMin = 0.40,
300		double ScSPrecisionMin = 0.6,
301		double ScSPrecisionMax = 1.0,
302		NumericMatrix KPs = NumericMatrix(0)
303		) {
304		RRng rng;
305	134x	qestats::RhoOptions opts;
306	134x	opts.test_set_baserate_inflation = testSetBaserateInflation;
307	134x	opts.ocs_length = OcSLength;
308	134x	opts.replicates = replicates;
309	134x	std::vector<qestats::KappaPrecision> kps;
310	134x	if (KPs.size() == 1 && KPs[0] == 0) {
311	68x	kps = qestats::rho_fast::generate_kps(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold,
312	34x	ScSPrecisionMin, ScSPrecisionMax, rng);
313		} else {
314	100x	kps = to_kps(KPs);
315		}
316	268x	return qestats::rho_fast::rho(x, OcSBaserate, testSetLength, kps, opts, rng);
317		}
318
319		// calcRho()'s computation (rhoK's default method): qestats::rho_standard.
320		// KPs, when given, is a matrix of precision (column 1) and kappa (column 2).
321		// [[Rcpp::export]]
322	159x	double calcRho_standard_c(
323		double x,
324		double OcSBaserate,
325		int testSetLength,
326		double testSetBaserateInflation,
327		int OcSLength,
328		int replicates,
329		double ScSKappaThreshold,
330		double ScSKappaMin,
331		double ScSPrecisionMin,
332		double ScSPrecisionMax,
333		Nullable<NumericMatrix> KPs = R_NilValue
334		) {
335		RRng rng;
336	159x	qestats::RhoOptions opts;
337	159x	opts.test_set_baserate_inflation = testSetBaserateInflation;
338	159x	opts.ocs_length = OcSLength;
339	159x	opts.replicates = replicates;
340	159x	std::vector<qestats::KappaPrecision> kps;
341	159x	if (KPs.isNull()) {
342	118x	kps = qestats::rho_standard::generate_kps(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold,
343	59x	ScSPrecisionMin, ScSPrecisionMax, rng);
344		} else {
345	100x	kps = to_kps(NumericMatrix(KPs.get()));
346		}
347	318x	return qestats::rho_standard::rho(x, OcSBaserate, testSetLength, kps, opts, rng);
348		}