# This is a Shiny application. You can run the application by clicking
# the 'Run App' button above.
#
#    http://shiny.rstudio.com/

# Load shiny library
library(shiny)
library(DT)

# UI
ui <- fluidPage(
  titlePanel("Bayesian probability calculator"),
  sidebarLayout(
    sidebarPanel(
      fileInput("file", "CSV File:",
                multiple = FALSE,
                accept = c("text/csv",
                           "text/comma-separated-values,text/plain",
                           ".csv")),
      #numericInput("nv", "Number of variants:", value = 1, min = 1),
      
      # New input for fraction (Numerator and Denominator)
      tagList(
        h4("Disease prevalence (rate):"),
        fluidRow(
          column(6, textInput("fe_numerator", "Numerator:", value = 0)),
          column(6, textInput("fe_denominator", "Denominator:", value = 0))
        ),
        actionButton("update_fe", "Update Frequency")
      ),
      
      # Numeric Input for Frequency of Disease
      numericInput("fe", "Disease prevalence (decimal):", value = 0.1, min = 0, max = 1, step = 0.01),
      
      # Radio buttons for probability type
      radioButtons("tipo_prob", "Probability:",
                   choices = list("MMC | Exposed" = 1,
                                  "MMC | Non-exposed" = 2,
                                  "Both" = 3),
                   selected = 1),
      
      # Buttons for calculating probabilities and downloading results
      actionButton("calculate", "Calculate Probabilities"),
      downloadButton("download_csv", "Download CSV"),
      downloadButton("download_txt", "Download TXT"),
      
      # Variants input (this is dynamic based on the number of variants)
      uiOutput("variant_inputs")
    ),
    
    mainPanel(
      h3("Bayes' Theorem Framework:"),
      withMathJax(
        p("This calculator computes the probability of Myelomeningocele (MMC) given exposure to variant alleles transmitted by the Transmission Disequilibrium Test (TDT) in family quartets."),
        h4("For individual variants:"),
        p("$$P(\\text{MMC} \\mid \\text{Exposed to TDT transmitted Variant } i) = \\frac{f_{\\text{MMC-sibs}}^{(i)} \\times P(\\text{MMC})}{f_{\\text{MMC-sibs}}^{(i)} \\times P(\\text{MMC}) + f_{\\text{Non-affected-sibs}}^{(i)} \\times (1 - P(\\text{MMC}))}$$"),
        
        # Explanatory text block, aligned left
        p("Where:"),
        
        # Using sprintf for LaTeX formatted text
        HTML(sprintf("\\(f_{\\text{MMC-sibs}}^{(i)}\\) frequency of the TDT transmitted alleles in MMC-affected sibs.")),
        p(),
        HTML(sprintf("\\(f_{\\text{Non-affected-sibs}}^{(i)}\\) frequency of the TDT transmitted alleles in Non-afffected sibs.")),
        p(),
        HTML(sprintf("\\(P(\\text{MMC})\\) Prevalence of MMC.")),
        p(),
        HTML(sprintf("\\(1-P(\\text{MMC})\\) Probability of not having MMC.")),
        p(),
        p(),
        h4("For variants combinations:"),
        p("$$P(\\text{MMC} \\mid C) = \\frac{\\prod_{i=1}^{n} f_{\\text{MMC-sibs}}(V_i) \\times P(\\text{MMC})}{\\prod_{i=1}^{n} f_{\\text{MMC-sibs}}(V_i) \\times P(\\text{MMC}) + \\prod_{i=1}^{n} f_{\\text{Non-affected-sibs}}(V_i) \\times (1 - P(\\text{MMC}))}$$"),
        
        # Explanatory text block
        p("Where:"),
        HTML(sprintf("\\(f_{\\text{MMC-sibs}}(V_i)\\) frequency of the TDT transmitted alleles for variant \\(V_i\\) in MMC-affected sibs.")),
        p(),
        HTML(sprintf("\\(f_{\\text{Non-affected-sibs}}(V_i)\\) frequency of the TDT transmitted alleles for variant \\(V_i\\) in Non-afffected sibs.")),
        p(),
        HTML(sprintf("\\(P(\\text{MMC})\\) prevalence of MMC.")),
        p(),
        HTML(sprintf("\\(C\\) combination of variants \\(\\{V_1, V_2, \\dots, V_n\\}\\) included.")),
        p()
      ),
      h4("RESULTS TABLE"),
      p("Note 1: Ratio = (1/Probability), represents the reciprocal of the calculated probability."),
      dataTableOutput("RESULTS")
    )
  )
)


