Package {couplr}

Title:	Optimal Pairing and Matching via Linear Assignment
Version:	1.6.1
Description:	Solves optimal pairing and matching problems using linear assignment algorithms. Provides implementations of the Hungarian method (Kuhn 1955) <doi:10.1002/nav.3800020109>, Jonker-Volgenant shortest path algorithm (Jonker and Volgenant 1987) <doi:10.1007/BF02278710>, Auction algorithm (Bertsekas 1988) <doi:10.1007/BF02186476>, cost-scaling (Goldberg and Kennedy 1995) <doi:10.1007/BF01585996>, scaling algorithms (Gabow and Tarjan 1989) <doi:10.1137/0218069>, push-relabel (Goldberg and Tarjan 1988) <doi:10.1145/48014.61051>, and Sinkhorn entropy-regularized transport (Cuturi 2013) <doi:10.48550/arxiv.1306.0895>. Designed for matching plots, sites, samples, or any pairwise optimization problem. Supports rectangular matrices, forbidden assignments, data frame inputs, batch solving, k-best solutions, and pixel-level image morphing for visualization. Includes automatic preprocessing with variable health checks, multiple scaling methods (standardized, range, robust), greedy matching algorithms, and comprehensive balance diagnostics for assessing match quality using standardized differences and distribution comparisons.
License:	MIT + file LICENSE
Language:	en-US
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	Rcpp (≥ 1.0.0), tibble (≥ 3.0.0), dplyr (≥ 1.0.0), rlang (≥ 0.4.0), generics (≥ 0.1.0), purrr (≥ 0.3.0)
Suggests:	testthat (≥ 3.0.0), htmlwidgets, e1071, R.utils, microbenchmark, withr, knitr, rmarkdown, bench, parallel, future (≥ 1.20.0), future.apply (≥ 1.8.0), parallelly (≥ 1.30.0), ggplot2, ggraph, tidygraph, magick, av, png, combinat, cobalt, MatchIt, marginaleffects, optmatch
LinkingTo:	Rcpp
SystemRequirements:	C++17
LazyData:	true
VignetteBuilder:	knitr
URL:	https://gillescolling.com/couplr/, https://github.com/gcol33/couplr
BugReports:	https://github.com/gcol33/couplr/issues
Config/testthat/edition:	3
Config/roxygen2/version:	8.0.0
Config/Needs/coverage:	xml2
NeedsCompilation:	yes
Packaged:	2026-08-22 23:04:39 UTC; Gilles Colling
Author:	Gilles Colling [aut, cre, cph]
Maintainer:	Gilles Colling <gilles.colling051@gmail.com>
Repository:	CRAN
Date/Publication:	2026-08-23 00:20:02 UTC

couplr: Optimal Pairing and Matching via Linear Assignment

Description

Solves optimal pairing and matching problems using linear assignment algorithms. Provides implementations of the Hungarian method (Kuhn 1955) doi:10.1002/nav.3800020109, Jonker-Volgenant shortest path algorithm (Jonker and Volgenant 1987) doi:10.1007/BF02278710, Auction algorithm (Bertsekas 1988) doi:10.1007/BF02186476, cost-scaling (Goldberg and Kennedy 1995) doi:10.1007/BF01585996, scaling algorithms (Gabow and Tarjan 1989) doi:10.1137/0218069, push-relabel (Goldberg and Tarjan 1988) doi:10.1145/48014.61051, and Sinkhorn entropy-regularized transport (Cuturi 2013) doi:10.48550/arxiv.1306.0895. Designed for matching plots, sites, samples, or any pairwise optimization problem. Supports rectangular matrices, forbidden assignments, data frame inputs, batch solving, k-best solutions, and pixel-level image morphing for visualization. Includes automatic preprocessing with variable health checks, multiple scaling methods (standardized, range, robust), greedy matching algorithms, and comprehensive balance diagnostics for assessing match quality using standardized differences and distribution comparisons.

Solves optimal pairing and matching problems using linear assignment algorithms. Designed for matching plots, sites, samples, or any pairwise optimization problem. Provides modern, tidy implementations of 'Hungarian', 'Jonker-Volgenant', 'Auction', and other LAP solvers.

Main functions

• lap_solve: Solve single assignment problems

• lap_solve_batch: Solve multiple problems efficiently

• lap_solve_kbest: Find k-best optimal solutions

Author(s)

Maintainer: Gilles Colling gilles.colling051@gmail.com (ORCID) [copyright holder]

Authors:

• Gilles Colling gilles.colling051@gmail.com (ORCID) [copyright holder]

See Also

Useful links:

• https://gillescolling.com/couplr/

• https://github.com/gcol33/couplr

• Report bugs at https://github.com/gcol33/couplr/issues

Solve an implicit assignment by generating the pairs it needs

Description

The backend for memory_mode = "implicit". Takes either a lazy_cost_spec, which states the complete problem without building it, or a materialized matrix, which states the same problem and lets the loop's answer be held against a complete solve of the same numbers.

Usage

.assignment_implicit(
  cost,
  maximize = FALSE,
  certify = TRUE,
  method = "auto",
  keep_per_row = .implicit_defaults()$keep_per_row,
  width = .implicit_defaults()$width,
  tol = .implicit_defaults()$tol,
  max_rounds = .implicit_defaults()$max_rounds
)

Value

A lap_solve_result carrying the duals the last master produced, the certificate for the complete problem, and search, the record of what the search cost.

Solve a lazy_cost_spec (memory_mode = "lazy" backend for assignment())

Description

Mirrors assignment()'s contract (same result shape) but computes costs on demand from the underlying feature data instead of a materialized matrix. Only "jv" and "auction" are supported – every other method is fundamentally dense (repeated full-matrix scans, or an algorithm not yet templated for a lazy cost source) and gets a clear error here rather than a silent dense fallback that would defeat the point of memory_mode = "lazy".

Usage

.assignment_lazy(cost, maximize = FALSE, method = "auto", auction_eps = NULL)

Compute automatic break points for CEM binning

Description

Compute automatic break points for CEM binning

Usage

.auto_breaks(x, n_bins)

Arguments

Value

Numeric vector of break points

Histogram of |std diff| via ggplot2

Description

Histogram of |std diff| via ggplot2

Usage

.autoplot_hist(object, threshold = 0.1)

Love plot via ggplot2

Description

Love plot via ggplot2

Usage

.autoplot_love(object, threshold = 0.1)

Variance ratio plot via ggplot2

Description

Variance ratio plot via ggplot2

Usage

.autoplot_variance(object)

Check a balance flow against the objective it is meant to encode

Description

Recomputes the matched sample's cardinality, its imbalance at every level and its total distance from the flow's own pair arcs, and compares the arc costs against P * (K - pairs) + sum_h eps_h * imbalance_h + sum_ij d_ij x_ij.

Usage

.balance_flow_audit(problem, index, flow, cost = NULL, tiers = index$tiers)

Value

A list holding the feasibility counts, the objective read both ways, and the gap between them.

The flow a matched set corresponds to

Description

Builds the flow vector a candidate matched set maps to: unit and pair arcs at one, the slack each category needs to fill its budget, and the transfers that carry its imbalance, each crossing at the lowest level its two cells share.

