Package 'couplr'

Description

Returns the character vector of method strings for which a trace function has been registered. Use any of these with lap_animate().

Usage

animated_methods()

Value

Character vector of registered method names.

Description

Turns a tidy assignment result back into a 0/1 assignment matrix.

Usage

as_assignment_matrix(x, n_sources = NULL, n_targets = NULL)

Arguments

Value

Integer matrix with 0 and 1 entries

Description

Constructs a matchit-class S3 object from a couplr result, enabling use with any function that accepts MatchIt objects (e.g., cobalt, marginaleffects).

Usage

as_matchit(
  result,
  left,
  right,
  formula = NULL,
  left_id = "id",
  right_id = "id",
  estimand = NULL,
  ...
)

Arguments

Value

An S3 object of class "matchit" with fields:

match.matrix

Match matrix (treated x controls)

treat

Named treatment vector (1/0)

weights

Matching weights

X

Covariate matrix

call

Original call

info

Metadata from couplr

Examples

## Not run: 
left <- data.frame(id = 1:5, age = c(25, 35, 45, 55, 65))
right <- data.frame(id = 6:15, age = runif(10, 20, 70))
result <- match_couples(left, right, vars = "age")
mi <- as_matchit(result, left, right)
# Now use with cobalt:
cobalt::bal.tab(mi)

## End(Not run)

Description

Solve the linear assignment problem (minimum- or maximum-cost matching) using several algorithms. Forbidden edges can be marked as NA or Inf.

Usage

assignment(
  cost,
  maximize = FALSE,
  method = c("auto", "jv", "hungarian", "munkres", "auction", "auction_gs",
    "auction_scaled", "sap", "ssp", "sap_dense", "csflow", "hk01", "bruteforce",
    "ssap_bucket", "cycle_cancel", "gabow_tarjan", "lapmod", "csa", "ramshaw_tarjan",
    "push_relabel", "network_simplex"),
  auction_eps = NULL,
  eps = NULL,
  memory_mode = "auto",
  certify = NULL,
  cardinality = c("complete", "maximum", "fixed"),
  n_matches = NULL,
  unmatched_penalty = NULL
)

Arguments

Details

method = "auto" selects an algorithm based on problem size and the costs:

• Very small (n <= 8 and m <= 8): "bruteforce" — exact enumeration

• Binary/constant costs: "hk01" — specialized for 0/1 costs

• Otherwise: "jv"

explain_dispatch() reports which rule fired and why. The other solvers are available by naming them explicitly.

Value

An object of class lap_solve_result, a list with elements:

• match — integer vector of length min(nrow(cost), ncol(cost)) giving the assigned column for each row (0 if unassigned).

• total_cost — numeric scalar, the objective value.

• status — character scalar drawn from solver_status_values(), computed from what the solver terminated on. "optimal" means the solver reached its own optimality condition with every row matched; it is not a checked proof. Use verify_assignment() for that.

• method_used — character scalar, the algorithm actually used.

• dispatch — list recording how method was chosen: the rule that fired under "auto", the condition that triggered it, and whether the method was named explicitly. See explain_dispatch().

• certificate — an assignment_certificate, present when one was checked. See certify.

• u, v — row and column duals in the orientation of cost, present when the solver produced optimal duals for cost: "jv" (and so "auto" on every problem it sends there), "hungarian", the lazy "jv" path, and memory_mode = "implicit", where they are the duals the last restricted master produced. verify_assignment() reads them instead of solving again. Absent under cardinality = "maximum" and "fixed", whose solver duals belong to the padded problem rather than to cost, and for "gabow_tarjan", whose duals are those of its scaled integer instance.

Under memory_mode = "implicit" the result also carries search: the columns the first round gave each row (seed_width), the pairs the candidate set ended up holding (candidate_edges) out of possible_edges, the pairs a cost was computed for (edges_evaluated), the round count, and rounds, one row per round of what the master held, what priced out and what each step cost.

Integer conversion for bit-scaling

"gabow_tarjan" is a bit-scaling algorithm and runs on integer costs. The conversion is part of the method, so the rule it follows and what it claims afterwards are both properties of the solve rather than of the caller's preparation.

Finite costs are shifted so the smallest is zero, multiplied by a scale factor s, and rounded to the nearest integer. s is set so that K * s * (max - min) stays at or below 10^13, where K is n + 1 on a square problem and 2 * min(n, m) + 1 on a rectangular one. K is the separation the algorithm needs between the optimum and a 1-optimal matching, and the bound holds every scaled cost clear of the sentinel that marks a forbidden pair while keeping the path sums the solver forms inside 64-bit arithmetic. Forbidden pairs take the sentinel and take part in neither the range nor the conversion.

A matrix whose finite entries are each within 1e-9 of an integer is taken as an integer matrix: s is one, the shift is an integer, and every cost reaches the solver as the integer nearest to it. Where those entries are integers the conversion is exact and the optimum is the optimum of the instance as supplied; where they are merely near one, each cost moves by at most 1e-9, so the matching returned costs at most 2 * min(n, m) * 1e-9 more than that optimum. Such a matrix is refused when its range exceeds 1.25 * 10^14 / K, since a scale of one is the only one available and no choice brings the instance inside the bound; the error names the limit and points at "jv" or "auction".

Costs that are not integers are solved on the rounded instance. Rounding moves each cost by at most 1 / (2 * s), so the matching returned costs at most min(n, m) * K * (max - min) / 10^13 more than the optimum of the matrix as supplied. On a square problem that is about n^2 * (max - min) / 10^13, which is 10^-7 of the range at n = 1000 and 10^-5 of it at n = 10000. Duals are divided by s and shifted back, so they belong to the original matrix; whether the matching is optimal for that matrix and not only for the rounded one is a question verify_assignment() answers rather than one this bound settles.

See Also

• lap_solve() — Tidy interface returning tibbles

• lap_solve_kbest() — Find k-best assignments ('Murty' algorithm)

• assignment_duals() — Extract dual variables for sensitivity analysis

• bottleneck_assignment() — Minimize maximum edge cost (minimax)

• sinkhorn() — Entropy-regularized optimal transport

Examples

cost <- matrix(c(4,2,5, 3,3,6, 7,5,4), nrow = 3, byrow = TRUE)
res  <- assignment(cost)
res$match; res$total_cost

Description

Solves the linear assignment problem and returns dual potentials (u, v) in addition to the optimal matching. The dual variables provide an optimality certificate and enable sensitivity analysis.

Usage

assignment_duals(cost, maximize = FALSE, certify = FALSE)

Arguments

Details

The dual variables satisfy the complementary slackness conditions:

• For minimization: u[i] + v[j] <= cost[i,j] for all (i,j)

• For any assigned pair (i,j): u[i] + v[j] = cost[i,j]

This implies that sum(u) + sum(v) = total_cost (strong duality).

Applications of dual variables:

• Optimality verification: Check that duals satisfy constraints

• Sensitivity analysis: Reduced cost c[i,j] - u[i] - v[j] shows how much an edge cost must decrease before it enters the solution

• Pricing in column generation: Use duals to price new columns

• Warm starting: Reuse duals when costs change slightly

Value

A list with class "assignment_duals_result" containing:

• match - integer vector of column assignments (1-based)

• total_cost - optimal objective value

• u - numeric vector of row dual variables (length n)

• v - numeric vector of column dual variables (length m)

• status - character, e.g. "optimal"

• certificate - an assignment_certificate, present only under certify = TRUE

See Also

assignment() for standard assignment without duals, verify_assignment() for the check certify = TRUE runs

Examples

cost <- matrix(c(4, 2, 5, 3, 3, 6, 7, 5, 4), nrow = 3, byrow = TRUE)
result <- assignment_duals(cost)

# Check optimality: u + v should equal cost for assigned pairs
for (i in 1:3) {
  j <- result$match[i]
  cat(sprintf("Row %d -> Col %d: u + v = %.2f, cost = %.2f\n",
              i, j, result$u[i] + result$v[j], cost[i, j]))
}

# Verify strong duality
cat("sum(u) + sum(v) =", sum(result$u) + sum(result$v), "\n")
cat("total_cost =", result$total_cost, "\n")

# Reduced costs (how much must cost decrease to enter solution)
reduced <- outer(result$u, result$v, "+")
reduced_cost <- cost - reduced
print(round(reduced_cost, 2))

Description

S3 method for augmenting matching results following the broom package conventions. This is a thin wrapper around join_matched() with sensible defaults for quick exploration.

Usage

## S3 method for class 'matching_result'
augment(x, left, right, ...)

Arguments

Details

