In CDU-data-science-team/pharmacyReporting: Generate Reports For Pharmacy Inventory Project

library(testthat)

# Load already included functions
pkgload::load_all(export_all = FALSE)

My function

#' Return the amount of a drug to be ordered
#' 
#' @param forecast dataframe with two columns- mean_demand and sd_demand
#' @param distribution distr6 object defining the distribution of order arrival
#' @param min_stock number level below which we want a probability < p_min of falling 
#' @param max_stock number level above which we want a prob < p_max of rising
#' @param p_min number probability of going into emergency stock
#' @param p_max number probability of ordering too much
#' @param inv_i int the inventory at the time of this (the i'th) order
#' @param delta_pref int the preferred order interval for drug in question
#' @param outstanding_orders int the number of units in outstanding orders
#' @return number giving amount to order
#' @export
drug_quantity <- function(forecast, distribution, min_stock, max_stock,
                          p_min, p_max, inv_i, delta_pref,
                          outstanding_orders) {

  MaxLoops = 100  # maximum number of loops permitted in each while loop below
  Q_i = 0 # initialise order quantity
  Delta_i = delta_pref # initialise time to next order
  flag_suf = 0 # initialise a flag indicating (=1) that Q_i is sufficent 
  flag_stor = 0 # initialise a flag indicating (=1) that Q_i is not too much 
  Forecast_quantiles <- make_quantiles(forecast)

  tmp1 <- c(0,0) # dummy storage for Q_i and flag_suf 
  dummy_counter <- 0

  while(flag_suf < 1 & dummy_counter < MaxLoops){ # repeat until Q_i sufficient

    dummy_counter <- dummy_counter + 1

    # returns whether Q_i sufficient and next Q_i to try

    tmp1 <- Q_enough_Q(
      distribution, inv_i, Q_i, outstanding_orders, Forecast_quantiles, 
      Delta_i, min_stock, p_min)

    Q_i = tmp1[1]

    flag_suf = tmp1[2]

  }

  # ERROR HANDLING ON POSSIBLE TIME-OUT ABOVE

  tmp2 <- c(Q_i, flag_stor)

  # returns whether Q_i satisfies storage condition and next Q_i to try
  tmp2 <- Q_toomuch_Q(forecast_q = Forecast_quantiles, 
                      o_orders = outstanding_orders,
                      choose_distribution = distribution,
                      current_q_i = Q_i,
                      max_stock = max_stock,
                      p_max = p_max,
                      inv_i = inv_i)

  Q_i <- tmp2[1]
  flag_stor <- tmp2[2]

  # if sufficient Q_i breaches storage, need to reduce and shorten Delta

  if(flag_stor < 1){

    flag_suf <- 0    # reset sufficiency flag to zero as about to reduce Q_i

    dummy_counter_2 <- 0
    # look for Q_i that meets storage condition 

    while(flag_stor<1 & dummy_counter_2 < MaxLoops){
      dummy_counter_2 <- dummy_counter_2 + 1

      # returns whether Q_i satisfies storage condition and next Q_i to try

      tmp2 <- Q_toomuch_Q(forecast_q = Forecast_quantiles, 
                          o_orders = outstanding_orders,
                          choose_distribution = distribution,
                          current_q_i = Q_i,
                          max_stock = max_stock,
                          p_max = p_max,
                          inv_i = inv_i) 

      Q_i <- tmp2[1]
      flag_stor <- tmp2[2]

    }

    tmp3 <- c(0,0)   # temporary storage of Delta_i and flag_suf
    dummy_counter_3 <- 0

    # now reduce Delta_i until reduced Q_i sufficient     

    while(flag_suf < 1 & dummy_counter_3 < MaxLoops){
      dummy_counter_3 <- dummy_counter_3 + 1
      tmp3 <- Q_enough_Delta(Forecast_quantiles,
                             choose_distribution = distribution,
                             d_i = Delta_i,
                             inv_i = inv_i,
                             current_q_i = Q_i,
                             min_stock, 
                             p_min, 
                             outstanding_orders = outstanding_orders)

      Delta_i <- tmp3[1]
      flag_suf <- tmp3[2]

    }
  }
  return(list("Q_i" = Q_i, 
              "Delta_i" = Delta_i))
}