Usage

.balance_flow_encode(matching, index, hier = index$hier)

Details

Returns NULL when the matched set is not balanced at the levels the design enforces exactly, since no flow in this network represents it.

Compile a balance design into a flow problem

Description

Compile a balance design into a flow problem

Usage

.balance_flow_problem(
  cost,
  hier,
  codes = NULL,
  tiers = NULL,
  arc_bounds = NULL
)

Value

A list with problem, a couplr_flow_problem, and index, holding the arc ranges and the units and cells behind each arc.

Read a solved balance flow

Description

Read a solved balance flow

Usage

.balance_flow_read(index, flow)

Value

A list with the matched pairs as left and right indices, the slack each finest category carried, and the transfers with the level they crossed at.

Node layout of a balance flow network

Description

Every base offset the network uses, and the accessors that turn a level and a category into a node id. No other function computes a node id.

Usage

.balance_node_layout(n_left, n_right, hierarchy)

Shared parallel block matching implementation

Description

Core logic for both optimal (LAP) and greedy parallel block matching. Called by match_blocks_parallel() and greedy_blocks_parallel().

Usage

.blocks_parallel(
  blocks,
  left,
  right,
  left_ids,
  right_ids,
  block_col,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  solver_fn,
  solver_params = list(),
  check_costs = FALSE,
  strict_no_pairs = FALSE,
  parallel = FALSE,
  replace = FALSE,
  ratio = 1L,
  sigma = NULL,
  memory_mode = "auto"
)

Arguments

Value

List with combined results from all blocks

Branch and bound over the moment constraints

Description

Searches the tree of inclusion decisions, bounding every node by its Lagrangian relaxation and keeping an incumbent that satisfies every stated constraint. Returns at any interruption with that incumbent and a global bound that is valid for the whole problem, never with an unproven claim of optimality.

Usage

.cardinality_branch_bound(
  problem,
  index = NULL,
  coefs = NULL,
  dual_steps = 20L,
  branch = c("unit", "pair"),
  node_limit = 500L,
  time_limit = Inf,
  should_stop = NULL,
  cost = NULL,
  tol = 1e-09
)

Arguments

Details

time_limit reaches the solver. A solve that runs out of budget stops between augmentations and comes back saying so, and the node it belonged to goes back on the frontier unopened, so the bound the search reports still covers the whole tree.

Value

A list of class cardinality_run.

One solve of a balance network

Description

Applies a node's arc bounds and multipliers, solves the network, certifies the flow against the problem it was solved on, audits it against the objective identity the design encodes, and reads the matched set back.

Usage

.cardinality_flow(
  problem,
  index = NULL,
  coefs = NULL,
  lambda = NULL,
  edits = NULL,
  cost = NULL,
  tol = 1e-09,
  warm = NULL,
  time_limit = Inf
)

Arguments

Value

A list with the solve status, the flow and potentials, the certificate, the audit, the matched set, the true objective objective and the relaxed objective relaxed the multipliers price. A solve that ran out of time comes back with status "interrupted", its flow and potentials, and nothing else: it proved neither an optimum nor the absence of one, so certifying and auditing it would be work spent on a number no caller may read.

Report for a pruned match

Description

The fields a cardinality_report carries that the prune loop computes, and NA for the ones it does not. The loop maximizes nothing and derives no bound, so best_possible and gap have no value to hold; NA is what they hold, rather than a number that would read as a bound.

Usage

.cardinality_heuristic_report(pairs, max_std_diff)

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

Arguments

Value

An object of class c("cardinality_heuristic", "cardinality_report").

Invisibly returns x.

Lagrangian bound for one node

Description

Solves the node's network at a sequence of multipliers, each step one .flow_solve() and one verify_flow(), and returns the largest certified relaxed optimum it reached. Multipliers move along the projected subgradient with step t_0 / (1 + k), warm-started from whatever lambda is handed in, which is the parent's best set during a search.

Usage

.cardinality_lagrangian(
  problem,
  coefs = NULL,
  lambda = NULL,
  steps = 20L,
  index = NULL,
  edits = NULL,
  incumbent = Inf,
  step0 = NULL,
  tol = 1e-09,
  cost = NULL,
  warm = NULL,
  deadline = Inf
)

Arguments

Details

Consecutive steps solve one network. The topology, the arc bounds and every non-pair cost are the ones the previous step solved; only the pair costs moved, by one multiplier step, so the previous step's flow and potentials are a near-optimal starting point for this one and are carried across. warm gives the first step the same footing, from whichever solve the caller has to hand.

Value

A list with bound, the multipliers that attained it, the relaxed solve at those multipliers, any moment-feasible solutions the ascent passed through, certified, whether the bound it reports came out of a certified solve, and warm, the last complete solve's flow and potentials. A status of "interrupted" means a solve ran out of time; the bound and the solutions reported alongside it came from the steps that finished, and are as valid as any others.

Balance pruning heuristic

Description

Runs a full optimal match, then deletes pairs until every matching variable sits inside the standardized-difference threshold or the iteration budget runs out.

Usage

.cardinality_prune(
  left,
  right,
  vars,
  max_std_diff,
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  auto_scale = FALSE,
  sigma = NULL,
  max_distance = Inf,
  calipers = NULL,
  left_id = NULL,
  right_id = NULL,
  method = "auto",
  max_iter = 100L,
  batch_fraction = 0.1
)

Arguments

Value

A matching_result carrying info$pruning_iterations and info$pairs_removed.

Report of a cardinality match

Description

Turns a completed search into the object callers read: the matched sample's size, the largest size the bound admits, the gap between them, and the state of every constraint the match was asked to meet.

Usage

.cardinality_report(run, specs = NULL, tol = 1e-09)

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

Arguments

Value

An object of class cardinality_report, a list with elements:

• n_matched, n_left_matched - matched pairs, and the left units they use.

• best_possible - the largest matched sample the bound admits.

• gap, gap_fraction - best_possible - n_matched, in matched units and as a share of best_possible.

• certified - TRUE only when the search settled, every solve its bound rests on was certified, the incumbent's own solve was certified and audited, and gap is zero.

• objective, bound - the incumbent's value of the network objective and the global lower bound on it. best_possible is read from bound.

• stopped_on, n_nodes, engine, status.

• constraints - one row per stated constraint, with what it asked for and what the matched sample achieved.

• balance - matched counts and imbalance per category at every level.

• tiers, precision_headroom, shift - the weights that order the objective, how many times its range fits inside the range a double orders exactly, and the constant taken off every distance.

• total_distance, pairs - the matched set itself.

Invisibly returns x.

Assemble a matching_result from a cardinality report

Description

Assemble a matching_result from a cardinality report

Usage

.cardinality_result(
  report,
  left,
  right,
  vars,
  left_ids,
  right_ids,
  distance = "euclidean",
  max_std_diff = NULL
)

Arguments

Value

A matching_result carrying cardinality, status, info$engine, and a certificate when the search certified optimality.

Balance-constrained cardinality match

Description

Builds the balance network for a distance matrix and a nested partition, searches it under whatever moment constraints are stated, and reports the matched sample with the bound that goes with it.

Usage