This method follows the augment() convention from the broom package, making it easy to integrate couplr into tidymodels workflows. It's equivalent to calling join_matched() with default parameters.

If the broom package is not loaded, you can use couplr::augment() to access this function.

Value

A tibble with matched pairs and original data (see join_matched())

Examples

left <- data.frame(
  id = 1:5,
  treatment = 1,
  age = c(25, 30, 35, 40, 45)
)

right <- data.frame(
  id = 6:10,
  treatment = 0,
  age = c(24, 29, 36, 41, 44)
)

result <- match_couples(left, right, vars = "age")
augment(result, left, right)

Description

Produces ggplot2-based balance assessment plots. Returns a ggplot object.

Usage

## S3 method for class 'balance_diagnostics'
autoplot(
  object,
  type = c("love", "histogram", "variance"),
  threshold = 0.1,
  ...
)

Arguments

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(42)
  left <- data.frame(id = 1:10, age = rnorm(10, 45, 10),
                     income = rnorm(10, 50000, 15000))
  right <- data.frame(id = 11:30, age = rnorm(20, 47, 10),
                       income = rnorm(20, 52000, 15000))
  result <- match_couples(left, right, vars = c("age", "income"))
  bal <- balance_diagnostics(result, left, right, vars = c("age", "income"))
  ggplot2::autoplot(bal)
}

Description

Produces ggplot2-based visualizations of matching distance distributions. Returns a ggplot object that can be further customized.

Usage

## S3 method for class 'matching_result'
autoplot(object, type = c("histogram", "density", "ecdf"), ...)

Arguments

Details

Use plot() for base graphics or autoplot() for ggplot2 output. The ggplot2 package must be installed.

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  left <- data.frame(id = 1:5, x = c(1, 2, 3, 4, 5))
  right <- data.frame(id = 6:10, x = c(1.1, 2.2, 3.1, 4.2, 5.1))
  result <- match_couples(left, right, vars = "x")
  ggplot2::autoplot(result)
  ggplot2::autoplot(result, type = "density")
}

Description

Plots p-value upper bounds against sensitivity parameter Gamma.

Usage

## S3 method for class 'sensitivity_analysis'
autoplot(object, alpha = 0.05, ...)

Arguments

Value

A ggplot object

Description

S3 method enabling cobalt::bal.tab() on couplr result objects. Requires the cobalt package to be installed.

Usage

## S3 method for class 'matching_result'
bal.tab(x, left, right, ...)

## S3 method for class 'full_matching_result'
bal.tab(x, left, right, ...)

## S3 method for class 'cem_result'
bal.tab(x, left, right, ...)

## S3 method for class 'subclass_result'
bal.tab(x, data = NULL, ...)

Arguments

Details

These methods convert couplr results to the format cobalt expects (a matchit-class object) and then delegate to cobalt's own bal.tab.matchit() method. The cobalt package must be installed but is not required for couplr to function.

Value

A cobalt balance table object

Description

Computes comprehensive balance statistics comparing the distribution of matching variables between left and right units in the matched sample.

Usage

balance_diagnostics(result, ...)

## S3 method for class 'matching_result'
balance_diagnostics(
  result,
  left,
  right,
  vars = NULL,
  left_id = "id",
  right_id = "id",
  ...
)

## S3 method for class 'full_matching_result'
balance_diagnostics(
  result,
  left,
  right,
  vars = NULL,
  left_id = "id",
  right_id = "id",
  ...
)

## S3 method for class 'cem_result'
balance_diagnostics(
  result,
  left,
  right,
  vars = NULL,
  left_id = "id",
  right_id = "id",
  ...
)

## S3 method for class 'subclass_result'
balance_diagnostics(result, data = NULL, vars = NULL, ...)

Arguments

Details

This function computes several balance metrics:

Standardized Difference: The difference in means divided by the pooled standard deviation. Values less than 0.1 indicate excellent balance, 0.1-0.25 good balance.

Variance Ratio: The ratio of variances (left/right). Values close to 1 are ideal; the conventional acceptable range is 0.5 to 2. For weighted estimators (full matching, CEM, subclassification) the means and variances are computed using the stratum weights.

KS Statistic: Kolmogorov-Smirnov test statistic comparing distributions. Lower values indicate more similar distributions.

Overall Metrics include mean absolute standardized difference across all variables, proportion of variables with large imbalance (|std diff| > 0.25), and maximum standardized difference.

Value

An S3 object of class balance_diagnostics containing:

var_stats

Tibble with per-variable balance statistics

overall

List with overall balance metrics

pairs

Tibble of matched pairs with variables

n_matched

Number of matched pairs

n_unmatched_left

Number of unmatched left units

n_unmatched_right

Number of unmatched right units

method

Matching method used

has_blocks

Whether blocking was used

block_stats

Per-block statistics (if blocking used)

Examples

# Create sample data
set.seed(123)
left <- data.frame(
  id = 1:10,
  age = rnorm(10, 45, 10),
  income = rnorm(10, 50000, 15000)
)
right <- data.frame(
  id = 11:30,
  age = rnorm(20, 47, 10),
  income = rnorm(20, 52000, 15000)
)

# Match
result <- match_couples(left, right, vars = c("age", "income"))

# Get balance diagnostics
balance <- balance_diagnostics(result, left, right, vars = c("age", "income"))
print(balance)

# Get balance table
balance_table(balance)

Description

Formats balance diagnostics into a clean table for display or export.

Usage

balance_table(balance, digits = 3)

Arguments

Value

A tibble with formatted balance statistics

Description

Finds an assignment that minimizes (or maximizes) the maximum edge cost in a perfect matching. Unlike standard LAP which minimizes the sum of costs, BAP minimizes the maximum (bottleneck) cost.

Usage

bottleneck_assignment(cost, maximize = FALSE)

Arguments

Details

The Bottleneck Assignment Problem (BAP) is a variant of the Linear Assignment Problem where instead of minimizing the sum of assignment costs, we minimize the maximum cost among all assignments (minimax objective).

Algorithm: Uses binary search on the sorted unique costs combined with 'Hopcroft-Karp' bipartite matching to find the minimum threshold that allows a perfect matching.

Complexity: O(E * sqrt(V) * log(unique costs)) where E = edges, V = vertices.

Applications:

• Task scheduling with deadline constraints (minimize latest completion)

• Resource allocation (minimize maximum load/distance)

• Network routing (minimize maximum link utilization)

• Fair division problems (minimize maximum disparity)

Value

A list with class "bottleneck_result" containing:

• match - integer vector of length nrow(cost) giving the assigned column for each row (1-based indexing)

• bottleneck - numeric scalar, the bottleneck (max/min edge) value

• status - character scalar, e.g. "optimal"

See Also

assignment() for standard LAP (sum objective), lap_solve() for tidy LAP interface

Examples

# Simple example: minimize max cost
cost <- matrix(c(1, 5, 3,
                 2, 4, 6,
                 7, 1, 2), nrow = 3, byrow = TRUE)
result <- bottleneck_assignment(cost)
result$bottleneck  # Maximum edge cost in optimal assignment

# Maximize minimum (fair allocation)
profits <- matrix(c(10, 5, 8,
                    6, 12, 4,
                    3, 7, 11), nrow = 3, byrow = TRUE)
result <- bottleneck_assignment(profits, maximize = TRUE)
result$bottleneck  # Minimum profit among all assignments

# With forbidden assignments
cost <- matrix(c(1, NA, 3,
                 2, 4, Inf,
                 5, 1, 2), nrow = 3, byrow = TRUE)
result <- bottleneck_assignment(cost)

Description

Maximizes the number of matched pairs subject to balance constraints, and reports how far the matched sample sits from the largest one the constraints admit.

Usage

cardinality_match(
  left,
  right,
  vars,
  fine = NULL,
  refined = NULL,
  refined_exact = 1L,
  moments = NULL,
  max_std_diff = Inf,
  engine = c("auto", "flow", "branch_bound", "heuristic"),
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  sigma = NULL,
  max_distance = Inf,
  calipers = NULL,
  left_id = NULL,
  right_id = NULL,
  time_limit = 30,
  node_limit = 500L,
  method = "auto",
  max_iter = 100L,
  batch_fraction = 0.1,
  memory_mode = "auto"
)

Arguments

Details

Cardinality matching (Zubizarreta 2012) asks for the largest matched sample that satisfies stated balance constraints, rather than the closest one. The constraints come in two kinds, and which kinds are stated decides how the problem is solved.