# set up as piece wise linear from input.

pwlcdf <- function(Forecast_quantiles, q_vals, num_q_vals, time_point, x){ 

  res <- numeric(length(x))

  for(i in seq(1:length(x))){

    low_cut <- Forecast_quantiles[time_point,1]
    high_cut <- Forecast_quantiles[time_point, num_q_vals]

    res[i][x[i] >= high_cut] = 1

    for(j in seq(1, num_q_vals - 1, 1)) {
      res[i][x[i] >= Forecast_quantiles[time_point,j] & 
               x[i] < Forecast_quantiles[time_point,j+1]] = 
        q_vals[j] + (q_vals[j+1] - q_vals[j]) * (x[i] - Forecast_quantiles[time_point, j]) / 
        (Forecast_quantiles[time_point,j+1] - Forecast_quantiles[time_point, j])

    }

    res[i][x[i] <= low_cut] = 0

  }
  return(res)
}

# set up as piece wise linear from input.

pwlquant <- function(Forecast_quantiles, q_vals, num_q_vals, time_point, p){ 

  res <- numeric(length(p))

  for (i in seq(1, length(p))){

    for(j in seq(1, num_q_vals - 1, 1)) {

      res[i][p[i] >= q_vals[j] & p[i] < q_vals[j + 1]] = Forecast_quantiles[time_point, j] + 
        (Forecast_quantiles[time_point, j + 1] - 
           Forecast_quantiles[time_point, j])*(p[i] - q_vals[j]) / (q_vals[j + 1] - q_vals[j])

    }

    res[i][p[i] >= q_vals[num_q_vals]] = Forecast_quantiles[time_point, num_q_vals]
  }
  return(res)
}