.cardinality_solve(
  left,
  right,
  cost,
  refined = NULL,
  exact = 1L,
  moments = NULL,
  max_std_diff = NULL,
  vars = NULL,
  dual_steps = 20L,
  branch = c("unit", "pair"),
  node_limit = 500L,
  time_limit = Inf,
  should_stop = NULL,
  tol = 1e-09
)

Arguments

Value

A cardinality_report.

Compute standardized differences for current pairs

Description

Compute standardized differences for current pairs

Usage

.compute_pair_balance(pairs, left, right, vars, left_ids, right_ids)

Shared blocked matching implementation

Description

Core logic for both optimal (LAP) and greedy blocked matching. Called by match_couples_blocked() (optimal and greedy paths).

Usage

.couples_blocked(
  left,
  right,
  left_ids,
  right_ids,
  block_col,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  solver_fn,
  solver_params = list(),
  check_costs = FALSE,
  strict_no_pairs = FALSE,
  parallel = FALSE,
  replace = FALSE,
  ratio = 1L,
  sigma = NULL,
  memory_mode = "auto"
)

Arguments

Value

List with pairs tibble, unmatched list, and info list.

Shared matching from precomputed distance object

Description

Core logic for both optimal (LAP) and greedy matching from distance objects. Called by match_couples_from_distance() (optimal and greedy paths).

Usage

.couples_from_distance(
  dist_obj,
  max_distance = Inf,
  calipers = NULL,
  ignore_blocks = FALSE,
  require_full_matching = FALSE,
  return_unmatched = TRUE,
  return_diagnostics = FALSE,
  solver_fn,
  solver_params = list(),
  check_costs = FALSE,
  strict_no_pairs = FALSE,
  method_label = "from_distance_object",
  extra_info = list(),
  diagnostics_fields = c("method", "n_matched", "total_distance")
)

Arguments

Value

A matching_result object with pairs, info, and optional diagnostics.

Replacement matching: each left picks its best right independently

Description

The compiled design gives every column capacity for every row, so the rows never compete and the optimum of the whole network is each row's own cheapest columns. plan$per_row is how many of them a row takes: the requested ratio, or the column count when there are fewer columns than that.

Usage

.couples_replace(
  cost_matrix,
  left,
  right,
  left_ids,
  right_ids,
  vars,
  ratio = 1L,
  plan
)

Value

List with pairs tibble, unmatched list, and info list.

Shared single matching implementation

Description

Core logic for both optimal (LAP) and greedy matching without blocking. Called by match_couples_single() (optimal and greedy paths).

Usage

.couples_single(
  left,
  right,
  left_ids,
  right_ids,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  solver_fn,
  solver_params = list(),
  check_costs = FALSE,
  strict_no_pairs = FALSE,
  replace = FALSE,
  ratio = 1L,
  sigma = NULL,
  memory_mode = "auto"
)

Arguments

Value

List with pairs tibble, unmatched list, and info list.

Summarise a per-variable balance table

Description

The three headline numbers every balance object reports, read off the per-variable table each of them builds.

Usage

.overall_balance(var_stats, n_vars)

Arguments

Value

List with mean_abs_std_diff, max_abs_std_diff, pct_large_imbalance and n_vars.

Compute per-pair differences on a single variable

Description

Compute per-pair differences on a single variable

Usage

.pair_var_diffs(pairs, left, right, var, left_ids, right_ids)

Compute Rosenbaum bounds via normal approximation

Description

Under hidden bias Gamma, each pair's probability of the treated unit having the larger value is between 1/(1+Gamma) and Gamma/(1+Gamma).

Usage

.rosenbaum_bounds(t_obs, ranks, n, gamma, alternative)

Value

List with p_upper and p_lower.

Large value for forbidden pairs

Description

A numeric constant used to mark forbidden pairs in cost matrices.

Usage

BIG_COST

Format

Numeric value (half of .Machine$double.xmax).

List methods that currently support animation

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.

Apply all constraints to cost matrix

Description

Main entry point for applying constraints.

Usage

apply_all_constraints(
  cost_matrix,
  left,
  right,
  vars,
  max_distance = Inf,
  calipers = NULL,
  forbidden = NULL
)

Value

Modified cost matrix with all constraints applied.

Apply caliper constraints

Description

Calipers impose per-variable maximum absolute differences.

Usage

apply_calipers(cost_matrix, left, right, calipers, vars)

Value

Modified cost matrix with forbidden pairs marked.

Apply maximum distance constraint

Description

Apply maximum distance constraint

Usage

apply_max_distance(cost_matrix, max_distance = Inf)

Value

Modified cost matrix with forbidden pairs marked.

Apply scaling to matching variables

Description

Apply scaling to matching variables

Usage

apply_scaling(left_mat, right_mat, method = "standardize")

Value

List with scaled left/right matrices and scaling parameters.

Apply weights to matching variables

Description

Apply weights to matching variables

Usage

apply_weights(mat, weights)

Value

Numeric matrix with columns weighted.

Convert assignment result to a binary matrix

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

Convert couplr Result to matchit Object

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)

Assign blocks using clustering

Description

Assign blocks using clustering

Usage

assign_blocks_cluster(left, right, block_vars, method, n_blocks, ...)

Value

List with modified left/right data frames (with block_id) and n_blocks_initial.

Assign blocks based on grouping variable(s)

Description

Assign blocks based on grouping variable(s)

Usage

assign_blocks_group(left, right, block_by)

Value

List with modified left/right data frames (with block_id) and n_blocks_initial.

Linear assignment solver

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", "csflow", "hk01", "bruteforce", "ssap_bucket",
    "cycle_cancel", "gabow_tarjan", "lapmod", "csa", "ramshaw_tarjan", "push_relabel",
    "orlin", "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/shape and data characteristics:

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

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

• Sparse (>50\

• Very rectangular (m >= 3n): "sap" — handles rectangular well

• Otherwise: "jv" — fastest general-purpose solver at every size

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.

Under memory_mode = "implicit" the result also carries u and v, the duals the last restricted master produced, and search: 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.

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

Solve assignment problem and return dual variables

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))

Augment Matching Results with Original Data (broom-style)

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)

Automatically encode categorical variables

Description

Converts categorical variables to numeric representations suitable for matching. Currently supports binary variables (0/1) and ordered factors.

Usage

auto_encode_categorical(left, right, var)

Arguments

Value

List with encoded left and right columns, plus encoding metadata

ggplot2 autoplot for balance diagnostics

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)
}

ggplot2 autoplot for matching results

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")
}

ggplot2 autoplot for sensitivity analysis

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

Balance Table for Matching Results (cobalt integration)

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

Balance Diagnostics for Matched Pairs

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)

Create Balance Table

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

Solve the Bottleneck Assignment Problem

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)

Build cost matrix for matching

Description

This is the main entry point for distance computation.

Usage

build_cost_matrix(
  left,
  right,
  vars,
  distance = "euclidean",
  weights = NULL,
  scale = FALSE,
  sigma = NULL,
  memory_mode = "auto",
  caller_supports_lazy = TRUE,
  caller_supports_implicit = caller_supports_lazy
)

Arguments

Value

Numeric matrix of distances with optional scaling/weights applied.

Calculate Variable-Level Balance Statistics

Description