Fine and refined covariate balance are representable in the matching network itself (Rosenbaum, Ross and Silber 2007; Pimentel, Kelz, Silber and Rosenbaum 2015): every category becomes a node, and a single min-cost flow solve returns the largest balanced sample together with a dual certificate of its optimality, at polynomial cost. max_std_diff defaults to Inf, stating no moment constraint, so a call asking for fine or refined balance alone takes this path and comes back certified:

cardinality_match(left, right, vars, fine = "region")

Linear moment constraints – a bound on a standardized or mean difference – cut across the network and are dualized instead. Their multipliers are updated per node, the Lagrangian value bounds the subtree, and the search branches on unit inclusion until the bound meets the incumbent or a budget runs out. A finite max_std_diff states one such constraint per variable and per direction, so it is what moves a call onto this path, and node_limit and time_limit are what end it when the bound does not close. The report says which happened: certified is TRUE only when the search settled and the gap is zero, and gap says how many matched units separate the answer from the bound when it is not. Every stopping path returns a matched set that satisfies every stated constraint, together with a bound valid for the whole problem.

The certificate carries no tolerance. Each node's bound is the dual objective of its Lagrangian relaxation evaluated in exact arithmetic, which weak duality makes a lower bound whatever the solver returned, and the pruning tests, the constraint values of a matched set and the cardinality read off the bound are decided exactly, over the distances and constraint coefficients as stored.

How long the search runs depends on whether the moment bounds bind. When the distance-minimizing match already satisfies them, which happens with a loose bound or with a right pool large enough that the closest partners are already balanced, the root node certifies and the call returns in milliseconds. When they bind, every node costs a flow solve, and although the root enters a feasible incumbent before the first node is opened, the bound above it closes slowly, so a call can spend its whole budget and still report a gap. Set node_limit and time_limit to a budget you are willing to spend, and read stopped_on and gap to see what the budget bought.

The objective is lexicographic. Matched cardinality comes first, total distance second, so a pair is never added at the cost of a match and never dropped to shorten one.

engine = "heuristic" runs a pruning loop instead: a full match, then repeated deletion of the pairs carrying the worst variable's imbalance. Its properties follow from the loop. It starts from a full match and deletes pairs, and never re-adds one, so it cannot recover from an over-aggressive deletion. It never maximizes cardinality, it only shrinks a starting match. It considers one variable per iteration, the current worst, so balance on the others is incidental. It computes no bound, so it cannot say how far its answer sits from the largest balanced sample; best_possible and gap are NA and certified is FALSE.

Value

A matching_result object. Beyond the fields every matching carries it holds:

• cardinality - the report: n_matched, best_possible, gap, gap_fraction, certified, stopped_on, n_nodes, and the state of every stated constraint. See print.cardinality_report().

• certificate - present when the search certified optimality, with arithmetic = "exact": the bounds, the comparisons between them and the matched sample's objective and constraint values it rests on are all decided exactly.

• potentials - the dual potentials of the network solve the matched sample came from, a list with elements left and right holding one value per unit, named by id, in distance terms: the reduced cost of a pair is its distance less the two, with the tier weights and any moment multipliers already folded in. Under engine = "flow" that solve is the whole problem, and the potentials certify the matched sample as its optimum. Under branch and bound it is the node the incumbent was read from, with the arcs that node fixed and the moment rows priced at its multipliers, and they certify the incumbent as that node's optimum; the bound on the whole problem is cardinality$best_possible. Present when that solve was certified.

• status - "optimal", "iteration_limit", or "heuristic".

• info$engine - the engine that answered.

info$pruning_iterations and info$pairs_removed are carried by engine = "heuristic".

References

Zubizarreta, J. R. (2012). Using mixed integer programming for matching in an observational study of kidney failure after surgery. Journal of the American Statistical Association, 107(500), 1360-1371.

Rosenbaum, P. R., Ross, R. N. and Silber, J. H. (2007). Minimum distance matched sampling with fine balance in an observational study of treatment for ovarian cancer. Journal of the American Statistical Association, 102(477), 75-83.

Pimentel, S. D., Kelz, R. R., Silber, J. H. and Rosenbaum, P. R. (2015). Large, sparse optimal matching with refined covariate balance in an observational study of the health outcomes produced by new surgeons. Journal of the American Statistical Association, 110(510), 515-527.

See Also

match_couples() for distance-minimizing matching, balance_diagnostics() for reading balance off any matched sample.

Examples

set.seed(42)
left <- data.frame(id = 1:20, x = rnorm(20), y = rnorm(20),
                   region = rep(c("A", "B"), length.out = 20))
right <- data.frame(id = 21:50, x = rnorm(30, 0.5), y = rnorm(30, 0.3),
                    region = rep(c("A", "B"), length.out = 30))

# Exact fine balance on region, no moment constraint: one flow solve,
# answered with a certificate.
fit <- cardinality_match(left, right, vars = c("x", "y"),
                         fine = "region")
fit$cardinality

# A standardized-difference bound as well: the same match, searched by
# branch and bound under a small node budget.
bb <- cardinality_match(left, right, vars = c("x", "y"),
                        fine = "region", max_std_diff = 0.1,
                        node_limit = 25L)
bb$cardinality

Description

Coarsens continuous variables into bins, then performs exact matching on the coarsened values. Units in strata containing both left and right units are kept; others are pruned. Matched units receive weights inversely proportional to stratum sizes to maintain balance.

Usage

cem_match(
  left,
  right,
  vars,
  cutpoints = NULL,
  n_bins = "sturges",
  grouping = NULL,
  keep = "all",
  left_id = "id",
  right_id = "id"
)

Arguments

Details

CEM algorithm:

1. Coarsen each numeric variable using cut with either user-specified breakpoints or automatic binning (Sturges, FD, or Scott rule)

2. Categorical variables in grouping are kept as-is

3. Create strata by concatenating all coarsened values

4. Drop strata with 0 left or 0 right units

5. Compute CEM weights: left units get weight 1, right units get weight n_left_in_stratum / n_right_in_stratum so that the total weight of right units in each stratum equals the number of left units

Value

An S3 object of class c("cem_result", "couplr_result") containing:

matched

Tibble with columns id, side, stratum, weight

strata_summary

Tibble with per-stratum counts

info

List with n_strata, n_matched_left, n_matched_right, n_pruned_left, n_pruned_right, method, vars

Examples

set.seed(42)
left <- data.frame(
  id = 1:20, age = rnorm(20, 40, 10),
  income = rnorm(20, 50000, 10000)
)
right <- data.frame(
  id = 21:60, age = rnorm(40, 42, 10),
  income = rnorm(40, 52000, 10000)
)
result <- cem_match(left, right, vars = c("age", "income"))
print(result)

Description

Precomputes a distance matrix between left and right datasets, allowing it to be reused across multiple matching operations with different constraints. This is particularly useful when exploring different matching parameters (max_distance, calipers, methods) without recomputing distances.

Usage

compute_distances(
  left,
  right,
  vars,
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  left_id = "id",
  right_id = "id",
  block_id = NULL,
  memory_mode = "auto"
)

Arguments

Details

This function computes distances once and stores them in a reusable object. The resulting distance_object can be passed to match_couples() (optionally with method = "greedy") instead of providing datasets and variables.

Benefits:

• Performance: Avoid recomputing distances when trying different constraints

• Exploration: Quickly test max_distance, calipers, or methods

• Consistency: Ensures same distances used across comparisons

With the default memory_mode = "auto"/"dense", the stored cost_matrix is a plain dense matrix; couplr has no sparse matrix representation, so forbidden pairs (Inf) still occupy a cell rather than being dropped from storage. With memory_mode = "lazy", cost_matrix is instead a lazy_cost_spec that recomputes distances on demand from the underlying feature data.

The distance_object stores the original datasets, allowing downstream functions like join_matched() to work seamlessly.

Value

An S3 object of class "distance_object" containing:

• cost_matrix: Numeric matrix of distances

• left_ids: Character vector of left IDs

• right_ids: Character vector of right IDs

• block_id: Block ID column name (if specified)

• metadata: List with computation details (vars, distance, scale, etc.)

• original_left: Original left dataset (for later joining)

• original_right: Original right dataset (for later joining)

Examples

# Compute distances once
left <- data.frame(id = 1:5, age = c(25, 30, 35, 40, 45), income = c(45, 52, 48, 61, 55) * 1000)
right <- data.frame(id = 6:10, age = c(24, 29, 36, 41, 44), income = c(46, 51, 47, 60, 54) * 1000)

dist_obj <- compute_distances(
  left, right,
  vars = c("age", "income"),
  scale = "standardize"
)

# Reuse for different matching strategies
result1 <- match_couples(dist_obj, max_distance = 0.5)
result2 <- match_couples(dist_obj, max_distance = 1.0)
result3 <- match_couples(dist_obj, method = "greedy", strategy = "sorted")

