# Load necessary library
#Original code by Valentina Tillman
library(mtbr)

# Define Parameters
AQL <- mtb_get_column("AQL")  # Acceptable Quality Level
alpha <- mtb_get_column("alpha")  # Producer's risk
RQL <- mtb_get_column("RQL")  # Rejectable Quality Level
beta <- mtb_get_column("beta")  # Consumer's risk

# ---- Validate inputs ----

# Validate AQL and RQL
if (is.null(AQL) || length(AQL) != 1 || is.na(AQL)) {
  warning("Define the AQL with a value in the first row of a column with the title AQL.")
}
if (is.null(RQL) || length(RQL) != 1 || is.na(RQL)) {
  warning("Define the RQL with a value in the first row of a column with the title RQL.")
}

# Validate alpha
if (is.null(alpha) || length(alpha) != 1 || is.na(alpha)) {
  alpha <- 0.05
  message("Missing or multiple values for alpha. The analysis uses 0.05. To specify alpha, enter a value in the first row of a column with the title alpha.")
}

# Validate beta
if (is.null(beta) || length(beta) != 1 || is.na(beta)) {
  beta <- 0.10
  message("Missing or multiple values for beta. The analysis uses 0.10. To specify beta, enter a value in the first row of a column with the title beta.")
}

#Simple comparison of AQL and RQL
if (AQL >= RQL) {
  stop("Specify an RQL that is greater than the AQL.")
}

# ---- Calculations ----
# Use algorithm to find n,c for the Poisson distribution
found <- FALSE
for (n in 10:5000) {
  lamA <- n*AQL
  lamR <- n*RQL
  for (c in 0:n) {
    p_accept_AQL <- ppois(c, lamA)
    p_accept_RQL <- ppois(c, lamR)
    if (p_accept_AQL >= 1 - alpha && p_accept_RQL <= beta) {
      found <- TRUE
      break
    }
  }
  if (found) break
}

if (found != TRUE) { 
  mtb_add_message ("No solution found. Results are for the largest values of n and c in the search.")
  }

# Generate OC curve values for number of defects per unit
p_vals <- seq(0.003, 0.117, by = 0.006)
lambda_vals <- p_vals * n
OC_vals <- ppois(c, lambda_vals)

# ---- Outputs ----
#Display analysis inputs
show_AQL <- sprintf("AQL: %s", sub("\\.?0+$", "", formatC(AQL, format = "f", digits = 6)))
mtb_add_message(show_AQL)
show_RQL <- sprintf("RQL: %s", sub("\\.?0+$", "", formatC(RQL, format = "f", digits = 6)))
mtb_add_message(show_RQL)
show_alpha <- sprintf("alpha: %s", sub("\\.?0+$", "", formatC(alpha, format = "f", digits = 6)))
mtb_add_message(show_alpha)
show_beta <- sprintf("beta: %s", sub("\\.?0+$", "", formatC(beta, format = "f", digits = 6)))
mtb_add_message(show_beta)

#Display calculated solutions
show_sample_size <- sprintf("Sample size: %d", n)
mtb_add_message(show_sample_size)
show_defects <- sprintf("Acceptance number: %d", c)
mtb_add_message(show_defects)

# OC Curve Table
oc_table <- data.frame("Lot_Defects_per_Unit" = p_vals, "Probability_of_Acceptance" = OC_vals)
head(oc_table, 10)

mytitle <- "OC Table" 
myheaders <- names(oc_table)
mtb_add_table(columns = oc_table, headers = myheaders, title = mytitle)