Calculate Variable-Level Balance Statistics

Usage

calculate_var_balance(
  left_vals,
  right_vals,
  var_name,
  w_left = NULL,
  w_right = NULL
)

Arguments

Value

List with balance statistics for this variable

Check if parallel processing is available

Description

Check if parallel processing is available

Usage

can_parallelize()

Value

Logical indicating if future package is available

Cardinality Matching

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
)

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.

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.

• 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

Coarsened Exact Matching

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)

Check cost distribution for problems

Description

Examines the distance matrix for common issues and provides helpful warnings.

Usage

check_cost_distribution(cost_matrix, threshold_zero = 1e-10, warn = TRUE)

Arguments

Value

List with diagnostic information

Check if full matching was achieved

Description

Check if full matching was achieved

Usage

check_full_matching(result)

Value

No return value; throws error if unmatched units exist.

Error on repeated id values

Description

Error on repeated id values

Usage

check_unique_ids(ids, id_col, side)

Value

The ids, invisibly unchanged, when they are unique.

Check variable health for matching

Description

Analyzes variables for common problems that can affect matching quality: constant columns, high missingness, extreme skewness, and outliers.

Usage

check_variable_health(
  left,
  right,
  vars,
  high_missingness_threshold = 0.5,
  low_variance_threshold = 1e-06
)

Arguments

Value

A list with class "variable_health" containing:

• summary: Tibble with per-variable diagnostics

• issues: List of detected issues with severity levels

• exclude_vars: Variables that should be excluded

• warnings: Human-readable warnings

Compute pairwise distance matrix

Description

Compute pairwise distance matrix

Usage

compute_distance_matrix(
  left_mat,
  right_mat,
  distance = "euclidean",
  sigma = NULL
)

Arguments

Value

Numeric matrix of pairwise distances (n_left x n_right).

Compute and Cache Distance Matrix for Reuse

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

Count valid pairs in cost matrix

Description

Count valid pairs in cost matrix

Usage

count_valid_pairs(cost_matrix)

Value

Integer count of valid (non-forbidden) pairs.

Get a themed emoji

Description

Get a themed emoji

Usage

couplr_emoji(
  type = c("error", "warning", "info", "success", "heart", "broken", "sparkles",
    "search", "chart", "warning_sign", "stop", "check")
)

Value

Character string with the emoji (or empty string if emoji disabled).

Info message with emoji

Description

Info message with emoji

Usage

couplr_inform(...)

Value

No return value, called for side effects (issues a message).

Couplr message helpers with emoji and humor

Description

Light, fun error/warning messages inspired by testthat, themed around coupling and matching. Makes errors less intimidating and more memorable.

Stop with a fun, themed error message

Description

Stop with a fun, themed error message

Usage

couplr_stop(..., call. = FALSE)

Value

No return value, throws an error.

Success message with emoji

Description

Success message with emoji

Usage

couplr_success(...)

Value

No return value, called for side effects (issues a message).

Warn with a fun, themed warning message

Description

Warn with a fun, themed warning message

Usage

couplr_warn(..., call. = FALSE)

Value

No return value, called for side effects (issues a warning).

What the design identifies

Description

Every matching design in the package weights the left side at 1 and the right side to reproduce the left distribution within a pair, group or stratum, so what it targets is the effect on the left population: the ATT, in the orientation the matching layer is written in, where left holds the treated units.

Usage

design_estimand(n_left, n_left_matched)

Details

focal_discarded is how many left units the design did not retain. Calipers, max_distance and coarsened strata all drop focal units, and the estimate then speaks about the retained ones rather than about the whole treated group, which is what as_matchit() reports when it hands the estimand to MatchIt and marginaleffects.

Value

Named list of fields to merge into a result's info.

Detect and validate blocking

Description

Detect and validate blocking

Usage

detect_blocking(left, right, block_id, ignore_blocks)

Value

List with use_blocking (logical) and block_col (character or NULL).

Diagnose distance matrix and suggest fixes

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

Invalid parameter error

Description

Invalid parameter error

Usage

err_invalid_param(param, value, expected)

Value

No return value, throws an error.

Missing data error

Description

Missing data error

Usage

err_missing_data(dataset = "left")

Value

No return value, throws an error.

Missing variables error

Description

Missing variables error

Usage

err_missing_vars(vars, dataset = "left")

Value

No return value, throws an error.

All pairs forbidden error

Description

All pairs forbidden error

Usage

err_no_valid_pairs(reason = NULL)

Value

No return value, throws an error.

Estimate dense cost-matrix memory footprint in megabytes

Description

The true peak footprint is not 8 bytes/cell: matrix(0, n, m) at the R level (8B) is followed by rcpp_to_cost_matrix()'s copy into a lap::CostMatrix (8B data + 4B mask), then prepare_for_solve()'s unconditional copy (another 12B), and a possible t() transpose copy (8B). Several of these are transient but can coexist during GC lag. A conservative multiplier avoids systematically under-warning; 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)

Example cost matrices for assignment problems

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)

Example assignment problem data frame

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)
  )

Explain which solver method = "auto" selects, and why

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))

Extract and standardize IDs from data frames

Description

The id a matching carries is the key every downstream verb joins on, so it is resolved once, here, and the same resolution answers match_couples() and join_matched(). id_col names the column to read; with no name given a column called id is used, then meaningful row names, then ids synthesized from prefix.

Usage

extract_ids(df, prefix = "id", id_col = NULL, warn_synthetic = FALSE)

Arguments

Details

Ids read from the data must be unique: a repeated value makes every id-keyed join downstream many-to-many, which expands rows and attaches one unit's covariates to another unit's pair. Synthesized ids are unique by construction and are not checked.

Value

Character vector of IDs.

Extract matching variables from data frame

Description

Extract matching variables from data frame

Usage

extract_matching_vars(df, vars)

Value

Numeric matrix of matching variables.

Filter blocks based on size and balance criteria

Description

Filter blocks based on size and balance criteria

Usage

filter_blocks(
  left,
  right,
  min_left,
  min_right,
  drop_imbalanced,
  imbalance_threshold
)

Value

List with filtered left/right data frames and dropped block info.

Full Matching

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 matching creates matched groups of variable size. Two algorithms are available:

Optimal (method = "optimal", default): Solves a min-cost max-flow problem that minimizes total distance across all group assignments simultaneously. Each left unit becomes a group center absorbing 1 to max_controls right units, with the globally optimal assignment found via Dijkstra's algorithm with Johnson potentials. When n_left > n_right, roles are transposed automatically.

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.

Weights are computed so that within each group, the total weight of right units equals the total weight of left units (which is 1). For a group with 1 left and k right units, each right unit receives weight 1/k.

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 proven. Present for method = "optimal" only.

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)

Standardize block ID column name

Description

Standardize block ID column name

Usage

get_block_id_column(df)

Value

Character string with column name, or NULL if not found.

Estimate available system RAM in megabytes

Description

Cross-platform, base-R-only (shells out; no new package dependency). Never errors: returns NA_real_ if detection fails or the platform is unrecognized, so callers must treat NA as "unknown" and fall back to a fixed threshold rather than skipping the guard entirely.

Usage

get_free_ram_mb()

Details