# All use the same precomputed distances

Description

Comprehensive diagnostics for a distance matrix with actionable suggestions.

Usage

diagnose_distance_matrix(
  cost_matrix,
  left = NULL,
  right = NULL,
  vars = NULL,
  warn = TRUE
)

Arguments

Value

List with diagnostic results and suggestions

Description

What the matrix itself costs a process, which is more than the 8 bytes a cell occupies: matrix(0, n, m) at the R level is copied into a lap::CostMatrix (8B data + 4B mask), and the two coexist while garbage collection lags. Building the matrix and stopping there has peaked at between 1.5 and 1.9 times the raw cell bytes from 5,000 to 20,000 units on the memory benchmark, so the default multiplier is above what has been measured rather than fitted to it. n/m are coerced to double before multiplying so the estimate itself can't overflow the way lap::CostMatrix's old int flat-index arithmetic did.

Usage

estimate_dense_matrix_mb(n, m, overhead_factor = 4)

Arguments

Details

This is the matrix, not the solve. estimate_dense_solve_mb() is what the memory guard reads; see there for why the two differ by more than a copy.

Value

Numeric scalar, the estimated footprint of the matrix in megabytes.

Examples

estimate_dense_matrix_mb(5000, 5000)

Description

The matrix is the allocation a caller can reason about and it is not what bounds a dense solve. Preparation copies, the solver's own workspace and the assignment structures it carries all sit on top of the matrix, and the peak is what decides whether the solve fits, so that is the quantity the guard compares against available RAM.

Usage

estimate_dense_solve_mb(n, m, solve_factor = 12)

Arguments

Details

The multiplier is read off the measurement rather than from an enumeration of the copies, which has proved to understate it. On the memory benchmark – one fresh R session per arm, peak resident set taken from outside and read against an idle session that loaded the same packages – a dense one-to-one solve has peaked at between 7 and 11 times the raw matrix bytes from 5,000 to 20,000 units, with the figure at a given size moving between runs. The default of 12 sits above every measured peak, which is the direction a guard should err in: it exists to refuse a solve that will not fit, not to predict where the peak lands.

Value

Numeric scalar, the estimated peak footprint of a dense solve in megabytes.

Examples

estimate_dense_solve_mb(5000, 5000)

Description

Small example datasets for demonstrating couplr functionality across different assignment problem types: square, rectangular, sparse, and binary.

Usage

example_costs

Format

A list containing four example cost matrices:

simple_3x3

A 3x3 cost matrix with costs ranging from 2-7. Optimal assignment: row 1 -> col 2 (cost 2), row 2 -> col 1 (cost 3), row 3 -> col 3 (cost 4). Total optimal cost: 9.

rectangular_3x5

A 3x5 rectangular cost matrix demonstrating assignment when rows < columns. Each of 3 rows is assigned to one of 5 columns; 2 columns remain unassigned. Costs range 1-6.

sparse_with_na

A 3x3 matrix with NA values indicating forbidden assignments. Use this to test algorithms' handling of constraints. Position (1,3), (2,2), and (3,1) are forbidden.

binary_costs

A 3x3 matrix with binary (0/1) costs, suitable for testing the HK01 algorithm. Diagonal entries are 0 (preferred), off-diagonal entries are 1 (penalty).

Details

These matrices are designed to test different aspects of LAP solvers:

simple_3x3: Basic functionality test. Any correct solver should find total cost = 9.

rectangular_3x5: Tests handling of non-square problems. The optimal solution assigns all 3 rows with minimum total cost.

sparse_with_na: Tests constraint handling. Algorithms must avoid NA positions while finding an optimal assignment among valid entries.

binary_costs: Tests specialized binary cost algorithms. The optimal assignment uses all diagonal entries (total cost = 0).

See Also

lap_solve, example_df

Examples

# Simple 3x3 assignment
result <- lap_solve(example_costs$simple_3x3)
print(result)
# Optimal: sources 1,2,3 -> targets 2,1,3 with cost 9

# Rectangular problem (3 sources, 5 targets)
result <- lap_solve(example_costs$rectangular_3x5)
print(result)
# All 3 sources assigned; 2 targets unassigned

# Sparse problem with forbidden assignments
result <- lap_solve(example_costs$sparse_with_na)
print(result)
# Avoids NA positions

# Binary costs - test HK01 algorithm
result <- lap_solve(example_costs$binary_costs, method = "hk01")
print(result)
# Finds diagonal assignment (cost = 0)

Description

A tidy data frame representation of assignment problems, suitable for use with grouped workflows and batch solving. Contains two independent 3x3 assignment problems in long format.

Usage

example_df

Format

A tibble with 18 rows and 4 columns:

sim

Simulation/problem identifier. Integer with values 1 or 2, distinguishing two independent assignment problems. Use with group_by(sim) for grouped solving.

source

Source node index. Integer 1-3 representing the row (source) in each 3x3 cost matrix.

target

Target node index. Integer 1-3 representing the column (target) in each 3x3 cost matrix.

cost

Cost of assigning source to target. Numeric values ranging from 1-7. Each source-target pair has exactly one cost entry.

Details

This dataset demonstrates couplr's data frame interface for LAP solving. The long format (one row per source-target pair) is converted internally to a cost matrix for solving.

Simulation 1: Costs from example_costs$simple_3x3

• Optimal assignment: (1->2, 2->1, 3->3)

• Total cost: 9

Simulation 2: Different cost structure

• Optimal assignment: (1->1, 2->3, 3->3) or equivalent

• Total cost: 4

See Also

lap_solve, lap_solve_batch, example_costs

Examples

library(dplyr)

# Solve both problems with grouped workflow
example_df |>
  group_by(sim) |>
  lap_solve(source, target, cost)

# Batch solving for efficiency
example_df |>
  group_by(sim) |>
  lap_solve_batch(source, target, cost)

# Inspect the data structure
example_df |>
  group_by(sim) |>
  summarise(
    n_pairs = n(),
    min_cost = min(cost),
    max_cost = max(cost)
  )

Description

assignment(method = "auto") picks a solver from a small set of rules evaluated against one pass over the cost matrix. This reports the decision without solving: which rule fired, the property that triggered it, the rules that were tested first and did not fire, and the internal representation the problem will be handed to the solver in.

Usage

explain_dispatch(
  cost,
  method = "auto",
  maximize = FALSE,
  cardinality = c("complete", "maximum", "fixed"),
  n_matches = NULL,
  unmatched_penalty = NULL
)

## S3 method for class 'dispatch_explanation'
print(x, ...)

Arguments

Details

The same rule table drives the dispatch itself, so the report cannot drift from the behaviour.

Value

An object of class dispatch_explanation, a list with elements:

• method - the solver that will run.

• explicit - TRUE when method was named rather than dispatched.

• rule, condition, reason - the rule that fired under "auto".

• considered - data frame of every rule tested, in order, with whether it fired.

• n, m, transposed - problem shape, and whether the solver sees the transpose. Rectangular problems are transposed internally so the solver always has at least as many columns as rows.

• probe - the single-pass summary the rules read.

Invisibly returns x.

See Also

assignment()

Examples

explain_dispatch(matrix(runif(400), 20, 20))
explain_dispatch(matrix(sample(0:1, 400, TRUE), 20, 20))

Description

Assigns every unit (left and right) to a matched group with variable ratios (1:k or k:1). Unlike 1:1 matching, full matching does not discard units, producing matched groups where each group contains at least one left and one right unit.

Usage

full_match(
  left,
  right,
  vars,
  distance = "euclidean",
  min_controls = 1,
  max_controls = Inf,
  caliper = NULL,
  caliper_sd = NULL,
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  sigma = NULL,
  left_id = "id",
  right_id = "id",
  method = "optimal",
  memory_mode = "auto"
)

Arguments

Details

full_match() builds matched groups of variable size. Under the default min_controls = 1 it solves full matching in the sense of Hansen and Klopfer (2006): a group holds either one left unit and several right ones or one right unit and several left ones, and both shapes may appear in the same solution. max_controls bounds the many side there, whichever side that is.

A lower bound above one admits only the one-to-many shape, because a group built around a single right unit holds exactly one of them and cannot meet a lower bound of two. Every group is then one left unit with between min_controls and max_controls right ones. Both bounds count right units whichever side holds more, so an instance with too few right units to give every left unit that many is refused as "infeasible" rather than answered with groups counted the other way round.

Two algorithms are available:

Optimal (method = "optimal", default): Solves a min-cost max-flow problem that minimizes total distance across all group assignments simultaneously, with the optimum found via Dijkstra's algorithm with Johnson potentials. At min_controls = 1 the network is an edge cover over the admissible pairs, which is what lets a group be centred on either side; above one it is one centre per group and the centres are the left units.

Greedy (method = "greedy"): A fast two-pass heuristic:

1. Each left unit picks its nearest eligible right unit

2. Remaining right units are assigned to their nearest already-matched left unit, respecting max_controls

This is faster but does not guarantee globally optimal results.

Every left unit carries weight 1, and the right units of a group share a total weight equal to the number of left units in that group, so the two sides of a group weigh the same. A group holding one left unit and k right units therefore gives each right unit 1/k; a group holding k left units and one right unit gives that right unit k. "greedy" builds one left unit per group and produces only the first shape. "optimal" at min_controls = 1 centres a group on whichever side is larger, so both arise; above one the centres are the left units and the first shape is the only one.

Value

An S3 object of class c("full_matching_result", "couplr_result") containing:

groups

Tibble with columns group_id, id, side ("left"/"right"), and weight

status

What the solver terminated on, one of "optimal" (every unit placed in a group meeting min_controls, at minimum total distance), "partial" (groups formed, some units left over, under a maximum-cardinality-then-minimum-cost objective), "infeasible" (no group meets the requested bounds), or "heuristic" (method = "greedy", which neither claims nor checks optimality). See solver_status_values.

info

List with n_groups, n_left, n_right, n_unmatched_left, n_unmatched_right, method, vars

unmatched

Left and right IDs that no group contains. Every unit is either a row of groups or an entry here.

potentials

Node potentials from the flow solve, a list with elements left and right holding one value per unit. They are one representative of the dual, in the gauge the solver fixes. Present for method = "optimal" only.

certificate

A flow_certificate from verify_flow, checking the solved flow and its potentials against the optimality conditions. status says what the solver terminated on; this says what was checked. Present for method = "optimal" only. Under memory_mode = "implicit" the check covers the pairs the flow was solved over, and omitted_proven_floor bounds the reduced cost of every pair it omitted; certified_optimal holds both together.

search

Under memory_mode = "implicit": the seed width, the pairs generated (candidate_edges) against the pairs there are (possible_edges), the distances computed (edges_evaluated) and one row per round.

Examples

set.seed(42)
left <- data.frame(id = 1:5, age = c(25, 35, 45, 55, 65))
right <- data.frame(id = 6:20, age = runif(15, 20, 70))
result <- full_match(left, right, vars = "age")
print(result)

Description

Extract method used from assignment result

Usage

get_method_used(x)

Arguments

Value

Character string indicating method used

Description

Extract total cost from assignment result

Usage

get_total_cost(x)

Arguments

Value

Numeric total cost

Description

A comprehensive example dataset for demonstrating couplr functionality across vignettes. Contains synthetic hospital staff scheduling data with nurses, shifts, costs, and preference scores suitable for assignment problems, as well as nurse characteristics for matching workflows.

Usage

hospital_staff

Format

A list containing eight related datasets:

basic_costs

A 10x10 numeric cost matrix for assigning 10 nurses to 10 shifts. Values range from approximately 1-15, where lower values indicate better fit (less overtime, matches skills, respects preferences). Use with lap_solve() for basic assignment.

preferences

A 10x10 numeric preference matrix on a 0-10 scale, where higher values indicate stronger nurse preference for a shift. Use with lap_solve(..., maximize = TRUE) to optimize preferences rather than minimize costs.

schedule_df

A tibble with 100 rows (10 nurses x 10 shifts) in long format for data frame workflows:

nurse_id

Integer 1-10. Unique identifier for each nurse.

shift_id

Integer 1-10. Unique identifier for each shift.

cost

Numeric. Assignment cost (same values as basic_costs).

preference

Numeric 0-10. Nurse preference score.

skill_match

Integer 0/1. Binary indicator: 1 if nurse skills match shift requirements, 0 otherwise.

nurses

A tibble with 10 rows describing nurse characteristics:

nurse_id

Integer 1-10. Links to schedule_df and basic_costs rows.

experience_years

Numeric 1-20. Years of nursing experience.

department

Character. Primary department: "ICU", "ER", "General", or "Pediatrics".

shift_preference

Character. Preferred shift type: "day", "evening", or "night".

certification_level

Integer 1-3. Certification level where 3 is highest (e.g., 1=RN, 2=BSN, 3=MSN).

shifts

A tibble with 10 rows describing shift requirements:

shift_id

Integer 1-10. Links to schedule_df and basic_costs cols.

department

Character. Department needing coverage.

shift_type

Character. Shift type: "day", "evening", or "night".

min_experience

Numeric. Minimum years of experience required.

min_certification

Integer 1-3. Minimum certification level.

weekly_df

A tibble for batch solving with 500 rows (5 days x 10 nurses x 10 shifts):

day

Character. Day of week: "Mon", "Tue", "Wed", "Thu", "Fri".

nurse_id

Integer 1-10. Nurse identifier.

shift_id

Integer 1-10. Shift identifier.

cost

Numeric. Daily assignment cost (varies by day).

preference

Numeric 0-10. Daily preference score.

Use with group_by(day) for solving each day's schedule.

nurses_extended

A tibble with 200 nurses for matching examples, representing a treatment group (e.g., full-time nurses):

nurse_id

Integer 1-200. Unique identifier.

age

Numeric 22-65. Nurse age in years.

experience_years

Numeric 0-40. Years of nursing experience.

hourly_rate

Numeric 25-75. Hourly wage in dollars.

department

Character. Primary department assignment.

certification_level

Integer 1-3. Certification level.

is_fulltime

Logical. TRUE for full-time status.

controls_extended

A tibble with 300 potential control nurses (e.g., part-time or registry nurses) for matching. Same structure as nurses_extended. Designed to have systematic differences from nurses_extended (older, less experience on average) to demonstrate matching's ability to create comparable groups.

Details

This dataset is used throughout the couplr documentation to provide a consistent, realistic example that evolves in complexity. It supports three use cases: (1) basic LAP solving with cost matrices, (2) batch solving across multiple days, and (3) matching workflows comparing nurse groups.

The dataset is designed to demonstrate progressively complex scenarios:

Basic LAP (vignette("getting-started")):

• basic_costs: Simple 10x10 assignment

• preferences: Maximization problem

• schedule_df: Data frame input, grouped workflows

• weekly_df: Batch solving across days

Algorithm comparison (vignette("algorithms")):

• Use basic_costs to compare algorithm behavior

• Modify with NA values for sparse scenarios

Matching workflows (vignette("matching-workflows")):

• nurses_extended: Treatment group (full-time nurses)

• controls_extended: Control pool (part-time/registry nurses)

• Match on age, experience, department for causal analysis

See Also

lap_solve for basic assignment solving, lap_solve_batch for batch solving, match_couples for matching workflows, vignette("getting-started") for introductory tutorial

Examples

# Basic assignment: assign nurses to shifts minimizing cost
lap_solve(hospital_staff$basic_costs)

# Maximize preferences instead
lap_solve(hospital_staff$preferences, maximize = TRUE)

# Data frame workflow
library(dplyr)
hospital_staff$schedule_df |>
  lap_solve(nurse_id, shift_id, cost)

# Batch solve weekly schedule
hospital_staff$weekly_df |>
  group_by(day) |>
  lap_solve(nurse_id, shift_id, cost)

# Matching workflow: match full-time to part-time nurses
match_couples(
  left = hospital_staff$nurses_extended,
  right = hospital_staff$controls_extended,
  vars = c("age", "experience_years", "certification_level"),
  auto_scale = TRUE
)

Description

Check if Object is a Distance Object

Usage

is_distance_object(x)

Arguments

Value

Logical: TRUE if x is a distance_object

Examples

left <- data.frame(id = 1:3, x = c(1, 2, 3))
right <- data.frame(id = 4:6, x = c(1.1, 2.1, 3.1))
dist_obj <- compute_distances(left, right, vars = "x")
is_distance_object(dist_obj)  # TRUE
is_distance_object(list())    # FALSE

Description

Check if object is a batch assignment result

Usage

is_lap_solve_batch_result(x)

Arguments

Value

Logical indicating if x is a batch assignment result

Description

Check if object is a k-best assignment result

Usage

is_lap_solve_kbest_result(x)

Arguments

Value

Logical indicating if x is a k-best assignment result

Description

Check if object is an assignment result

Usage

is_lap_solve_result(x)

Arguments

Value

Logical indicating if x is an assignment result

Description

Creates an analysis-ready dataset by joining matched pairs with variables from the original left and right datasets. This eliminates the need for manual joins and provides a convenient format for downstream analysis.

