From f17c71e7bb766d979d7b3d4b6d3274e9de0eb0dc Mon Sep 17 00:00:00 2001
From: Sean Kent <41379385+skent259@users.noreply.github.com>
Date: Mon, 19 Feb 2024 10:23:26 -0600
Subject: [PATCH] Reduce `n_pairs` in `brsmatch()` when too many are specified
 (#18)

* Fix `brsmatch()` to reduce the number of pairs if too many for optimization
* Add testing to check in multiple scenarios
* Increment version number to 0.2.0.9002
* Update brsmatch documentation and NEWS.md
---
 DESCRIPTION                    |  2 +-
 NEWS.md                        |  1 +
 R/brsmatch.R                   | 86 +++++++++++++++++++++++++++-------
 man/brsmatch.Rd                |  6 +--
 tests/testthat/test-brsmatch.R | 46 ++++++++++++++++++
 5 files changed, 119 insertions(+), 22 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index f4263bc..c292265 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: rsmatch
 Title: Matching Methods for Time-Varying Observational Studies
-Version: 0.2.0.9001
+Version: 0.2.0.9002
 Authors@R: c(
     person("Sean", "Kent", , "skent259@gmail.com", role = c("aut", "cre", "cph"),
            comment = c(ORCID = "0000-0001-8697-9069")),
diff --git a/NEWS.md b/NEWS.md
index f635acc..f7b388b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,6 @@
 # rsmatch (development version)
 * Update `brsmatch()` and `coxpsmatch()` to handle NA rows via removing them
+* Fix `brsmatch()` to iteratively reduce the number of pairs if too many are specified
 
 # rsmatch 0.2.0
 
diff --git a/R/brsmatch.R b/R/brsmatch.R
index b209843..c558a76 100644
--- a/R/brsmatch.R
+++ b/R/brsmatch.R
@@ -11,9 +11,9 @@
 #'
 #' Note that when using exact matching, the `n_pairs` are split roughly in
 #' proportion to the number of treated subjects in each exact matching group.
-#' This has a possibility of failing  when `n_pairs` is large.  If this happens
-#' to you, we suggest manually performing exact matching, for example with
-#' `split()`, and selecting `n_pairs` for each group interactively.
+#' If you would like to control `n_pairs` exactly, we suggest manually
+#' performing exact matching, for example with `split()`, and selecting
+#' `n_pairs` for each group interactively.
 #'
 #' @param n_pairs The number of pairs desired from matching.
 #' @param data A data.frame or similar containing columns matching the `id,
@@ -201,21 +201,11 @@ brsmatch <- function(
   model <- .rsm_optimization_model(n_pairs, edges, bal, optimizer, verbose, balance)
 
   .print_if(verbose, "Preparing to run optimization model")
-  if (optimizer == "gurobi") {
-    res <- gurobi::gurobi(model, params = list(OutputFlag = 1 * verbose))
-    matches <- res$x[grepl("f", model$varnames)]
-  } else if (optimizer == "glpk") {
-    res <- Rglpk::Rglpk_solve_LP(
-      model$obj,
-      model$mat,
-      model$dir,
-      model$rhs,
-      types = model$types,
-      max = model$max,
-      control = list(verbose = verbose, presolve = TRUE)
-    )
-    matches <- res$solution[grepl("f", model$varnames)]
-  }
+  res <- .solve_or_reduce_pairs(n_pairs, model, optimizer, verbose)
+  matches <- switch(optimizer,
+    "gurobi" = res$x[grepl("f", model$varnames)],
+    "glpk" = res$solution[grepl("f", model$varnames)]
+  )
 
   matched_ids <- edges[matches == 1, c("trt_id", "all_id"), drop = FALSE]
   return(matched_ids)
@@ -592,3 +582,63 @@ brsmatch <- function(
   }
   return(model)
 }
+
+
+#' Solve brsmatch model even if too many pairs specified
+#'
+#' If the `n_pairs` is too large, the model will be infeasible, and will return
+#' a status code indicating this. This function will iteratively reduce the
+#' number of pairs until the model becomes solvable, then will return the
+#' solution with a warning.
+#'
+#' NOTE: gurobi functionality is untested, as I have since lost my license.
+#' Code is based on the documentation at
+#' https://www.gurobi.com/documentation/current/refman/r_grb.html
+#'
+#' @inheritParams brsmatch
+#' @param model The model output from `.rsm_optimization_model()`
+#'
+#' @return The result from [Rglpk::Rglpk_solve_LP] after possible n_pair
+#'   reduction.
+#'
+#' @noRd
+.solve_or_reduce_pairs <- function(n_pairs, model, optimizer, verbose) {
+  n_pairs_solve <- n_pairs
+
+  while (TRUE) {
+    if (optimizer == "gurobi") {
+      res <- gurobi::gurobi(model, params = list(OutputFlag = 1 * verbose))
+    } else if (optimizer == "glpk") {
+      res <- Rglpk::Rglpk_solve_LP(
+        model$obj,
+        model$mat,
+        model$dir,
+        model$rhs,
+        types = model$types,
+        max = model$max,
+        control = list(verbose = verbose, presolve = TRUE)
+      )
+    }
+    solved <- switch(optimizer,
+      "gurobi" = res$status == "OPTIMAL",
+      "glpk" = res$status == 0,
+    )
+    if (solved) {
+      break
+    }
+
+    n_pairs_solve <- n_pairs_solve - 1
+    # n pairs only appears in the first two model RHS constraints
+    model$rhs[1:2] <- c(n_pairs_solve, -n_pairs_solve)
+  }
+
+  if (n_pairs_solve != n_pairs) {
+    rlang::warn(
+      paste(
+        "Number of pairs reduced from", n_pairs, "to",
+        n_pairs_solve, "to create a solveable model."
+      )
+    )
+  }
+  return(res)
+}
diff --git a/man/brsmatch.Rd b/man/brsmatch.Rd
index 8e3670e..c1a7d7b 100644
--- a/man/brsmatch.Rd
+++ b/man/brsmatch.Rd
@@ -80,9 +80,9 @@ individual is matched to one other individual.
 \details{
 Note that when using exact matching, the \code{n_pairs} are split roughly in
 proportion to the number of treated subjects in each exact matching group.
-This has a possibility of failing  when \code{n_pairs} is large.  If this happens
-to you, we suggest manually performing exact matching, for example with
-\code{split()}, and selecting \code{n_pairs} for each group interactively.
+If you would like to control \code{n_pairs} exactly, we suggest manually
+performing exact matching, for example with \code{split()}, and selecting
+\code{n_pairs} for each group interactively.
 }
 \examples{
 if (requireNamespace("Rglpk", quietly = TRUE)) {
diff --git a/tests/testthat/test-brsmatch.R b/tests/testthat/test-brsmatch.R
index d00a674..5d1d4c7 100644
--- a/tests/testthat/test-brsmatch.R
+++ b/tests/testthat/test-brsmatch.R
@@ -442,3 +442,49 @@ test_that("brsmatch works when some input are NA", {
   # can't hold everyone's hand all the time...
 })
 
+
+test_that("brsmatch works when too many pairs specified (#11)", {
+  brsmatch_on_oasis <- function(n_pairs, ...) {
+    brsmatch(
+      n_pairs, oasis,
+      id = "subject_id", time = "visit", trt_time = "time_of_ad",
+      ...
+    )
+  }
+
+  check_for_glpk()
+
+  # With `balance = FALSE`
+  too_many <- brsmatch_on_oasis(n_pairs = 14, balance = FALSE) %>%
+    expect_warning("Number of pairs reduced")
+  enough <- brsmatch_on_oasis(n_pairs = 13, balance = FALSE)
+  expect_equal(too_many, enough)
+
+  # Way too many pairs
+  way_too_many <- brsmatch_on_oasis(n_pairs = 55, balance = FALSE) %>%
+    expect_warning("Number of pairs reduced")
+  expect_equal(way_too_many, too_many)
+
+  # With `balance = TRUE`
+  too_many <- brsmatch_on_oasis(n_pairs = 14, balance = TRUE) %>%
+    expect_warning("Number of pairs reduced")
+  enough <- brsmatch_on_oasis(n_pairs = 13, balance = TRUE)
+  expect_equal(too_many, enough)
+
+  # With `exact_match` variable
+  too_many <- brsmatch_on_oasis(n_pairs = 14, balance = FALSE, exact_match = c("m_f")) %>%
+    expect_warning("Number of pairs reduced")
+  enough <- brsmatch_on_oasis(n_pairs = 13, balance = FALSE, exact_match = c("m_f"))
+  expect_equal(too_many, enough)
+  # The warning message isn't ideal, but its okay for now
+
+
+  # With gurobi
+  check_for_gurobi()
+
+  too_many <- brsmatch_on_oasis(n_pairs = 14, balance = FALSE, options = list(optimizer = "gurobi")) %>%
+    expect_warning("Number of pairs reduced")
+  enough <- brsmatch_on_oasis(n_pairs = 13, balance = FALSE, options = list(optimizer = "gurobi"))
+  expect_equal(too_many, enough)
+
+})