"Available" means memory an allocation can obtain without swapping, which on every platform is more than the untouched free list: Linux reports it directly as MemAvailable, and on macOS it is the free, inactive and speculative pages together, since the kernel keeps almost nothing on the free list and reclaims the rest on demand.

Value

Numeric scalar (MB of available RAM), or NA_real_ if undetermined.

Extract method used from assignment result

Description

Extract method used from assignment result

Usage

get_method_used(x)

Arguments

Value

Character string indicating method used

Extract total cost from assignment result

Description

Extract total cost from assignment result

Usage

get_total_cost(x)

Arguments

Value

Numeric total cost

Greedy match blocks in parallel

Description

Greedy match blocks in parallel

Usage

greedy_blocks_parallel(
  blocks,
  left,
  right,
  left_ids,
  right_ids,
  block_col,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  strategy,
  parallel = FALSE
)

Arguments

Value

List with combined results from all blocks

Re-export of dplyr::group_by

Description

Re-export of dplyr::group_by

Value

See group_by.

Check if data frame has blocking information

Description

Check if data frame has blocking information

Usage

has_blocks(df)

Value

Logical indicating whether data frame has block ID column.

Check if any valid pairs exist

Description

Check if any valid pairs exist

Usage

has_valid_pairs(cost_matrix)

Value

Logical indicating whether any valid pairs exist.

Hospital staff scheduling example dataset

Description

A comprehensive example dataset for demonstrating couplr functionality across vignettes. Contains 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
)

Low match rate info

Description

Low match rate info

Usage

info_low_match_rate(n_matched, n_left, pct)

Value

No return value, called for side effects (issues a message or warning).

Check if Object is a Distance Object

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

Check if object is a batch assignment result

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

Check if object is a k-best 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

Check if object is an 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

Join Matched Pairs with Original Data

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
)

Animate an assignment algorithm step-by-step

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)

Solve linear assignment problems

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)

Solve multiple assignment problems efficiently

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)

Find k-best optimal assignments

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)

Solve 1-D Line Assignment Problem

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)

Calipers of a lazy cost spec, keyed by variable name

Description

The C++ lazy cost source takes its calipers as a named list of thresholds, while the spec stores them as records carrying an index into spec$vars.

Usage

lazy_cost_spec_calipers(spec)

Value

Named list of numeric thresholds, one per caliper.

Precompute the Mahalanobis inverse covariance matrix for a lazy cost spec

Description

Mirrors compute_distance_matrix()'s pooled within-group covariance logic exactly (R/matching_distance.R) – computed once in R rather than reimplemented in C++, so the two code paths can't drift apart.

Usage

lazy_cost_spec_inv_cov(spec)

Value

p x p inverse covariance matrix, or NULL if distance != "mahalanobis".

The memory mode a lazy cost specification was built for

Description

A spec built before mode existed, or by hand, is a lazy one: solving every pair is what the class has always meant.

Usage

lazy_cost_spec_mode(spec)

Value

"lazy" or "implicit".

Compute paired (not cross) distances for specific matched pairs

Description

Given matched row/column index pairs (as produced by a solve), recomputes each pair's distance directly from left_mat/right_mat. This is cheap regardless of n_left/n_right: the number of matched pairs never exceeds min(n_left, n_right), so this never approaches the O(n*m) cost the lazy path exists to avoid. Mirrors compute_distance_matrix()'s per-metric formulas exactly, but pairwise rather than all-pairs.

Usage

lazy_pair_distances(spec, matched_rows, matched_cols)

Value

Numeric vector of length length(matched_rows).

Build a cost matrix from source / target / cost columns

Description

The long-format door of lap_solve() and lap_solve_kbest(). Cells no row names are forbidden, which is what makes an absent pair an absent edge rather than a zero-cost one.

Usage

long_to_cost_matrix(source_vals, target_vals, cost_vals, forbidden = NA)

Value

List with the matrix and the source / target level vectors its row and column indices stand for.

Mark forbidden pairs

Description

Generic function to mark specific pairs as forbidden.

Usage

mark_forbidden_pairs(cost_matrix, forbidden_indices)

Value

Modified cost matrix with forbidden pairs marked.

Apply the forbidden sentinel to a cost matrix

Description

Every front door that documents a forbidden argument masks with this one, so lap_solve(), lap_solve_batch() and lap_solve_kbest() read the same sentinel the same way. NA and Inf cells are forbidden to the solvers already, so forbidden = NA is the identity.

Usage

mask_forbidden(cost_matrix, forbidden = NA)

Value

The cost matrix with sentinel cells replaced by Inf.

Match blocks in parallel

Description

Match blocks in parallel

Usage

match_blocks_parallel(
  blocks,
  left,
  right,
  left_ids,
  right_ids,
  block_col,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  method,
  parallel = FALSE
)

Arguments

Value

List with combined results from all blocks

Match two datasets into couples

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.

Under memory_mode = "implicit" it also carries certificate, the checked proof of optimality (see verify_assignment()), 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")

Match with blocking (multiple problems)

Description

Match with blocking (multiple problems)

Usage

match_couples_blocked(
  left,
  right,
  left_ids,
  right_ids,
  block_col,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  method,
  strategy = "row_best",
  parallel = FALSE,
  replace = FALSE,
  ratio = 1L,
  sigma = NULL,
  memory_mode = "auto"
)

Value

List with pairs tibble and matching info.

Match from Precomputed Distance Object

Description

Internal function to handle matching when a distance_object is provided

Usage

match_couples_from_distance(
  dist_obj,
  max_distance = Inf,
  calipers = NULL,
  ignore_blocks = FALSE,
  require_full_matching = FALSE,
  method = "auto",
  strategy = "row_best",
  return_unmatched = TRUE,
  return_diagnostics = FALSE,
  check_costs = TRUE
)

Value

A matching_result object with pairs, info, and optional diagnostics.

Match without blocking (single problem)

Description

Match without blocking (single problem)

Usage

match_couples_single(
  left,
  right,
  left_ids,
  right_ids,
  vars,
  distance,
  weights,
  scale,
  max_distance,
  calipers,
  method,
  strategy = "row_best",
  check_costs = TRUE,
  replace = FALSE,
  ratio = 1L,
  sigma = NULL,
  memory_mode = "auto",
  certify = NULL
)

Value

List with pairs tibble and matching info.

Extract Analysis-Ready Data from Matching Results

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)

Match across a range of one design choice

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

Create blocks for stratified matching

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)

Construct a lazy cost specification

Description

mode is the memory mode that resolved to this specification: "lazy", solved over every pair, or "implicit", solved by generating the pairs the answer turns out to need.

Usage

new_lazy_cost_spec(
  left_mat,
  right_mat,
  distance,
  sigma,
  weights,
  vars,
  mode = c("lazy", "implicit")
)

Value

An object of class "lazy_cost_spec".

Parallel lapply using future

Description

Parallel lapply using future

Usage

parallel_lapply(X, FUN, ..., parallel = FALSE)

Arguments

Value

List of results

Pixel-level image morphing (final frame only)

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)
  }
}

Pixel-level image morphing (animation)

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)
  }
}

Plot method for balance diagnostics

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)

Plot method for matching results

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)

Plot method for sensitivity analysis (base graphics)

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)

