options(stringsAsFactors = FALSE)

# library("clusterSim")

list.of.packages <- c("PRROC", "e1071")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]

if(length(new.packages)) install.packages(new.packages)

library("e1071")

library("PRROC")

source("./utils.r")

# Confusion matrix rates

confusion_matrix_rates <- function (actual_labels, predicted_values, keyword)

{

    fg_test <- predicted_values[actual_labels==1]

    bg_test <- predicted_values[actual_labels==0]

    pr_curve_test <- pr.curve(scores.class0 = fg_test, scores.class1 = bg_test, curve = F)

    # plot(pr_curve_test)

    # print(pr_curve_test)

    prc_auc <- pr_curve_test$auc.integral

    cat("\nPR AUC (integral) \t", prc_auc, "\n", sep="")    

    # cat("PRC AUC (Davis & Goadrich) ", pr_curve_test$auc.davis.goadrichl, "\n", sep="")

    roc_curve_test <- roc.curve(scores.class0 = fg_test, scores.class1 = bg_test, curve = F)

    # plot(pr_curve_test)

    # print(roc_curve_test)

    roc_auc <- roc_curve_test$auc

    cat("ROC AUC \t\t", roc_auc, "\n\n", sep="")

    predicted_values_binary <- as.numeric(predicted_values)

    predicted_values_binary[predicted_values_binary>=threshold]=1

    predicted_values_binary[predicted_values_binary<threshold]=0

    actual <- actual_labels

    predicted <- predicted_values_binary

  TP <- sum(actual == 1 & predicted == 1)

  TN <- sum(actual == 0 & predicted == 0)

  FP <- sum(actual == 0 & predicted == 1)

  FN <- sum(actual == 1 & predicted == 0)

  cat("\nTOTAL:\n\n")

  cat(" FN = ", (FN), " / ", (FN+TP), "\t (truth == 1) & (prediction < threshold)\n");

  cat(" TP = ", (TP), " / ", (FN+TP),"\t (truth == 1) & (prediction >= threshold)\n\n");

  cat(" FP = ", (FP), " / ", (FP+TN), "\t (truth == 0) & (prediction >= threshold)\n");

  cat(" TN = ", (TN), " / ", (FP+TN), "\t (truth == 0) & (prediction < threshold)\n\n");

  sum1 <- TP+FP; sum2 <-TP+FN ; sum3 <-TN+FP ; sum4 <- TN+FN;

  denom <- as.double(sum1)*sum2*sum3*sum4 # as.double to avoid overflow error on large products

  if (any(sum1==0, sum2==0, sum3==0, sum4==0)) {

    denom <- 1

  }

  mcc <- ((TP*TN)-(FP*FN)) / sqrt(denom)

  f1_score <- 2*TP / (2*TP + FP + FN)

  accuracy <- (TN+TP) / (TN + TP + FP + FN)

  recall <- TP / (TP + FN)

  specificity <- TN / (TN + FP)

  cat("\n\n",keyword,"\t MCC \t F1_score \t accuracy \t TP_rate \t TN_rate \t PR AUC \t ROC AUC\n")

  cat(keyword,"      ", signed_dec_two(mcc), " \t ", dec_two(f1_score), " \t ", dec_two(accuracy), " \t ", dec_two(recall), " \t ", dec_two(specificity),  "\t\t ", dec_two(prc_auc), "\t\t", dec_two(roc_auc),  "\n\n")

#   cat("\nMCC = ", dec_two(mcc), "\n\n", sep="")

#   

#   cat("f1_score = ", dec_two(f1_score), "\n", sep="")

#   cat("accuracy = ", dec_two(accuracy), "\n", sep="")

#   

#   cat("\n")

#   cat("true positive rate = recall = ", dec_two(recall), "\n", sep="")

#   cat("true negative rate = specificity = ", dec_two(specificity), "\n", sep="")

#   cat("\n")

}

# Matthews correlation coefficient

mcc <- function (actual, predicted)

{

  # Compute the Matthews correlation coefficient (MCC) score

  # Jeff Hebert 9/1/2016

  # Geoffrey Anderson 10/14/2016 

  # Added zero denominator handling.

  # Avoided overflow error on large-ish products in denominator.

  #

  # actual = vector of true outcomes, 1 = Positive, 0 = Negative

  # predicted = vector of predicted outcomes, 1 = Positive, 0 = Negative

  # function returns MCC

  TP <- sum(actual == 1 & predicted == 1)

  TN <- sum(actual == 0 & predicted == 0)

  FP <- sum(actual == 0 & predicted == 1)

  FN <- sum(actual == 1 & predicted == 0)

  #TP;TN;FP;FN # for debugging

  sum1 <- TP+FP; sum2 <-TP+FN ; sum3 <-TN+FP ; sum4 <- TN+FN;

  denom <- as.double(sum1)*sum2*sum3*sum4 # as.double to avoid overflow error on large products

  if (any(sum1==0, sum2==0, sum3==0, sum4==0)) {

    denom <- 1

  }

  mcc <- ((TP*TN)-(FP*FN)) / sqrt(denom)

  cat("\nMCC = ", dec_two(mcc), "\n\n", sep="")

  return(mcc)

}