Q_enough_Delta <- function(forecast_q, choose_distribution, d_i, inv_i, current_q_i,
                           min_stock, p_min, outstanding_orders){

  # determines whether Q_i is sufficient and, if not, suggests next Delta_i to try 
  # function used when Q_i has been reduced to meet storage condition

  epsilon <- 0.01
  flag1 <- 0
  lmax <- length(choose_distribution)
  Delta_i <- d_i
  # prob that, if next  order delivered at time y, Q_i for this 
  # order will be insufficient to meet demand to then
  P_Q_insuff <- c(rep(1,Delta_i+lmax+1)) 
  # initialise to zero - probability that the next order is delivered at time y
  Prob_y <-c(rep(0,Delta_i+lmax+1)) 
  term_y <- c(rep(0,Delta_i+lmax+1)) 
  Q_i <- current_q_i
  num_q_vals <- ncol(forecast_q)
  Forecast_quantiles <- forecast_q

  # set up vector of quantile levels - evenly spaced for now
  q_vals <- c(seq(0,1, 1 / (num_q_vals - 1))) 
  Q_out <- sum(outstanding_orders)

  # now set probabilities for next delivery arriving at time y 
  # and Q_i being insufficient if that is the case

  for (y in 1:Delta_i){
    Prob_y[y] = 0

    # probability that demand up to and including day y eats into emergency stock
    P_Q_insuff[y] <- 1 - pwlcdf(Forecast_quantiles, q_vals, num_q_vals, y, 
                                inv_i + Q_out + Q_i - min_stock) 
  }

  # we only need to go just beyond (Delta_i + longest lead time)

  for (y in seq(Delta_i + 1, Delta_i + lmax, 1)) { 

    Prob_y[y] = choose_distribution[y-Delta_i]

    # probability that demand up to and including day y eats into emergency stock

    P_Q_insuff[y] <- 1 - pwlcdf(Forecast_quantiles, q_vals, num_q_vals, y, 
                                inv_i + Q_out + Q_i - min_stock) 

  } # loop over second set of values of y

  phi <- sum(term_y)

  if (phi > p_min) { # Q_i insufficient given current value of Delta_i

    sc <-  (p_min / phi) * (1 - epsilon) # get reduction required in phi
    y_peak <- which.max(term_y) # find biggest term in phi

    # find target probability for next delivery being ordered at y

    P_target <- Prob_y[y_peak] * sc  

    y_hat <- y_peak 
    test <- Prob_y[y_hat]

    # effectively we are bringing the order of the next forward until we get 
    # desired reduction in the (currently) biggest contribution 
    # to the prob of Q_i being insufficient 

    while(test > P_target){   

      y_hat <- y_hat + 1  
      test <- Prob_y[y_hat]
    }

    Delta_i <- Delta_i -(y_hat - y_peak)

  } else {

    flag1 = 1        
  }
  if (Delta_i < 1) {
    print("CANNOT BE SOLVED WITH CURRENT INPUTS - REVISE STORAGE CONSTRAINT, 
                          EMERGENCY STOCK OR TOLERANCES")
  } 
  res <- c(Delta_i, flag1)

  return(res)
}

Q_enough_Q <- function(lead_time_dis, inv_i, current_q_i, outstanding_orders, 
                       Forecast_quantiles, Delta_i, min_stock, p_min){

  # This function assesses whether Q_i will be sufficient and, if not, 
  # suggests next Q_i to try  

  epsilon <- 0.01
  flag1 = 0 # initialise as insufficient
  phi = 1 # this is the probability that order is insufficient

  lmax <- length(lead_time_dis)

  # this is the prob that, if next  order delivered at time y, Q_i for this 
  # order will be insufficient to meet demand to then

  P_Q_insuff <- c(rep(1,Delta_i + lmax + 1)) 
  # initialise to zero - this is the probability that the next order is 
  # delivered at time y

  Prob_y <-c(rep(0,Delta_i+lmax+1)) 
  term_y <- c(rep(0,Delta_i+lmax+1))

  num_q_vals <- ncol(Forecast_quantiles)
  # set up vector of quantile levels - evenly spaced for now
  q_vals <- c(seq(0,1, 1 / (num_q_vals - 1))) 

  Q_i <- current_q_i
  Q_out <- sum(outstanding_orders) # quantity associated with outstanding orders

  for (y in 1 : Delta_i){

    Prob_y[y] = 0

    # probability that demand up to and including day y eats into emergency stock
    P_Q_insuff[y] <- 1 - pwlcdf(Forecast_quantiles,
                                q_vals, num_q_vals, y,
                                inv_i + Q_out + Q_i - min_stock) 
  }  
  # we only need to go just beyond (Delta_i + longest lead time)
  for (y in seq(Delta_i + 1, Delta_i + lmax, 1)) {

    Prob_y[y] = lead_time_dis[y - Delta_i] # wrong code Martin sent
    # Prob_y[y] = lead_time_dis$cdf(y-Delta_i)

    # returns probability that demand up to and including day y eats into emergency stock
    P_Q_insuff[y] <- 1 - pwlcdf(Forecast_quantiles, q_vals, num_q_vals, y, 
                                inv_i + Q_out + Q_i - min_stock) 

  } # loop over second set of values of y

  # joint probability delivery of next order occurs at y and that q_i insufficient if so

  term_y <- Prob_y * P_Q_insuff  

  # sums over all possible delivery times y to give probability that Q_i insufficient

  phi <- sum(term_y) 

  if (phi > p_min) { # if Q_i not sufficient

    # get scaling factor required to bring prob phi under tolerance value - 
    # with slight over-adjustment to stop any asymptotic behaviour

    sc <-  (p_min / phi) * (1 - epsilon) 
    y_peak <- which.max(term_y)  # find biggest term in phi

    # get target for reduced contribution from biggest term

    P_target <- P_Q_insuff[y_peak] * sc 

    # get from forecast the demand associated with that target probability

    B <- pwlquant(Forecast_quantiles, q_vals, num_q_vals, y_peak, (1 - P_target))

    # set Q_i to be sufficient reduce bring biggest term in phi by sc - 
    # amount determined above

    Q_i <- B + min_stock - inv_i - Q_out

  } else {

    flag1 = 1        

  }

  res <- c(Q_i,flag1)

  return(res)
}

Q_toomuch_Q <- function(forecast_q, o_orders, choose_distribution, current_q_i,
                        max_stock, p_max, inv_i) {

  # function indicates if Q_i meets storage condition and, if not, gives next Q_i to try  

  epsilon <- 0.01
  flag2 = 0
  Q_i <- current_q_i
  Forecast_quantiles <- forecast_q

  num_q_vals <- ncol(Forecast_quantiles)
  # set up vector of quantile levels - evenly spaced for now
  q_vals <- c(seq(0,1, 1 / (num_q_vals - 1))) 


  # NEW FOR OCTOBER 2021  

  lmax <- length(choose_distribution)   

  # p if next order delivered at time y
  # Q_i for this order will be insufficient to meet demand to then

  P_Q_toomuch <- c(rep(1,lmax)) 

  # initialise to zero - p that the next order is delivered at time y
  Prob_x <-c(rep(0,lmax)) 
  term_x <- c(rep(0,lmax))  

  #  P_Q_toomuch <- c(rep(1,Forecast_length)) # initialise as too much
  #  Prob_x <-c(rep(0,Forecast_length)) # initialise to zero
  #  term_x <- c(rep(0,Forecast_length))  

  Q_out <- o_orders

  # NEW FOR OCTOBER 2021  

  Prob_x <- choose_distribution
  #  for (x in seq(1,Forecast_length,1)){

  for (x in 1 : lmax){

    # to reflect change in lead_time_dis from a Distr6 entity to a vector           
    #    Prob_x[x] <- lead_time_dis$cdf(x)- lead_time_dis$cdf(x-1)  # probability that this delivery will be on day x
    #    P_Q_toomuch[x] <-  pnorm((Inv_i+Q_out+Q_i-Storage_constraint),
    #      Forecast$mean_cumulative_demand[x],Forecast$sd_cumulative[x], lower.tail = TRUE)    
    #    P_Q_toomuch[x] <- Fcast[x]$cdf(Inv_i+Q_out+Q_i-Storage_constraint,lower.tail = TRUE )   

    # gives probability that, if this delivery arrives at time x, it will 
    # breach storage constraint
    P_Q_toomuch[x] <- pwlcdf(
      Forecast_quantiles, q_vals, num_q_vals, x, 
      inv_i + Q_out + Q_i - max_stock)
  }

  term_x <- Prob_x * P_Q_toomuch  # joint probability that delivery occurs at x and that Q_i breaches storage if so

  gam <- sum(term_x) # summming over all possible x, this is prob that delivery will breach storage on arrival

  if (gam > p_max) { # compare to tolerance 

    sc <-  (p_max / gam) * (1 - epsilon) # get reduction required in prob of breach
    x_peak <- which.max(term_x) # find x with biggest contribution to prob of breach

    P_target <- P_Q_toomuch[x_peak] * sc   # get target for reduced contribution from biggest term

    #    B <- qnorm(P_target,Forecast$mean_cumulative_demand[x_peak],Forecast$sd_cumulative[x_peak],lower.tail = TRUE)
    #    B <- Fcast[x_peak]$quantile(P_target,lower.tail = TRUE)  
    B <- pwlquant(Forecast_quantiles,q_vals,num_q_vals,x_peak,P_target)   # get from forecast the demand associated with target prob

    # set Q_i to meet that demand-  this is not guaranteed to be enough of a reduction.

    Q_i <- B - inv_i - Q_out + max_stock   
  } else {

    flag2 = 1        

  }

  res <- c(Q_i,flag2)

  return(res)
}

#' Forecast and run inventory management algorithm, returning reorder quantities
#' by site and drug
#'
#' @param site String. Site code (selected from dynamic UI listing all sites)
#' @param supplier String. Supplier (changes weekly, selected in Shiny interface)
#' @param product dataframe. Contents of product_sup_profile, loaded in 
#' app_server.R
#' @param w_order dataframe. Contents of w_order_log_df1, loaded in 
#' app_server.R
#'
#' @return dataframe containing drug name, stock level, days to order, etc.
#' @export
#'
inventory_reorder <- function(site, supplier, product, w_order, requis, holidays,
                              updateProgress = NULL){

  product <- product %>% 
    dplyr::filter(Site == site, Supplier_name == supplier) %>% 
    dplyr::arrange(Drug_name)

  for(chem in 1:nrow(product)[1]){

    starting_stock_level1 <- initial_product_list %>% 
      dplyr::filter(Drug_code == product$Drug_code[chem])

    starting_stock_level <- starting_stock_level1$stocklvl[1]

    w_trans_log_df2a <- w_trans_log_df1 %>% 
      dplyr::filter(Drug_code == product$Drug_code[chem], Site == site, 
                    Date > as.Date("2022-02-25")) %>% 
      dplyr::summarise(Total_Qty = sum(Qty))


    w_order_log_df3 <- w_order_log_df2 %>% 
      dplyr::filter(Drug_code == product$Drug_code[chem] & Site == site & 
                      Kind == "R" & DateReceived > as.Date("2022-02-25") & 
                      DateReceived < Sys.Date()) %>% 
      dplyr::summarise(Rec_Qty = (sum(QtyRec) * product$Packsize[chem]))


    stock_level <- starting_stock_level + 
      w_order_log_df3$Rec_Qty - w_trans_log_df2a
    product$stocklvl[chem] <- as.numeric(stock_level)
    product$stocklvl <- utils::type.convert(product$stocklvl)

    product$stocklvl[chem] <- ifelse(product$stocklvl[chem] < 0, 0, 
                                     product$stocklvl[chem])

    date_last_ordered <- w_order_log_df2 %>% 
      dplyr::filter(Drug_code == product$Drug_code[chem], Site == site,
                    Kind %in% c("D", "O", "I")) %>% 
      dplyr::arrange(desc(DateOrdered))

    product$lastordered[chem] <- date_last_ordered$DateOrdered[1]

  }

  # Settings which aren't yet provided within the code
  max_storage_capacity <-  100000

  # add in other waste / expiry / adjustment codes
  waste_adjust_codes <- c("Adj", "ADJ","COMSP","EXP", "HADJ","HEXP", "HWAST", 
                          "MOCK",  "RWAST", "TEST", "TRG", "WAST", "WASTE", 
                          "XXXX")

  # list of excluded suppliers for site 100

  cross_site_order_exclusions <- c("HIGH", "MILL", "PHWRC", "RAMP")

  # takes two reactive inputs- site and supplier

  order_list <- product %>% 
    dplyr::filter(Site == site, Supplier_name == supplier) %>% 
    dplyr::arrange(Drug_name)

  # Add in order cycle - only needed for shadow testing
  order_list <- dplyr::left_join(order_list, product_ord_cycle,
                                 by = c("Site" = "Site", "Drug_code" = "Drug_code",
                                        "Supplier_name" = "Supplier_name"))

  # Identify new products which haven't been assigned an AAH ordercycle and 
  # add in ordercycle details - only needed for shadow testing
  order_list$Order_cycle[is.na(order_list$Order_cycle)] <- "empty"
  order_list$Order_cycle <- ifelse(order_list$Order_cycle == "empty", "AAHb", 
                                   order_list$Order_cycle)

  # select orders placed within the last 2 years
  order_log <- w_order %>% 
    subset(DateOrdered > Sys.Date() - 730) %>% 
    dplyr::filter(Site == site)

  # this is only needed for site 100
  if(site == 100){
    order_log <- order_log %>%
      dplyr::filter(!Supplier_name %in% cross_site_order_exclusions)
  }

  purrr::pmap_dfr(order_list, function(Drug_code, ProductID, ...){

    comment <-  "None"

    # cat(Drug_code)

    # Filter product for relevant product

    product_info <- product %>% 
      dplyr::filter(.data$Drug_code == .env$Drug_code, Site == site)

    # calculate leadtimes

    ord_log <- order_log %>% 
      dplyr::filter(Kind %in% c("D", "O", "I") & 
                      .data$Drug_code == .env$Drug_code) %>% 
      dplyr::select(OrderNum, Drug_code, DateOrdered, QtyOrd)

    if(nrow(ord_log) <= 1){

      return(
        data.frame(drug = Drug_code,
                   order = NA,
                   actual_order_qty = 0,
                   days_to_order = NA,
                   Drug_Name = product_info$Drug_name[1],
                   Supplier_Tradename = product_info$SupplierTradeName[1],
                   PackSize = product_info$Packsize[1],
                   Units = product_info$PrintForm[1],
                   Stock_Level = product_info$stocklvl[1],
                   Last_Order_Date = product_info$lastordered[1],
                   Min_stock_risk = NA,
                   Min_stock_level = NA,
                   Comments = "No history - order manually")
      )
    }

    # Get rid of multiple lines for same order & received date
    rec_log <- order_log %>% 
      dplyr::filter(Site == site & Kind == "R" & 
                      .data$Drug_code == .env$Drug_code) %>% 
      dplyr::group_by(OrderNum, Drug_code, Supplier_name, 
                      DateOrdered, DateReceived, Site, Kind) %>%  
      dplyr::summarise(QtyRec = sum(QtyRec)) %>% 
      dplyr::ungroup()

    if(nrow(rec_log) <= 1){

      return(
        data.frame(drug = Drug_code,
                   order = NA,
                   actual_order_qty = 0,
                   days_to_order = NA,
                   Drug_Name = product_info$Drug_name[1],
                   Supplier_Tradename = product_info$SupplierTradeName[1],
                   PackSize = product_info$Packsize[1],
                   Units = product_info$PrintForm[1],
                   Stock_Level = product_info$stocklvl[1],
                   Last_Order_Date = product_info$lastordered[1],
                   Min_stock_risk = NA,
                   Min_stock_level = NA,
                   Comments = "No history - order manually")
      )
    }

    ord_rec_log <- dplyr::left_join(ord_log, rec_log, 
                                    by = c("OrderNum" = "OrderNum", 
                                           "Drug_code" = "Drug_code", 
                                           "DateOrdered" = "DateOrdered")) %>% 
      tidyr::replace_na(list(Supplier_name = supplier,
                             DateReceived = as.Date("1999-12-30"),
                             Site = site,
                             Kind = "R",
                             QtyRec = 0)) 

    ord_rec_log <- ord_rec_log %>% 
      dplyr::mutate(Outstanding_order = QtyOrd - QtyRec)

    orders_outstanding <- ord_rec_log %>%
      dplyr::group_by(Site, Supplier_name, Drug_code, OrderNum, DateOrdered, QtyOrd) %>%
      dplyr::summarise(TotalRecQty = sum(QtyRec)) %>%
      dplyr::mutate(Outstanding_order = QtyOrd - TotalRecQty) %>%
      dplyr::ungroup()

    leadtime <- ord_rec_log %>% 
      dplyr::filter(DateReceived != "1999-12-30") %>% 
      dplyr::rowwise() %>%
      dplyr::mutate(leadtime = 
                      n_weekdays(DateOrdered, DateReceived, holidays)) %>% 
      subset(leadtime <= 14)

    if(nrow(leadtime) == 0){

      return(
        data.frame(drug = Drug_code,
                   order = NA,
                   actual_order_qty = 0,
                   days_to_order = NA,
                   Drug_Name = product_info$Drug_name[1],
                   Supplier_Tradename = product_info$SupplierTradeName[1],
                   PackSize = product_info$Packsize[1],
                   Units = product_info$PrintForm[1],
                   Stock_Level = product_info$stocklvl[1],
                   Last_Order_Date = product_info$lastordered[1],
                   Min_stock_risk = NA,
                   Min_stock_level = NA,
                   Comments = "No history - order manually")
      )
    }

    leadmin <- min(leadtime$leadtime)
    leadmax <- max(leadtime$leadtime)
    leadmode <- calc_mode(leadtime$leadtime)

    # set distribution for delivery lead time
    lead_time_distribution <- distr6::Triangular$new(
      lower = leadmin - 0.5, 
      upper = leadmax + 0.5, # Have had to make upper +0.5 rather than -0.5 in  
      # original code otherwise if min, max & mode is 1 distr6 doesn't work 
      # as upper is less than mode
      mode = leadmode)

    lead_time_dis <- lead_time_distribution$pdf(seq(leadmin - 0.5, leadmax + 0.5, 1))

    # find outstanding order quantity for orders placed within last 14 days
    outstanding_order_qty  <- orders_outstanding %>%
      subset(DateOrdered > Sys.Date() - 14) %>%
      dplyr::mutate(out_ord_qty = 
                      Outstanding_order * product_info$Packsize[1]) %>% 
      dplyr::summarise(total_out_ord_qty = sum(out_ord_qty))

    # See if any outstanding orders from last 60 days add comment to shiny table
    outstanding_order_comment <- orders_outstanding %>%
      subset(DateOrdered > Sys.Date() - 60) %>%
      dplyr::filter(Outstanding_order > 0) %>%
      dplyr::arrange(desc(DateOrdered))

    if(nrow(outstanding_order_comment) == 0){
      comment <- "None"
    }else{
      comment <- paste("Outstanding order from", 
                       outstanding_order_comment$DateOrdered[1], sep = " ")
    }

    # find outstanding requisitions which need fulfilling

    outstanding_requis <- requis %>%
      dplyr::filter(Site == site, .data$Drug_code == .env$Drug_code)

    if(nrow(outstanding_requis) == 0){
      outstanding_requis_quantity <- 0L

    }else{
      outstanding_requis <- outstanding_requis %>%
        subset(DateOrdered > Sys.Date() - 14) %>%
        dplyr::summarise(total_out_requis =
                           sum(Outstanding))
      outstanding_requis_quantity <- outstanding_requis$total_out_requis[1]
      comment <- "Outstanding requisition"
    }

    # Filter out issues to adjustment and waste codes
    # to consider in future deleting issues done in error i.e. 
    # where issue is followed by a return to stock

    demand_data <- w_trans_log_df1 %>% 
      dplyr::filter(Site == site,
                    .data$Drug_code == .env$Drug_code,
                    !Ward %in% waste_adjust_codes,
                    Qty >= 0,
                    Date < order_date,
                    Date < Sys.Date()) %>%
      rbind(data.frame(WTranslogID = "NA", Date = Sys.Date(), 
                       Drug_code = Drug_code,  
                       Qty = 0, Ward = "NA", Site = site)) %>% 
      dplyr::arrange(Date) %>% 
      dplyr::group_by(Date) %>% 
      dplyr::summarise(Total_Qty = sum(Qty)) %>% 
      dplyr::ungroup() %>% 
      tidyr::complete(Date = 
                        seq.Date(min(Date), 
                                 max(Date), by="day")) %>% 
      tidyr::replace_na(list(Total_Qty = 0))

    # if no orders in this range jump out of function

    if(sum(demand_data$Total_Qty) == 0){

      return(
        data.frame(drug = Drug_code,
                   order = NA,
                   actual_order_qty = 0,
                   days_to_order = NA,
                   Drug_Name = product_info$Drug_name[1],
                   Supplier_Tradename = product_info$SupplierTradeName[1],
                   PackSize = product_info$Packsize[1],
                   Units = product_info$PrintForm[1],
                   Stock_Level = product_info$stocklvl[1],
                   Last_Order_Date = product_info$lastordered[1],
                   Min_stock_risk = NA,
                   Min_stock_level = NA,
                   Comments = "No history - order manually")
      )
    }

    # Calculate min_stock_level
    min_stock_level <- demand_data %>% 
      subset(Date > Sys.Date() - 365)

    min_stock_level <- stats::quantile(min_stock_level$Total_Qty, 
                                       .98, na.rm = TRUE) %>% 
      ceiling()

    # Work out p_min from order scheduling

    product_ordercycle <- product_ord_cycle %>%
      dplyr::filter(.data$Drug_code == .env$Drug_code)

    sup_ordercycle <- sup_order_sched %>%
      dplyr::filter(Supplier == supplier & Order_cycle == product_ordercycle$Order_cycle[1])

    ar <- (lubridate::interval(sup_ordercycle$Order_cycle_start_date[1], Sys.Date()) %>%
             as.numeric('days')) %% 14

    time_til_next_order <- ord_sch %>%
      dplyr::filter(Day == sup_order_sched$Order_day[1] & x == ar)

    if(nrow(time_til_next_order) == 0){

      return(
        data.frame(drug = Drug_code,
                   order = NA,
                   actual_order_qty = 0,
                   days_to_order = NA,
                   Drug_Name = product_info$Drug_name[1],
                   Supplier_Tradename = product_info$SupplierTradeName[1],
                   PackSize = product_info$Packsize[1],
                   Units = product_info$PrintForm[1],
                   Stock_Level = product_info$stocklvl[1],
                   Last_Order_Date = product_info$lastordered[1],
                   Min_stock_risk = NA,
                   Min_stock_level = NA,
                   Comments = "No history - order manually")
      )
    }

    if(time_til_next_order$Delta_p[1] > 7){
      risk_of_min_stock <- 0.15
      risk_of_exceeding_max_stock <- 0.05

    }else{
      risk_of_min_stock <- 0.4
      risk_of_exceeding_max_stock <- 0.05
    }

    # Run forecast
    demand_data <- make_tsibble(demand_data, frequency = "Daily")

    daily_forecast <- forecast_series(demand_data, length(lead_time_dis) + 14,
                                      frequency = "Daily")

    actual_forecast <- daily_forecast %>% 
      dplyr::filter(.model == "ARIMA")

    step <- drug_quantity(
      forecast = actual_forecast,
      distribution = lead_time_dis,
      min_stock = min_stock_level,
      max_stock = max_storage_capacity,
      p_min = risk_of_min_stock,
      p_max = risk_of_exceeding_max_stock,
      inv_i = product_info$stocklvl[1],
      delta_pref = time_til_next_order$Delta_p[1],
      outstanding_orders = outstanding_order_qty$total_out_ord_qty[1])

    comment <- ifelse(comment != "None", comment, "None")

    to_return <- data.frame(drug = Drug_code,
                            order = ceiling(step$Q_i / product_info$Packsize[1] +
                                              outstanding_requis_quantity),
                            actual_order_qty = ceiling(
                              step$Q_i / product_info$Packsize[1] +
                                outstanding_requis_quantity),
                            days_to_order = step$Delta_i,
                            Drug_Name = product_info$Drug_name[1],
                            Supplier_Tradename = product_info$SupplierTradeName[1],
                            PackSize = product_info$Packsize[1],
                            Units = product_info$PrintForm[1],
                            Stock_Level = product_info$stocklvl[1],
                            Last_Order_Date = product_info$lastordered[1],
                            Min_stock_risk = risk_of_min_stock,
                            Min_stock_level = min_stock_level,
                            Comments = comment)

    # updateProgress(detail = paste0("Drug: ", to_return$drug, ", order:",
    #                                to_return$order))

    # Populate order quantity & time 'til next order in output table
    return(to_return)

  })

}

test_that("inventory_functions works", {

  load("pharmacy.rda")

  lower_lead = 1.5 # discrete min - 0.5
  upper_lead = 10.5 # discrete max - 0.5
  mode_lead = 4 # discrete mode

  # set distribution for delivery lead time

  lead_time_distribution <- distr6::Triangular$new(lower = lower_lead, 
                                                   upper = upper_lead,
                                                   mode = mode_lead)

  lead_time_dis <- lead_time_distribution$pdf(
    seq(lower_lead - 0.5, upper_lead + 0.5, 1))

  lead_time_dis <- lead_time_distribution$pdf(
    seq(lower_lead - 0.5, upper_lead + 0.5, 1))

  test_data <- pharmacy %>%
    dplyr::filter(Site1 == "Site C", NSVCode == "Drug A")

  daily_data <- make_tsibble(test_data, frequency = "Daily")

  test_forecast <- forecast_series(daily_data, length(lead_time_dis) + 14,
                                   frequency = "Daily")

  test_stock <- drug_quantity(forecast = test_forecast,
                              outstanding_orders = 100,
                              distribution = lead_time_dis,
                              min_stock = 100,
                              max_stock = 10000,
                              p_min = 0.01,
                              p_max = 0.05,
                              inv_i = 150,
                              delta_pref = 14)

  testthat::expect_equal(test_stock$Q_i, 10091.85, tolerance = 1)
  testthat::expect_equal(test_stock$Delta_i, 5)


})

# Run but keep eval=FALSE to avoid infinite loop
# Execute in the console directly
fusen::inflate(flat_file = "dev/flat_inventory_functions.Rmd", 
               vignette_name = NA)

CDU-data-science-team/pharmacyReporting documentation built on March 24, 2023, 2:32 p.m.