Preprocess matching variables with automatic checks and scaling

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

Print Method for Balance Diagnostics

Description

Print Method for Balance Diagnostics

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print Method for CEM Results

Description

Print Method for CEM Results

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print Method for Flow Problems

Description

Print Method for Flow Problems

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print a design path

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.

Print Method for Distance Objects

Description

Print Method for Distance Objects

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print Method for Full Matching Results

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.

Print method for batch assignment results

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.

Print method for k-best assignment results

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.

Print method for assignment results

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.

Print method for matching results

Description

Print method for matching results

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print method for matchmaker results

Description

Print method for matchmaker results

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print method for preprocessing result

Description

Print method for preprocessing result

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print method for sensitivity analysis

Description

Print method for sensitivity analysis

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print Method for Subclassification Results

Description

Print Method for Subclassification Results

Usage

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

Arguments

Value

Invisibly returns the input object x.

Print method for variable health

Description

Print method for variable health

Usage

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

Arguments

Value

Invisibly returns the input object x.

Propensity Score Matching

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)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

augment()

Resolve a requested memory_mode to a concrete decision

Description

Resolve a requested memory_mode to a concrete decision

Usage

resolve_memory_mode(
  n,
  m,
  memory_mode = c("auto", "dense", "lazy", "implicit"),
  solver_supports_lazy = FALSE,
  solver_supports_implicit = FALSE,
  ram_fraction = 0.5,
  fallback_threshold_mb = 4000
)

Arguments

Value

"dense", "lazy" or "implicit".

Restore original parallel plan

Description

Restore original parallel plan

Usage

restore_parallel(parallel_state)

Arguments

Value

No return value, called for side effects (restores parallel plan).

Rosenbaum Sensitivity Analysis

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

Details

Rosenbaum (2002, Chapter 4) bounds quantify how much hidden bias (an unobserved confounder) would be needed to explain away the observed treatment effect. The sensitivity parameter Gamma represents the maximum ratio of treatment odds between two matched units:

• Gamma = 1: No hidden bias (standard Wilcoxon test)

• Gamma = 2: One unit could be twice as likely to receive treatment due to an unobserved factor

The function computes upper and lower bounds on the p-value of the Wilcoxon signed-rank test under each level of hidden bias. A finding is "insensitive to bias" if p_upper remains below 0.05 even at large Gamma.

Value

An S3 object of class sensitivity_analysis containing:

results

Tibble with columns: gamma, t_stat, p_upper, p_lower

n_pairs

Number of matched pairs analyzed

critical_gamma

Smallest gamma at which p_upper > 0.05

alternative

Direction of test

References

Rosenbaum, P.R. (2002). Observational Studies, 2nd edition. Springer.

Examples

set.seed(42)
left <- data.frame(id = 1:20, x = rnorm(20), outcome = rnorm(20, 1, 1))
right <- data.frame(id = 21:40, x = rnorm(20), outcome = rnorm(20, 0, 1))
result <- match_couples(left, right, vars = "x")
sens <- sensitivity_analysis(result, left, right, outcome_var = "outcome")
print(sens)

Setup parallel processing with future

Description

Setup parallel processing with future

Usage

setup_parallel(parallel = FALSE, n_workers = NULL)

Arguments

Value

List with original plan and whether we set up parallelization

'Sinkhorn-Knopp' optimal transport solver

Description

Compute an entropy-regularized optimal transport plan using the 'Sinkhorn-Knopp' algorithm. Unlike other LAP solvers that return a hard 1-to-1 assignment, this returns a soft assignment (doubly stochastic matrix).

Usage

sinkhorn(
  cost,
  lambda = 10,
  tol = 1e-09,
  max_iter = 1000,
  r_weights = NULL,
  c_weights = NULL
)

Arguments

Details

The 'Sinkhorn-Knopp' algorithm solves the entropy-regularized optimal transport problem:

P^* = \arg\min_P \langle C, P \rangle - \frac{1}{\lambda} H(P)

subject to row sums = r_weights and column sums = c_weights.

The entropy term H(P) encourages spread in the transport plan. As lambda -> Inf, the solution approaches the standard (unregularized) optimal transport.

Key differences from standard LAP solvers:

• Returns a soft assignment (probabilities) not a hard 1-to-1 matching

• Supports unequal marginals (weighted distributions)

• Differentiable, making it useful in ML pipelines

• Very fast: O(n^2) per iteration with typically O(1/tol^2) iterations

Use sinkhorn_to_assignment() to round the soft assignment to a hard matching.

Value

A list with elements:

• transport_plan — numeric matrix, the optimal transport plan P. Row sums approximate r_weights, column sums approximate c_weights.

• cost — the transport cost <C, P> (without entropy term).

• u, v — scaling vectors (P = diag(u) * K * diag(v) where K = exp(-lambda*C)).

• converged — logical, whether the algorithm converged.

• iterations — number of iterations used.

• lambda — the regularization parameter used.

References

Cuturi, M. (2013). 'Sinkhorn Distances': Lightspeed Computation of Optimal Transport. Advances in Neural Information Processing Systems, 26.

See Also

assignment() for hard 1-to-1 matching, sinkhorn_to_assignment() to round soft assignments.

Examples

cost <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)

# Soft assignment with default parameters
result <- sinkhorn(cost)
print(round(result$transport_plan, 3))

# Sharper assignment (higher lambda)
result_sharp <- sinkhorn(cost, lambda = 50)
print(round(result_sharp$transport_plan, 3))

# With custom marginals (more mass from row 1)
result_weighted <- sinkhorn(cost, r_weights = c(0.5, 0.25, 0.25))
print(round(result_weighted$transport_plan, 3))

# Round to hard assignment
hard_match <- sinkhorn_to_assignment(result)
print(hard_match)

Round 'Sinkhorn' transport plan to hard assignment

Description

Convert a soft transport plan from sinkhorn() to a hard 1-to-1 assignment using greedy rounding.

Usage

sinkhorn_to_assignment(result)

Arguments

Details

Greedy rounding iteratively assigns each row to its most probable column, ensuring no column is assigned twice. This may not give the globally optimal hard assignment; for that, use the transport plan as a cost matrix with assignment().

Value

Integer vector of column assignments (1-based), same format as assignment().

See Also

sinkhorn()

Examples

cost <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)
result <- sinkhorn(cost, lambda = 20)
hard_match <- sinkhorn_to_assignment(result)
print(hard_match)

Solver status values

Description

The closed set of values a couplr solve result's status field can take.

Usage

solver_status_values()

Details

"optimal"

The solver terminated on its own optimality condition and every row it was asked to match is matched.

"partial"

Feasible, but fewer rows are matched than requested, under a maximum-cardinality-then-minimum-cost objective.

"infeasible"

No assignment of the requested cardinality exists.

"eps_optimal"

An auction solve terminated at a caller-supplied epsilon. The result is within n * epsilon of optimal, which is a bound, not a proof of optimality.

"iteration_limit"

The solver stopped on an iteration cap rather than on optimality. The solution is feasible; its optimality is unproven.

"interrupted"

The solver stopped on a time budget or a user interrupt, part way through. What it had placed respects every capacity, and is short of what the problem asked it to place, so it is neither an answer nor evidence that no answer exists.