# Server
server <- function(input, output, session) {
  # Reactive values
  rv <- reactiveValues(variants = NULL, results = NULL, fe_value = NULL)  # Add fe_value to store the updated prevalence
  
  # Observe when the user uploads a CSV file
  observeEvent(input$file, {
    req(input$file)
    data <- read.csv(input$file$datapath)
    nv <- nrow(data)  # Automatically set the number of variants based on the file
    rv$nv <- nv  # Store the number of variants in a reactive value
    
    output$variant_inputs <- renderUI({
      input_list <- lapply(1:rv$nv, function(i) {  # Use the internally determined number of variants
        tagList(
          h4(paste0("Variant ", i)),
          textInput(paste0("name_", i), "Name of the variant:", value = data[i, "Name.of.the.variant"]),
          numericInput(paste0("fee_", i), "Alleles in MMC-sibs | Exposed:", value = data[i, "Transmitted.Alleles.in.MMC.Affected.Exposed"]),
          numericInput(paste0("fne_", i), "Alleles in MMC-sibs | Non-exposed:", value = data[i, "Transmitted.Alleles.in.MMC.Affected.Non.exposed"]),
          numericInput(paste0("fnee_", i), "Alleles in Non affected Sibs | Exposed:", value = data[i, "Transmitted.Alleles.in.Non.affected.Sibs.Exposed"]),
          numericInput(paste0("fnne_", i), "Alleles in Non affected Sibs | Non-exposed:", value = data[i, "Transmitted.Alleles.in.Non.affected.Sibs.Non.exposed"])
        )
      })
      do.call(tagList, input_list)
    })
  })
  
  # Convert fraction input (numerator/denominator) into a decimal and update the numericInput for frequency
  observeEvent(input$update_fe, {
    # Convert the input fraction to a decimal
    fe_numerator <- as.numeric(input$fe_numerator)
    fe_denominator <- as.numeric(input$fe_denominator)
    
    if (fe_denominator == 0) {
      showNotification("Denominator cannot be zero", type = "error")
      return(NULL)
    }
    
    # Calculate the frequency as a fraction
    fe_value <- fe_numerator / fe_denominator
    
    # Store the updated prevalence in reactiveValues
    rv$fe_value <- fe_value  # Store the new prevalence for later use
    
    # Update the numericInput 'fe' with the calculated value
    updateNumericInput(session, "fe", value = fe_value)
  })
  
  # When the calculate button is pressed
  observeEvent(input$calculate, {
    # Create a progress bar
    progress <- shiny::Progress$new(session)
    progress$set(message = "Calculating probabilities", value = 0)
    
    # Get the prevalence from reactiveValues, fall back to the numeric input if it's NULL
    fe <- rv$fe_value
    if (is.null(fe)) {
      fe <- input$fe  # Fallback to input$fe if the reactive value is not set
    }
    
    # Get variant inputs
    nv <- rv$nv  # Use the internally determined number of variants
    vars <- lapply(1:nv, function(i) {
      list(
        name = input[[paste0("name_", i)]],
        fee = input[[paste0("fee_", i)]],
        fne = input[[paste0("fne_", i)]],
        fnee = input[[paste0("fnee_", i)]],
        fnne = input[[paste0("fnne_", i)]]
      )
    })
    
    # Store variants
    rv$variants <- vars
    
    # Calculate the total number of combinations
    total_combinations <- sum(sapply(1:length(vars), function(x) choose(length(vars), x)))
    
    # Initialize the counter for completed combinations
    completed_combinations <- 0
    
    # Calculate and store probabilities
    rv$results <- calculate_probabilities(vars, fe, input$tipo_prob, progress, total_combinations)
    
    # Close the progress bar
    progress$close()
  })
  
  # Display the results in a data table
  output$RESULTS <- renderDataTable({
    if (!is.null(rv$results)) {
      rv$results
    } else {
      data.frame(Result = "No results to display.")
    }
  })
  
  # Download results as TXT
  output$download_txt <- downloadHandler(
    filename = function() {
      paste("results", ".txt", sep = "")
    },
    content = function(file) {
      if (!is.null(rv$results)) {
        write.table(rv$results, file, sep = "\t")
      }
    }
  )
  
  # Download results as CSV
  output$download_csv <- downloadHandler(
    filename = "results.csv",
    content = function(file) {
      if (!is.null(rv$results)) {
        write.csv(rv$results, file)
      }
    }
  )
}