Usage

join_matched(result, ...)

## S3 method for class 'matching_result'
join_matched(
  result,
  left,
  right,
  left_vars = NULL,
  right_vars = NULL,
  left_id = "id",
  right_id = "id",
  suffix = c("_left", "_right"),
  include_distance = TRUE,
  include_pair_id = TRUE,
  include_block_id = TRUE,
  ...
)

## S3 method for class 'full_matching_result'
join_matched(result, left, right, left_id = "id", right_id = "id", ...)

## S3 method for class 'cem_result'
join_matched(result, left, right, left_id = "id", right_id = "id", ...)

## S3 method for class 'subclass_result'
join_matched(result, data = NULL, ...)

Arguments

Details

This function simplifies the common workflow of joining matched pairs with original data. Instead of manually merging result$pairs with left and right datasets, join_matched() handles the joins automatically and applies consistent naming conventions.

When variables appear in both left and right datasets, suffixes are appended to distinguish them (e.g., "age_left" and "age_right"). This makes it easy to compute differences or use both values in models.

Value

A tibble with one row per matched pair, containing:

• pair_id: Sequential pair identifier (if include_pair_id = TRUE)

• left_id: ID from left dataset

• right_id: ID from right dataset

• distance: Matching distance (if include_distance = TRUE)

• block_id: Block identifier (if blocking used and include_block_id = TRUE)

• Variables from left dataset (with left suffix)

• Variables from right dataset (with right suffix)

Examples

# Basic usage
left <- data.frame(
  id = 1:5,
  treatment = 1,
  age = c(25, 30, 35, 40, 45),
  income = c(45000, 52000, 48000, 61000, 55000)
)

right <- data.frame(
  id = 6:10,
  treatment = 0,
  age = c(24, 29, 36, 41, 44),
  income = c(46000, 51500, 47500, 60000, 54000)
)

result <- match_couples(left, right, vars = c("age", "income"))
matched_data <- join_matched(result, left, right)
head(matched_data)

# Specify which variables to include
matched_data <- join_matched(
  result, left, right,
  left_vars = c("treatment", "age", "income"),
  right_vars = c("age", "income"),
  suffix = c("_treated", "_control")
)

# Without distance or pair_id
matched_data <- join_matched(
  result, left, right,
  include_distance = FALSE,
  include_pair_id = FALSE
)

Description

Produce an interactive bipartite-graph animation showing how a chosen linear-assignment algorithm transforms the matching over time. The result is an htmlwidget suitable for use in R Markdown, Quarto, pkgdown vignettes, Shiny apps, or standalone HTML output.

Usage

lap_animate(
  cost,
  method = "hungarian",
  maximize = FALSE,
  width = NULL,
  height = NULL,
  elementId = NULL,
  ...
)

Arguments

Details

This is a teaching interface. It runs a slower R reference implementation that emits a state trace at every step, then plays it back in the browser. For production solving, use assignment() or lap_solve() which call the fast C++ backends.

Value

An htmlwidget object.

Animated methods

Animation support is added incrementally. Call animated_methods() to see which method strings currently have a registered trace.

See Also

assignment() for production solving, lap_solve() for the tidy interface.

Examples

## Not run: 
cost <- matrix(c(4, 2, 5,
                 3, 3, 6,
                 7, 5, 4), nrow = 3, byrow = TRUE)
lap_animate(cost, method = "hungarian")

## End(Not run)

Description

Provides a tidy interface for solving the linear assignment problem using 'Hungarian' or 'Jonker-Volgenant' algorithms. Supports rectangular matrices, NA/Inf masking, and data frame inputs.

Usage

lap_solve(
  x,
  source = NULL,
  target = NULL,
  cost = NULL,
  maximize = FALSE,
  method = "auto",
  forbidden = NA
)

Arguments

Value

A tibble with columns:

• source: row/source indices

• target: column/target indices

• cost: cost of each assignment

• total_cost: total cost (attribute)

Examples

# Matrix input
cost <- matrix(c(4, 2, 5, 3, 3, 6, 7, 5, 4), nrow = 3)
lap_solve(cost)

# Data frame input
library(dplyr)
df <- tibble(
  source = rep(1:3, each = 3),
  target = rep(1:3, times = 3),
  cost = c(4, 2, 5, 3, 3, 6, 7, 5, 4)
)
lap_solve(df, source, target, cost)

# With NA masking (forbidden assignments)
cost[1, 3] <- NA
lap_solve(cost)

# Grouped data frames
df <- tibble(
  sim = rep(1:2, each = 9),
  source = rep(1:3, times = 6),
  target = rep(1:3, each = 3, times = 2),
  cost = runif(18, 1, 10)
)
df |> group_by(sim) |> lap_solve(source, target, cost)

Description

Solve many independent assignment problems at once. Supports lists of matrices, 3D arrays, or grouped data frames. Optional parallel execution via n_threads.

Usage

lap_solve_batch(
  x,
  source = NULL,
  target = NULL,
  cost = NULL,
  maximize = FALSE,
  method = "auto",
  n_threads = 1,
  forbidden = NA
)

Arguments

Value

A tibble with columns:

• problem_id: identifier for each problem

• source: source indices for assignments

• target: target indices for assignments

• cost: cost of each assignment

• total_cost: total cost for each problem

• method_used: algorithm used for each problem

Examples

# List of matrices
costs <- list(
  matrix(c(1, 2, 3, 4), 2, 2),
  matrix(c(5, 6, 7, 8), 2, 2)
)
lap_solve_batch(costs)

# 3D array
arr <- array(runif(2 * 2 * 10), dim = c(2, 2, 10))
lap_solve_batch(arr)

# Grouped data frame
library(dplyr)
df <- tibble(
  sim = rep(1:5, each = 9),
  source = rep(1:3, times = 15),
  target = rep(1:3, each = 3, times = 5),
  cost = runif(45, 1, 10)
)
df |> group_by(sim) |> lap_solve_batch(source, target, cost)

# Parallel execution (requires n_threads > 1)
lap_solve_batch(costs, n_threads = 2)

Description

Returns the top k optimal (or near-optimal) assignments using 'Murty' algorithm. Useful for exploring alternative optimal solutions or finding robust assignments.

Usage

lap_solve_kbest(
  x,
  k = 3,
  source = NULL,
  target = NULL,
  cost = NULL,
  maximize = FALSE,
  method = "murty",
  single_method = "jv",
  forbidden = NA
)

Arguments

Value

A tibble with columns:

• rank: ranking of solutions (1 = best, 2 = second best, etc.)

• solution_id: unique identifier for each solution

• source: source indices

• target: target indices

• cost: cost of each edge in the assignment

• total_cost: total cost of the complete solution

Examples

# Matrix input - find 5 best solutions
cost <- matrix(c(4, 2, 5, 3, 3, 6, 7, 5, 4), nrow = 3)
lap_solve_kbest(cost, k = 5)

# Data frame input
library(dplyr)
df <- tibble(
  source = rep(1:3, each = 3),
  target = rep(1:3, times = 3),
  cost = c(4, 2, 5, 3, 3, 6, 7, 5, 4)
)
lap_solve_kbest(df, k = 3, source, target, cost)

# With maximization
lap_solve_kbest(cost, k = 3, maximize = TRUE)

Description

Solves the linear assignment problem when both sources and targets are ordered points on a line. Uses efficient O(n*m) dynamic programming for rectangular problems and O(n) sorting for square problems.

Usage

lap_solve_line_metric(x, y, cost = "L1", maximize = FALSE)

Arguments

Details

This is a specialized solver that exploits the structure of 1-dimensional assignment problems where costs depend only on the distance between points on a line. It is much faster than general LAP solvers for this special case.

The algorithm works as follows:

Square case (n == m): Both vectors are sorted and matched in order: x[1] -> y[1], x[2] -> y[2], etc. This is optimal for any metric cost function on a line.

Rectangular case (n < m): Uses dynamic programming to find the optimal assignment that matches all n sources to a subset of the m targets, minimizing total distance. The DP recurrence is:

dp[i][j] = min(dp[i][j-1], dp[i-1][j-1] + cost(x[i], y[j]))

This finds the minimum cost to match the first i sources to the first j targets.

Complexity:

• Time: O(n*m) for rectangular, O(n log n) for square

• Space: O(n*m) for DP table

Value

A list with components:

• match: Integer vector of length n with 1-based column indices

• total_cost: Total cost of the assignment

Examples

# Square case: equal number of sources and targets
x <- c(1.5, 3.2, 5.1)
y <- c(2.0, 3.0, 5.5)
result <- lap_solve_line_metric(x, y, cost = "L1")
print(result)