"heuristic"

Produced by a heuristic. Optimality is neither claimed nor checked.

Value

Character vector of the permitted status values.

See Also

verify_assignment() to turn a status into a checked certificate.

Examples

solver_status_values()

Calculate Standardized Difference

Description

Computes the standardized mean difference between two groups. This is a key metric for assessing balance in matched samples.

Usage

standardized_difference(x1, x2, pooled = TRUE, w1 = NULL, w2 = NULL)

Arguments

Details

Standardized difference = (mean1 - mean2) / pooled_sd where pooled_sd = sqrt((sd1^2 + sd2^2) / 2)

Common thresholds: less than 0.1 is excellent balance, 0.1-0.25 is good balance, 0.25-0.5 is acceptable balance, and greater than 0.5 is poor balance.

Value

Numeric value representing the standardized difference

Subclassification on Propensity Score

Description

Divides units into K strata based on quantiles of the propensity score, then computes within-stratum weights for treatment effect estimation. This is a simple, transparent approach to propensity score adjustment that allows visual inspection of balance within each subclass.

Usage

subclass_match(
  formula = NULL,
  data = NULL,
  treatment = NULL,
  n_subclasses = 5L,
  ps = NULL,
  ps_model = NULL,
  estimand = "ATT"
)

Arguments

Details

The algorithm:

1. Estimate propensity scores via logistic regression (or use pre-computed scores)

2. Divide the propensity score distribution into K quantile-based strata

3. For each stratum, check overlap (both treated and control units present)

4. Compute within-stratum weights based on the target estimand:

   • ATT: Treated units get weight 1; control units get weight n_treated_in_stratum / n_control_in_stratum

   • ATE: Both groups get weight proportional to stratum size relative to total sample

   • ATC: Control units get weight 1; treated units get weight n_control_in_stratum / n_treated_in_stratum

Value

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

matched

Tibble with columns id, side, subclass, ps, weight

subclass_summary

Tibble with per-subclass statistics: counts, mean PS, and overlap status

info

List with n_subclasses, estimand, n_left, n_right, method, vars

Examples

set.seed(42)
n <- 200
data <- data.frame(
  id = 1:n,
  age = rnorm(n, 40, 10),
  income = rnorm(n, 50000, 15000)
)
data$treatment <- rbinom(n, 1, plogis(-2 + 0.05 * data$age))
result <- subclass_match(treatment ~ age + income, data, treatment = "treatment")
print(result)

Perfect balance success message

Description

Perfect balance success message

Usage

success_good_balance(mean_std_diff)

Value

No return value, called for side effects (issues a message).

Suggest scaling method based on variable characteristics

Description

Analyzes variable distributions and suggests appropriate scaling methods.

Usage

suggest_scaling(left, right, vars)

Arguments

Value

A character string with the suggested scaling method: "standardize", "range", "robust", or "none"

Summarize block structure

Description

Summarize block structure

Usage

summarize_blocks(left, right, block_vars = NULL)

Value

Tibble with block_id, n_left, n_right, and optional variable means.

Summary method for balance diagnostics

Description

Summary method for balance diagnostics

Usage

## S3 method for class 'balance_diagnostics'
summary(object, ...)

Arguments

Value

A list containing summary statistics (invisibly)

Summary Method for Distance Objects

Description

Summary Method for Distance Objects

Usage

## S3 method for class 'distance_object'
summary(object, ...)

Arguments

Value

Invisibly returns the input object.

Get summary of k-best results

Description

Extract summary information from k-best assignment results.

Usage

## S3 method for class 'lap_solve_kbest_result'
summary(object, ...)

Arguments

Value

A tibble with one row per solution containing:

• rank: solution rank

• solution_id: solution identifier

• total_cost: total cost of the solution

• n_assignments: number of assignments in the solution

Summary method for matching results

Description

Summary method for matching results

Usage

## S3 method for class 'matching_result'
summary(object, ...)

Arguments

Value

A list containing summary statistics (invisibly)

Summary method for sensitivity analysis

Description

Summary method for sensitivity analysis

Usage

## S3 method for class 'sensitivity_analysis'
summary(object, ...)

Arguments

Value

A list containing summary statistics (invisibly)

Swap left/right in a lazy cost spec

Description

A cheap metadata field-swap (left_mat <-> right_mat, n_left <-> n_right), unlike the dense path's t() matrix copy. Calipers/max_distance are unaffected: a caliper's var_index refers to a matching VARIABLE (a column shared by both sides), not a left/right unit index, so it does not need to change when the roles of left/right are swapped.

Usage

transpose_lazy_cost_spec(spec)

Update Constraints on Distance Object

Description

Apply new constraints to a precomputed distance object without recomputing the underlying distances. This is useful for exploring different constraint scenarios quickly.

Usage

update_constraints(dist_obj, max_distance = Inf, calipers = NULL)

Arguments

Details

This function creates a new distance_object with modified constraints applied to the cost matrix. The original distance_object is not modified.

Constraints:

• max_distance: Sets cost to Inf for pairs exceeding this threshold

• calipers: Per-variable restrictions (e.g., calipers = list(age = 5))

The function returns a new object rather than modifying in place, following R's copy-on-modify semantics.

Value

A new distance_object with updated cost_matrix

Examples

left <- data.frame(id = 1:5, age = c(25, 30, 35, 40, 45))
right <- data.frame(id = 6:10, age = c(24, 29, 36, 41, 44))
dist_obj <- compute_distances(left, right, vars = "age")

# Apply constraints
constrained <- update_constraints(dist_obj, max_distance = 2)
result <- match_couples(constrained)

Check if emoji should be used

Description

Check if emoji should be used

Usage

use_emoji()

Value

Logical indicating whether emoji should be used.

Validate calipers parameter

Description

Validate calipers parameter

Usage

validate_calipers(calipers, vars)

Value

Validated calipers (list or named numeric), or NULL if none.

Validate and prepare cost data

Description

Internal helper that ensures a numeric, non-empty cost matrix.

Usage

validate_cost_data(x, forbidden = NA)

Arguments

Value

Numeric cost matrix

Validate matching inputs

Description

Validate matching inputs

Usage

validate_matching_inputs(left, right, vars = NULL)

Value

Invisibly returns TRUE if validation passes; otherwise throws an error.

Validate weights parameter

Description

Validate weights parameter

Usage

validate_weights(weights, vars)

Value

Numeric vector of validated weights.

Verify that an assignment is optimal

Description

Checks a solved assignment against the linear-programming optimality conditions and returns the result of each check. Unlike the status field on a solve result, which records what the solver terminated on, this is a proof: certified_optimal is TRUE only when every condition holds.

Usage

verify_assignment(x, cost = NULL, duals = NULL, maximize = FALSE, tol = 1e-09)

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

Arguments

Details

The check needs dual variables. If x carries them (as assignment_duals() results do), they are used. Otherwise they are obtained by solving cost with assignment_duals(), which costs a second solve. Either way the duals are verified, not trusted: dual feasibility is checked over every admissible pair, so duals that do not certify anything cause the verification to fail rather than pass.