# Calculate probabilities function
calculate_probabilities <- function(vars, fe, tipo_prob, progress, total_combinations) {
  
  # Initialize a data frame to store the results
  results <- data.frame(
    Names = character(),
    Probability = numeric(),
    Ratio = character(),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  
  # Functions to calculate probabilities
  get_prob_comb_enf_exp <- function(v, fe) {
    fne <- 1 - fe
    num <- prod(sapply(v, function(var) var$fee)) * fe
    denom <- num + prod(sapply(v, function(var) var$fnee)) * fne
    num / denom
  }
  
  get_prob_comb_enf_no_exp <- function(v, fe) {
    fne <- 1 - fe
    num <- prod(sapply(v, function(var) var$fne)) * fe
    denom <- num + prod(sapply(v, function(var) var$fnne)) * fne
    num / denom
  }
  
  total_combinations <- 2^length(vars) - 1
  completed_combinations <- 0
  
  for (n_comb in 1:length(vars)) {
    comb_indices <- combn(length(vars), n_comb)
    
    for (i in 1:ncol(comb_indices)) {
      comb <- comb_indices[, i]
      selected_vars <- lapply(comb, function(index) vars[[index]])
      combination_name <- paste0(sapply(selected_vars, function(var) var$name), collapse = " & ")
      
      if (tipo_prob == 1) {
        # Only calculate P(MMC | Exposed)
        prob_enf_exp <- get_prob_comb_enf_exp(selected_vars, fe)
        Ratio_exp <- ifelse(prob_enf_exp > 0, paste0("(1/", round(1/prob_enf_exp), ")"), "...")
        
        results <- rbind(results, data.frame(
          Names = paste0("P(MMC | Exposed) ", combination_name),
          Probability = prob_enf_exp,
          Ratio = Ratio_exp,
          stringsAsFactors = FALSE
        ))
      } else if (tipo_prob == 2) {
        # Only calculate P(MMC | Non-exposed)
        prob_enf_no_exp <- get_prob_comb_enf_no_exp(selected_vars, fe)
        Ratio_no_exp <- ifelse(prob_enf_no_exp > 0, paste0("(1/", round(1/prob_enf_no_exp), ")"), "...")
        
        results <- rbind(results, data.frame(
          Names = paste0("P(MMC | Non-exposed) ", combination_name),
          Probability = prob_enf_no_exp,
          Ratio = Ratio_no_exp,
          stringsAsFactors = FALSE
        ))
      } else if (tipo_prob == 3) {
        # Calculate both probabilities
        prob_enf_exp <- get_prob_comb_enf_exp(selected_vars, fe)
        prob_enf_no_exp <- get_prob_comb_enf_no_exp(selected_vars, fe)
        or_val <- prob_enf_exp / fe
        
        # Handling the Ratio format for exposed and non-exposed
        Ratio_exp <- ifelse(prob_enf_exp > 0, paste0("(1/", round(1/prob_enf_exp), ")"), "...")
        Ratio_no_exp <- ifelse(prob_enf_no_exp > 0, paste0("(1/", round(1/prob_enf_no_exp), ")"), "...")
        
        results <- rbind(results,
                         data.frame(
                           Names = paste0("P(MMC | Exposed) ", combination_name),
                           Probability = prob_enf_exp,
                           Ratio = Ratio_exp,
                           stringsAsFactors = FALSE
                         ),
                         data.frame(
                           Names = paste0("P(MMC | Non-exposed) ", combination_name),
                           Probability = prob_enf_no_exp,
                           Ratio = Ratio_no_exp,
                           stringsAsFactors = FALSE
                         ),
                         data.frame(
                           Names = paste0("Ratio [P(MMC | Exposed)/ P(MMC)] ", combination_name),
                           Probability = or_val,
                           Ratio = "-",
                           stringsAsFactors = FALSE
                         )
        )
      }
      
      completed_combinations <- completed_combinations + 1
      progress$set(value = completed_combinations / total_combinations)
    }
  }
  
  return(results)
}


# Run the app
shinyApp(ui = ui, server = server)