# Rectangular case: more targets than sources
x <- c(1.0, 3.0, 5.0)
y <- c(0.5, 2.0, 3.5, 4.5, 6.0)
result <- lap_solve_line_metric(x, y, cost = "L2")
print(result)

# With unsorted inputs (will be sorted internally)
x <- c(5.0, 1.0, 3.0)
y <- c(4.5, 0.5, 6.0, 2.0, 3.5)
result <- lap_solve_line_metric(x, y, cost = "L1")
print(result)

Description

Performs one-to-one matching between two datasets. Supports blocking, distance constraints, and various distance metrics.

Usage

match_couples(
  left,
  right = NULL,
  vars = NULL,
  left_id = NULL,
  right_id = NULL,
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  max_distance = Inf,
  calipers = NULL,
  block_id = NULL,
  ignore_blocks = FALSE,
  require_full_matching = FALSE,
  method = "auto",
  strategy = c("row_best", "sorted", "pq"),
  return_unmatched = TRUE,
  return_diagnostics = FALSE,
  parallel = FALSE,
  replace = FALSE,
  ratio = 1L,
  check_costs = TRUE,
  sigma = NULL,
  memory_mode = "auto",
  certify = NULL
)

Arguments

Details

With method set to a LAP solver (the default "auto", or "jv", "hungarian", ...) it finds the matching that minimizes total distance among all feasible matchings. With method = "greedy" it uses a fast greedy strategy (selected by strategy) that does not guarantee the optimal total distance but scales to very large datasets.

Value

A list with class "matching_result" containing:

• pairs: Tibble of matched pairs with distances

• unmatched: List of unmatched left and right IDs

• info: Matching diagnostics and metadata

• status: One of solver_status_values(), computed from what the solve achieved. "optimal" when every left unit found a partner under an optimal method, "partial" when constraints left some unmatched, "heuristic" when a greedy method ran, either because it was asked for or because the constrained path fell back to it, and "infeasible" when nothing could be matched.

• potentials: The dual potentials of the design's linear program, a list with elements left and right holding one value per unit, named by id. The reduced cost of a pair is its distance minus the two potentials: at least zero on every admissible pair the matching left out, zero on every pair of a 1:1 matching, and at most zero on the pairs of a k:1 or with-replacement matching, whose pair arcs carry at most one unit. A left unit's potential on the k:1 design is the largest of its replicas'. With replacement the LP separates by row: a left unit's potential is the distance to the farthest partner it took and every right unit's is zero. A unit with no admissible partner, or in a block with none, reads NA. Absent for method = "greedy" and for a constrained problem answered by maximum cardinality, which is not the LP the potentials belong to. On the 1:1 design they are the duals verify_assignment() checks, and a method that returns none has them computed by assignment_duals().

Under memory_mode = "implicit" it also carries certificate, the checked statement of optimality (see verify_assignment(), which names the arithmetic it was decided in), and search: the pairs the loop generated out of the pairs the problem states, the pairs a cost was computed for, and one row per round of what each round did. An infeasible answer carries witness instead, naming the units that could not be matched and the partners they have between them.

Examples

# Basic matching
left <- data.frame(id = 1:5, x = c(1, 2, 3, 4, 5), y = c(2, 4, 6, 8, 10))
right <- data.frame(id = 6:10, x = c(1.1, 2.2, 3.1, 4.2, 5.1), y = c(2.1, 4.1, 6.2, 8.1, 10.1))
result <- match_couples(left, right, vars = c("x", "y"))
print(result$pairs)

# With constraints
result <- match_couples(left, right, vars = c("x", "y"),
                        max_distance = 1,
                        calipers = list(x = 0.5))

# With blocking
left$region <- c("A", "A", "B", "B", "B")
right$region <- c("A", "A", "B", "B", "B")
blocks <- matchmaker(left, right, block_type = "group", block_by = "region")
result <- match_couples(blocks$left, blocks$right, vars = c("x", "y"))

# Fast greedy matching for large datasets
result <- match_couples(left, right, vars = c("x", "y"),
                        method = "greedy", strategy = "sorted")

Description

A generic function that converts any couplr matching result into a single analysis-ready data frame with weights, subclass, and distance columns. This is the couplr equivalent of MatchIt's match.data().

Usage

match_data(result, ...)

## S3 method for class 'matching_result'
match_data(result, left, right, left_id = "id", right_id = "id", ...)

## S3 method for class 'full_matching_result'
match_data(result, left, right, left_id = "id", right_id = "id", ...)

## S3 method for class 'cem_result'
match_data(result, left, right, left_id = "id", right_id = "id", ...)

## S3 method for class 'subclass_result'
match_data(result, data = NULL, ...)

Arguments

Details

The output format is compatible with downstream packages like cobalt, WeightIt, and marginaleffects. The stacked (long) format with treatment and weights columns is the standard layout expected by these tools.

Value

A tibble with all original variables plus standardized columns:

id

Unit identifier

treatment

1 for left/treated, 0 for right/control

weights

Matching weights

subclass

Matched group/stratum identifier

distance

Matching distance (where applicable)

Examples

set.seed(42)
left <- data.frame(id = 1:5, age = c(25, 35, 45, 55, 65))
right <- data.frame(id = 6:15, age = runif(10, 20, 70))
result <- match_couples(left, right, vars = "age")
md <- match_data(result, left, right)
head(md)

Description

Solves the matching once for each value of one argument and returns what happened at each: how many units were matched, at what total distance, how balanced the matched sample is, and the certificate saying the matching is optimal for that value. The points are solved as one sequence rather than one at a time, each starting from the matching the previous value found.

Usage

match_path(
  left,
  right,
  vars,
  left_id = NULL,
  right_id = NULL,
  vary = "max_distance",
  values,
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  calipers = NULL,
  sigma = NULL,
  certify = TRUE,
  keep_per_row = .implicit_defaults()$keep_per_row,
  width = .implicit_defaults()$width,
  tol = .implicit_defaults()$tol,
  max_rounds = .implicit_defaults()$max_rounds
)

Arguments

Details

values must ascend. Each point is solved from the point before it, which works because a wider value only ever adds pairs to choose from; a descending sweep takes pairs away from a matching already using them, and is refused rather than reordered.

Value

An object of class couplr_path: ⁠$path⁠, one row per point, ⁠$balance⁠, one row per point per variable, and the match vector, certificate, round record and Hall witness for each of them.

Examples

set.seed(1)
left <- data.frame(id = 1:20, x = rnorm(20), y = rnorm(20))
right <- data.frame(id = 1:60, x = rnorm(60), y = rnorm(60))
path <- match_path(left, right, vars = c("x", "y"),
                   vary = "max_distance", values = c(0.5, 1, 2, Inf))
path$path
path$balance

Description

Constructs blocks (strata) for matching, using either grouping variables or clustering algorithms. Returns the input data frames with block IDs assigned, along with block summary statistics.

Usage

matchmaker(
  left,
  right,
  block_type = c("none", "group", "cluster"),
  block_by = NULL,
  block_vars = NULL,
  block_method = "kmeans",
  n_blocks = NULL,
  min_left = 1,
  min_right = 1,
  drop_imbalanced = FALSE,
  imbalance_threshold = Inf,
  return_dropped = TRUE,
  ...
)

Arguments

Details

This function does NOT perform matching - it only creates the block structure. Use match_couples() to perform matching within blocks.

Value

A list with class "matchmaker_result" containing:

• left: Left data frame with block_id column added

• right: Right data frame with block_id column added

• block_summary: Summary statistics for each block

• dropped: Information about dropped blocks (if any)

• info: Metadata about blocking process

Examples

# Group blocking
left <- data.frame(id = 1:10, region = rep(c("A", "B"), each = 5), x = rnorm(10))
right <- data.frame(id = 11:20, region = rep(c("A", "B"), each = 5), x = rnorm(10))
blocks <- matchmaker(left, right, block_type = "group", block_by = "region")
print(blocks$block_summary)

# Clustering
blocks <- matchmaker(left, right, block_type = "cluster",
                     block_vars = "x", n_blocks = 3)

Description

Computes optimal pixel assignment from A to B and returns the final transported frame (without intermediate animation frames).

Usage

pixel_morph(
  imgA,
  imgB,
  n_frames = 16L,
  mode = c("color_walk", "exact", "recursive", "color_match"),
  lap_method = "jv",
  maximize = FALSE,
  quantize_bits = 5L,
  downscale_steps = 0L,
  alpha = 1,
  beta = 0,
  patch_size = 1L,
  upscale = 1,
  show = interactive()
)

Arguments

Details