Optimal duals are shared by all optimal solutions of a linear program, so a matching from one solver can be certified against duals from another. That is what makes it possible to certify solvers that return no duals of their own.

Value

An object of class assignment_certificate, a list with elements:

• certified_optimal — logical, the conclusion. TRUE only when every condition below holds.

• primal_feasible — logical; no column claimed twice, no forbidden pair matched, no index out of range.

• dual_feasible — logical; c_ij - u_i - v_j >= -tol over every admissible pair, and, when there are more columns than rows, v_j <= tol for every column.

• complementary_slackness — logical; both halves.

• cs_matched_tight, cs_unmatched_free — the two halves separately.

• primal_objective, dual_objective, duality_gap — numeric.

• min_reduced_cost, worst_i, worst_j — the most violated pair, if any.

• max_matched_slack, max_v_unmatched, max_v — the quantities the slackness and sign conditions bound.

• n_matched, n_rows, n_cols, transposed, tolerance.

Invisibly returns x.

See Also

assignment(), assignment_duals(), solver_status_values()

Examples

set.seed(1)
cost <- matrix(runif(100), 10, 10)
verify_assignment(assignment(cost), cost)

# A rectangular problem, where the condition on unmatched columns bites.
# Passing the duals result reuses its duals instead of solving again.
rect <- matrix(runif(120), 6, 20)
verify_assignment(assignment_duals(rect), rect)

Verify that a flow is optimal

Description

Checks a flow and a set of node potentials against the linear-programming optimality conditions for the minimum-cost flow problem, and returns the result of each check. Unlike the status field on a solve result, which records what the solver terminated on, this is a proof: certified_optimal is TRUE only when every condition holds.

Usage

verify_flow(x, problem = NULL, potential = NULL, tol = 1e-09)

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

Arguments

Details

The problem is

  min  sum_a cost(a) f(a)
  s.t. (flow out of v) - (flow into v) = supply(v)   at every node v
       lower(a) <= f(a) <= upper(a)                  on every arc a

and, writing cbar(a) = cost(a) + pi(tail(a)) - pi(head(a)) for the reduced cost under the potentials pi, the conditions checked are primal feasibility, cbar(a) >= -tol(a) on every arc that can still take flow, and cbar(a) <= tol(a) on every arc carrying more than its lower bound. Their objective form is checked too: the duality gap is the sum of the slackness violations weighted by ⁠|cbar|⁠, so it is where a violation too small to trip the per-arc tolerance still accumulates.

tol is relative. Each arc is compared against

  tol(a) = tol * max(1, |cost(a)|, |pi(tail(a))|, |pi(head(a))|)

because cbar(a) is computed from those three numbers, and its last bits are worth the largest of them times the machine epsilon. A design that stacks lexicographic tier weights, as cardinality_match() does to rank cardinality above balance above distance, reaches potentials in the millions, where one unit in the last place is around 1e-9 and an exactly optimal flow cannot meet an absolute 1e-9. The scale never falls below 1, so a problem whose costs and potentials are of order 1 is checked against tol itself. The widest tolerance any comparison used is reported as dual_tolerance.

The check needs potentials. If x carries them, as a solve result does, they are used. Otherwise they are obtained by solving problem, which costs a second solve. Either way they are verified, not trusted: an arc priced below the tolerance makes the verification fail rather than pass.

Optimal potentials are shared by all optimal solutions of a linear program, so a flow from one solver can be certified against potentials from another.

Value

An object of class flow_certificate, a list with elements:

• certified_optimal - logical, the conclusion. TRUE only when every condition below holds.

• primal_feasible - logical; every arc inside its bounds and every node's net flow equal to its supply.

• n_capacity_violations, n_conservation_violations, max_conservation_error - what primal feasibility failed on.

• dual_feasible - logical; no arc that can still take flow prices below its own -tol(a).

• complementary_slackness - logical; no arc above its lower bound prices above its own tol(a).

• n_cs_violations, min_residual_reduced_cost, worst_arc - the smallest reduced cost over the residual graph and the arc attaining it, which is a violation when it falls below that arc's -tol(a). worst_arc is 0 when no arc can either take or give up flow.

• dual_tolerance - the widest tol(a) any comparison was made against, so the verdict names the resolution it was reached at.

• primal_objective, dual_objective, duality_gap - numeric.

• tolerance - the relative tol as supplied.

Invisibly returns x.

See Also

verify_assignment(), solver_status_values()

Examples

# Two supply nodes shipping to two demand nodes, stated directly as a flow.
prob <- list(
  n_nodes = 4,
  supply  = c(2, 1, -2, -1),
  arcs = data.frame(
    tail  = c(1, 1, 2, 2),
    head  = c(3, 4, 3, 4),
    lower = c(0, 0, 0, 0),
    upper = c(2, 2, 2, 2),
    cost  = c(1, 3, 2, 1)
  )
)

# Both of node 1's units go to node 3, node 2's unit goes to node 4.
verify_flow(c(2, 0, 0, 1), prob, potential = c(0, 0, 1, 1))

# The same flow against potentials that certify nothing.
verify_flow(c(2, 0, 0, 1), prob, potential = c(0, 0, 0, 0))

Available megabytes implied by a block of vm_stat output

Description

Inactive and speculative pages are reclaimed on demand, so they count as available to an allocation in the same sense as Linux's MemAvailable; macOS keeps almost nothing on the free list, so the free count alone understates what a large matrix can actually obtain.

Usage

vm_stat_available_mb(vm, page_size)

Arguments

Details

Split from get_free_ram_mb() so the page-size and page-class handling can be checked without a macOS host.

Value

Numeric scalar (MB available), or NA_real_ if unparseable.

Page size behind a block of vm_stat output

Description

vm_stat counts in pages and states its own page size in the header line ("page size of 16384 bytes"). Apple Silicon pages are 16K and Intel pages are 4K, so the size is read rather than assumed; ⁠sysctl hw.pagesize⁠ is the fallback when the header cannot be parsed.

Usage

vm_stat_page_size(vm)

Arguments

Value

Page size in bytes, or NA_real_ if undetermined.

All distances identical warning

Description

All distances identical warning

Usage

warn_constant_distance(value)

Value

No return value, called for side effects (issues a warning).

Constant variable warning

Description

Constant variable warning

Usage

warn_constant_var(var)

Value

No return value, called for side effects (issues a warning).

Extreme cost ratio warning

Description

Extreme cost ratio warning

Usage

warn_extreme_costs(p95, p99, ratio, problem_vars = NULL)

Value

No return value, called for side effects (issues a warning).

Many forbidden pairs warning

Description

Many forbidden pairs warning

Usage

warn_many_forbidden(pct_forbidden, n_valid, n_left)

Value

No return value, called for side effects (issues a warning).

Too many zeros warning

Description

Too many zeros warning

Usage

warn_many_zeros(pct, n_zeros)

Value

No return value, called for side effects (issues a warning).

Parallel package missing warning (reuse from matching_parallel.R)

Description

Parallel package missing warning (reuse from matching_parallel.R)

Usage

warn_parallel_unavailable()

Value

No return value, called for side effects (issues a warning).

High distance matches warning

Description

High distance matches warning

Usage

warn_poor_quality(pct_poor, threshold)

Value

No return value, called for side effects (issues a warning).