Transport-Only Semantics

This function returns a SHARP, pixel-perfect transport of A's pixels to positions determined by the assignment to B.

Key Points:

• Assignment computed using: cost = alpha * color_dist + beta * spatial_dist

• B's COLORS influence assignment but DO NOT appear in output

• Result has A's colors arranged to match B's layout

• No motion blur (unlike intermediate frames in animation)

See pixel_morph_animate for detailed explanation of assignment vs rendering semantics.

Permutation Warnings

Assignment is guaranteed to be a bijection (permutation) ONLY when:

• downscale_steps = 0 (no resolution changes)

• mode = "exact" with patch_size = 1

With downscaling or tiled modes, assignment may have:

• Overlaps: Multiple source pixels map to same destination (last write wins)

• Holes: Some destinations never filled (remain transparent)

If assignment is not a bijection (due to downscaling or tiling), a warning will be issued. The result may contain:

• Overlapped pixels (multiple sources -> one destination)

• Transparent holes (some destinations unfilled)

For guaranteed pixel-perfect results, use:

  pixel_morph(A, B, mode = "exact", downscale_steps = 0)

Value

magick image object of the final transported frame

See Also

pixel_morph_animate for animated version

Examples

if (requireNamespace("magick", quietly = TRUE)) {
  imgA <- system.file("extdata/icons/circleA_40.png", package = "couplr")
  imgB <- system.file("extdata/icons/circleB_40.png", package = "couplr")
  if (nzchar(imgA) && nzchar(imgB)) {
    result <- pixel_morph(imgA, imgB, n_frames = 4, show = FALSE)
  }
}

Description

Creates an animated morph by computing optimal pixel assignment from image A to image B, then rendering intermediate frames showing the transport.

Usage

pixel_morph_animate(
  imgA,
  imgB,
  n_frames = 16L,
  fps = 10L,
  format = c("gif", "webp", "mp4"),
  outfile = NULL,
  show = interactive(),
  mode = c("color_walk", "exact", "recursive", "color_match"),
  lap_method = "jv",
  maximize = FALSE,
  quantize_bits = 5L,
  downscale_steps = 0L,
  alpha = 1,
  beta = 0,
  patch_size = 1L,
  upscale = 1
)

Arguments

Details

Assignment vs Rendering Semantics

CRITICAL: This function has two separate phases with different semantics:

Phase 1 - Assignment Computation:

The assignment is computed by minimizing:

  cost(i,j) = alpha * color_distance(A[i], B[j]) + 
              beta * spatial_distance(pos_i, pos_j)

This means B's COLORS influence which pixels from A map to which positions.

Phase 2 - Rendering (Transport-Only):

The renderer uses ONLY A's colors:

• Intermediate frames: A's pixels move along paths with motion blur

• Final frame: A's pixels at their assigned positions (sharp, no blur)

• B's colors NEVER appear in the output

Result: You get A's colors rearranged to match B's geometry/layout.

What This Means

• B influences WHERE pixels go (via similarity in cost function)

• B does NOT determine WHAT COLORS appear in output

• Final image has A's palette arranged to mimic B's structure

Parameter Guidance

For pure spatial rearrangement (ignore B's colors in assignment):

  pixel_morph_animate(A, B, alpha = 0, beta = 1)

For color-similarity matching (default):

  pixel_morph_animate(A, B, alpha = 1, beta = 0)

For hybrid (color + spatial):

  pixel_morph_animate(A, B, alpha = 1, beta = 0.2)

Permutation Guarantees

Assignment is guaranteed to be a bijection (permutation) ONLY when:

• downscale_steps = 0 (no resolution changes)

• mode = "exact" with patch_size = 1

With downscaling or tiled modes, assignment may have:

• Overlaps: Multiple source pixels map to same destination (last write wins)

• Holes: Some destinations never filled (remain transparent)

A warning is issued if overlaps/holes are detected in the final frame.

Value

Invisibly returns a list with animation object and metadata:

Examples

if (requireNamespace("magick", quietly = TRUE)) {
  imgA <- system.file("extdata/icons/circleA_40.png", package = "couplr")
  imgB <- system.file("extdata/icons/circleB_40.png", package = "couplr")
  if (nzchar(imgA) && nzchar(imgB)) {
    outfile <- tempfile(fileext = ".gif")
    pixel_morph_animate(imgA, imgB, outfile = outfile, n_frames = 4, show = FALSE)
  }
}

Description

Produces a Love plot (dot plot) of standardized differences.

Usage

## S3 method for class 'balance_diagnostics'
plot(x, type = c("love", "histogram", "variance"), threshold = 0.1, ...)

Arguments

Value

The balance_diagnostics object (invisibly)

Description

Produces a histogram of pairwise distances from a matching result.

Usage

## S3 method for class 'matching_result'
plot(x, type = c("histogram", "density", "ecdf"), ...)

Arguments

Value

The matching_result object (invisibly)

Description

Plot method for sensitivity analysis (base graphics)

Usage

## S3 method for class 'sensitivity_analysis'
plot(x, alpha = 0.05, ...)

Arguments

Value

The sensitivity_analysis object (invisibly)

Description

Main preprocessing function that orchestrates variable health checks, categorical encoding, and automatic scaling selection.

Usage

preprocess_matching_vars(
  left,
  right,
  vars,
  auto_scale = TRUE,
  scale_method = "auto",
  check_health = TRUE,
  remove_problematic = TRUE,
  verbose = TRUE
)

Arguments

Value

A list with class "preprocessing_result" containing:

• left: Preprocessed left data frame

• right: Preprocessed right data frame

• vars: Final variable names (after exclusions)

• health: Variable health diagnostics

• scaling_method: Selected scaling method

• excluded_vars: Variables that were excluded

• warnings: List of warnings issued

Description

Print Method for Balance Diagnostics

Usage

## S3 method for class 'balance_diagnostics'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print Method for CEM Results

Usage

## S3 method for class 'cem_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Shows how many pairs the whole sweep generated out of how many exist, and then the point table: one row per value, with the matching it found and what finding it cost.

Usage

## S3 method for class 'couplr_path'
print(x, ...)

Arguments

Value

x, invisibly.

Description

Print Method for Distance Objects

Usage

## S3 method for class 'distance_object'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print Method for Full Matching Results

Usage

## S3 method for class 'full_matching_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Prints a summary and the table of results for a batch of assignment problems solved with lap_solve_batch().

Usage

## S3 method for class 'lap_solve_batch_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for k-best assignment results

Usage

## S3 method for class 'lap_solve_kbest_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Nicely prints a lap_solve_result object, including the assignments, total cost, and method used.

Usage

## S3 method for class 'lap_solve_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for matching results

Usage

## S3 method for class 'matching_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for matchmaker results

Usage

## S3 method for class 'matchmaker_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for preprocessing result

Usage

## S3 method for class 'preprocessing_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for sensitivity analysis

Usage

## S3 method for class 'sensitivity_analysis'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print Method for Subclassification Results

Usage

## S3 method for class 'subclass_result'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Print method for variable health

Usage

## S3 method for class 'variable_health'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

Description

Matches treated and control units based on estimated propensity scores. Fits a logistic regression model (or accepts a pre-fitted one), computes logit propensity scores, and calls match_couples() with a caliper on the logit scale.

Usage

ps_match(
  formula = NULL,
  data = NULL,
  treatment = NULL,
  ps_model = NULL,
  caliper_sd = 0.2,
  method = "auto",
  replace = FALSE,
  ratio = 1L,
  ...
)

Arguments

Details

The propensity score is the probability of treatment assignment conditional on observed covariates. Matching is performed on the logit of the propensity score (Rosenbaum and Rubin 1985), which provides better distributional properties than matching on the raw probability scale.

The default caliper of 0.2 SD of logit(PS) is recommended by Austin (2011) as removing approximately 98% of bias.

Value

A matching_result object with additional propensity score info in result$info$ps_model and result$info$caliper_value.

Examples

set.seed(42)
n <- 100
data <- data.frame(
  id = seq_len(n),
  treated = rbinom(n, 1, 0.4),
  age = rnorm(n, 50, 10),
  income = rnorm(n, 50000, 15000)
)
result <- ps_match(treated ~ age + income, data = data, treatment = "treated")
print(result)

Description

Assesses how sensitive a matched comparison is to hidden bias using Rosenbaum bounds on the Wilcoxon signed-rank statistic.

Usage

sensitivity_analysis(
  result,
  left,
  right,
  outcome_var,
  gamma = seq(1, 3, by = 0.25),
  alternative = c("greater", "less", "two.sided"),
  left_id = "id",
  right_id = "id"
)

Arguments