greedy_couples() is removed; greedy matching is
now match_couples(method = "greedy"). The two
functions duplicated ~130 lines of identical scaffolding (validation,
scaling, id extraction, blocking dispatch, metadata) over the same
shared engine. They are now one front door: match_couples()
gains a method = "greedy" value and a
strategy argument (“row_best”, “sorted”, “pq”). Replace
greedy_couples(x, strategy = "sorted") with
match_couples(x, method = "greedy", strategy = "sorted").
The result object and info$method == "greedy" are
unchanged.pixel_morph() and pixel_morph_animate()
gain a mode = "color_match" option: pixels sharing a
quantized colour are matched spatially and any remainder falls back to
identity. A lighter-weight alternative to the default
"color_walk" palette LAP.sensitivity_analysis() no longer scrambles
matched pairs. Outcomes were assembled with two independent
merge() calls, each sorted by its own key, so the pair
difference subtracted outcomes from mismatched pairs and every
downstream quantity (Wilcoxon T+, Rosenbaum bounds, critical gamma) was
computed on a scrambled pairing. Outcomes are now looked up by ID,
preserving the row-wise pair correspondence (#4).
subclass_match(estimand = "ATE") weights
corrected. ATE subclass weights carried an extra factor of the
stratum size, over-weighting large subclasses quadratically. A treated
unit in subclass k now carries (n_k / N) / n_t as intended;
ATT and ATC were already correct (#5).
balance_diagnostics() now applies stratum
weights for full matching, CEM, and subclassification. The
weights were computed and discarded, so standardized differences were
unweighted for the very estimators whose balance is achieved through
weighting. Weighted mean, variance, and standardized difference are now
used on both sides. Also: the variance ratio is now a true ratio of
variances (matching the conventional 0.5-2 bounds) rather than a ratio
of standard deviations, and the unmatched right-unit count no longer
goes negative under ratio > 1 / replace (it
counts distinct matched right units, not pair rows) (#6).
ssap_bucket no longer silently rounds
fractional costs to a wrong optimum. The integer-scaling step
tried only multipliers {1, 10, 100, 1000} and, failing
those, rounded at 1000 – flipping which permutation was
optimal on costs needing more than three decimals. It now searches
ascending powers of ten with a fixed (scale-independent) integrality
tolerance and refuses the problem, redirecting to
method = "jv" or "auction", when no bounded
integer scaling is exact. The animation mirror
trace_ssap_bucket() applies the same rule (#19).
lap_solve_line_metric(maximize = TRUE) now
returns the true maximum-weight matching. The DP always built
the sorted (minimum-cost) pairing and merely negated the total; on a
line the maximum-weight matching is the anti-monotone pairing. The DP
now runs against the descending target ordering and returns that
assignment and its true total (#8).
gabow_tarjan returns a perfect matching when
the diagonal is forbidden. The C_max == 0 fast
path assigned the diagonal without checking feasibility, returning an
empty matching when the diagonal cells were forbidden but a perfect
matching existed. It now finds a maximum-cardinality matching over the
allowed edges via augmenting paths (#9).
compute_distances(auto_scale = TRUE) now
scales. It read a nonexistent field (making vars
NULL) and disabled scaling under the belief it had already
happened. It now reads the selected variables and forwards the chosen
scaling method to the cost builder (#7).
Pre-fitted propensity models predict on the supplied
data. ps_match() and subclass_match()
called predict() without newdata =, so a
ps_model fitted on a differently ordered or subset frame
attached scores to the wrong rows. They now pass
newdata = data (#10).
match_couples(ratio > 1) falls back to a
partial match on infeasibility. The ratio > 1
path called the solver directly and hard errored when constraints
forbade every edge of some unit, whereas the 1:1 path returned a partial
matching. Both paths now share the same partial-feasibility / greedy
fallback (#10).
lap_solve() honors the
forbidden sentinel for matrix input, and
lap_solve_batch() preserves singleton-dimension
orientation. The matrix path silently ignored a
non-NA forbidden; it now masks matching cells
as forbidden. A 3-D array slice with a singleton row or column dimension
was dropped to a vector and transposed; slices are now reshaped
explicitly (#11).
Guarded the C++ solvers against silently wrong results and crashes at extreme scale. None of these affect ordinary inputs; they add error paths and 64-bit arithmetic where 32-bit overflow or a fixed tolerance could previously produce a wrong “optimal” or a crash (#13):
Overflow / narrowing. ssap_bucket
errors clearly when cost magnitudes exceed what the integer-bucket
solver can represent (rather than overflowing the sentinel or allocating
an enormous bucket queue); the network-simplex iteration bound and
gabow_tarjan’s bit-scaling range are computed and checked
in 64-bit; gabow_tarjan also rejects costs that collide
with its forbidden sentinel; the brute-force solver caps total
enumeration work instead of running unbounded in the number of
columns.
Large n*m indexing. Flat
cost/kernel indexing in prepare_cost_matrix,
solve_sinkhorn, and the auction epsilon is done in 64-bit;
network_simplex and lapmod reject problems
whose arc / entry counts would overflow a 32-bit index.
Tolerances and status.
solve_munkres scales its zero tolerance with the cost
magnitude (a fixed 1e-12 could make a solvable large-cost
matrix throw); solve_csa scales non-integer costs to
integers before the epsilon-scaling auction, so its optimality guarantee
(which assumes integer costs) also holds for real-valued inputs with
near-tied assignments; full_matching now reports
infeasible when the group capacity is below the number of
units instead of silently dropping units as optimal;
solve_sinkhorn reports the correct iteration count on
non-convergence.
hk01 fallback. The pure solve_hk01
now falls back to the exact weighted solver (solve_csflow)
when the zero-cost subgraph of a {0,1} matrix has no
perfect matching, instead of erroring – matching the Rcpp path that
assignment(method = "hk01") already used.
Bounds. The internal
morph_pixel_level helpers assert their pixel / assignment
buffer sizes, matching the exported wrappers.
csa shipped path now carries the fixes it
was tested for. The Rcpp entry point for
method = "csa" ran a separate copy of the solver that never
received the integer-scaling fix above, so
assignment(method = "csa") could still return a suboptimal
matching on fractional costs. It now delegates to the single pure
solve_csa implementation exercised by the C++ tests. That
implementation also gained square padding for rectangular problems,
which it previously solved greedily (and suboptimally). The Rcpp
*_impl wrappers now share one
rcpp_to_cost_matrix / lap_result_to_rcpp
conversion pair instead of per-file copies (#15).
Remaining solvers now ship the tested
implementation. Following csa, the Rcpp entry
points for sap/ssp, csflow,
cycle_cancel, push_relabel,
ssap_bucket, bruteforce,
bottleneck, hk01,
network_simplex, and the three auction
variants each ran a second copy of the algorithm that had drifted from
the pure lap::solve_* exercised by the C++ tests. They now
delegate to that single pure implementation, so the shipped path and the
tested path are identical. Each pure copy was checked against brute
force over randomised integer, fractional, rectangular, maximize, and
forbidden-edge inputs before its wrapper was pointed at it.
network_simplex thereby picks up the pure copy’s
O(n^2) pivot bound (the shipped copy used the slower
O(arcs * nodes) bound).
auction, auction_gs, and
auction_scaled now return the exact optimum. The
basic and Gauss-Seidel auctions used a single fixed epsilon, which
leaves a duality-gap slack of up to n * eps and returned
suboptimal matchings on closely-spaced costs
(assignment(method = "auction") could disagree with
jv); auction_scaled additionally threw on some
feasible rectangular problems with forbidden edges under
maximize. All three now run one shared epsilon-scaling core
that scales epsilon down to a tiny final value, recovering the exact
assignment. lap_solve_auction_gs() keeps its
bids diagnostic.
hk01 maximize. The pure
solve_hk01 flipped maximize by negation,
turning a {0,1} matrix into {0,-1}, which its
palette check no longer recognised as binary – so it threw on feasible
binary maximization problems. It now flips via cmax - c,
preserving the {0,1} palette so the fast path and the
solve_csflow fallback engage.
Greedy matching wrappers
(greedy_matching, _sorted,
_row_best, _pq) now delegate to the pure
lap::greedy_matching_*. To keep the shipped behaviour
identical, the pure copies gained the two tolerances the Rcpp copies had
and they lacked: they skip the large-finite BIG sentinel
the matching layer uses for forbidden edges (so a row whose only
remaining options are forbidden is left unmatched rather than paired to
a forbidden column), and they accept n > m by returning
a partial matching instead of erroring. Verified byte-identical to the
previous wrappers over 400 randomised cases spanning integer ties,
NA/BIG-forbidden edges, and rectangular
shapes. The three per-strategy Rcpp exports
(greedy_matching_sorted / _row_best /
_pq) were folded into the single
greedy_matching(strategy = ...) dispatcher they duplicated;
greedy_couples(strategy = ...) remains the user-facing
verb.
test-statistical-recovery.R): sensitivity pair
alignment, prefitted-PS row alignment, propensity-matching imbalance
reduction, ATE subclass weight values, weighted-balance means,
known-effect recovery across seeds, and nominal coverage of matched-pair
confidence intervals (#14).cpp_tests/tests/test_ground_truth.cpp) that compares every
optimal pure solver against brute-force enumeration over thousands of
integer, fractional, rectangular, maximize, and forbidden-edge matrices
(bottleneck against a brute-force minimax). This is the gate that
decides whether a solver’s Rcpp wrapper may delegate to the pure copy,
and it is what surfaced the auction and hk01
bugs above.Fixes two solver paths that could stall indefinitely on
match_couples() inputs with max_distance,
calipers, or other forbidden-edge constraints. These stalls caused the
M1mac and linux-arm64 additional CRAN checks for 1.4.0 to hit the
1.5-hour test timeout.
Forbidden-cell marker is now Inf instead of
a large finite value. apply_max_distance(),
apply_calipers(), and mark_forbidden_pairs()
previously wrote a large finite BIG_COST into forbidden
cells. The Jonker-Volgenant and small-n SSP solvers treated
BIG_COST as a regular expensive edge and could degenerate
on sparse, near-square inputs instead of short-circuiting on
infeasibility. Switched to Inf so the C++ solvers’
non-finite check fires.
Auto-dispatch no longer routes sparse inputs through SSP
for small n. Previously lap_solve()
with method = "auto" selected "sap"
(lap_solve_ssp) for sparse matrices with
n <= 100. SSP has its own worst-case stall on
near-square, highly-sparse cost matrices. All sparse inputs now go
through lapmod regardless of size.
match_couples() now drops fully-forbidden
rows/columns before LAP. match_couples() and
.couples_from_distance() route through a new internal
.solve_with_partial_feasibility() helper. It removes rows
and columns with no allowed edges before the LAP call and falls back to
greedy_matching() if the optimal solver still cannot find a
perfect matching on the feasibility-pruned submatrix (Hall’s-condition
violation). Dropped rows/columns are returned as unmatched, preserving
the partial-matching semantics that tests with tight
max_distance / caliper constraints already
expected.
jv_core: drop the same-pass reprocess
in AUGMENTING ROW REDUCTION. The reprocess could revisit a
freshly-reduced row in the same pass and delay convergence on degenerate
inputs without changing the final assignment.lap_animate() now covers every method that
assignment() accepts. Ten new step-by-step traces
ship: auction_gs, ramshaw_tarjan,
ssap_bucket, hk01, csflow,
cycle_cancel, push_relabel, csa,
orlin, network_simplex.
animated_methods() returns all 20 method strings.testthat suite
(tests/testthat/test-trace-parity.R) on a battery of small
cost matrices including forbidden cells. Each frame’s matching is
validated for in-range entries, no double-bookings, and no use of
forbidden edges; the final-frame total is compared to the C++ oracle
within tolerance.R/trace_helpers_frame.R (make_frame(),
make_meta(), prepare_cost_work(),
matching_total_cost(), validate_cost_input())
and R/trace_helpers_mcf.R (min-cost-flow graph, residual
edges, Dijkstra with Johnson potentials, Bellman-Ford, negative-cycle
finder, push/extract). Used by all min-cost-flow traces.prepare_cost_matrix.cpp: entries equal
to +Inf were treated as regular very-large costs rather
than forbidden, which made cmax become Inf and
silently skipped the maximize flip. Result:
assignment(method = X, maximize = TRUE) on matrices
containing Inf returned the minimizing answer for
any solver routing through prepare_cost_matrix_impl
(auction, auction_scaled, sap,
csflow, hk01, bruteforce). Now
NA and any non-finite value are marked forbidden
consistently.lap_solve_orlin and
lap_solve_network_simplex_wrapper: the R-side
wrapper used work[is.na(work)] <- Inf which missed the
-Inf produced by negating +Inf in maximize
mode, letting forbidden cells slip through as extreme-cost real edges.
Fixed to work[!is.finite(work)] <- Inf.network_simplex initial spanning tree:
the greedy initialiser in ns_init.h built a
partial matching (any row that couldn’t claim a fresh column
was left unmatched) and connected unmatched columns to row
maximize,
assignment(method = "network_simplex") returned an
infeasible result with one row unmatched. Fixed by adding an
augmenting-path repair after the greedy pass: every still-unmatched row
runs BFS for an augmenting path on the allowed-edge bipartite graph,
extending the initial matching to a perfect matching whenever one
exists.method = "hungarian" now uses the shortest-augmenting-path
solver shared with JV; the original O(n^4) Munkres implementation
remains available as method = "munkres". At n = 2000 the
new Hungarian runs orders of magnitude faster than 1.3.2.auction and
auction_gs. Cleaner inner loop; no behaviour
change.solve_auction_scaled collapsed into a thin
wrapper over scaled_params (~200 lines removed);
behaviour identical.r > bn pruning is the
algorithm, not a wart); added the 6n pruning heuristic from p.9.paper/benchmark-table.csv and
paper/scaling-results.csv re-measured on the current
development machine for n <= 2000 (per-method table) and n_total
<= 2000 (cross-package table). Larger-n rows in both files are
carried over from the previous machine and not directly comparable.test-lap-solve-batch-coverage.R. Debian r-devel, local
r-release, and local R CMD check --as-cran all pass; the
crash did not reproduce off win-builder.Config/testthat/parallel: true removed from DESCRIPTION)
to eliminate cross-file worker-state leakage as a possible cause of the
win-builder crash.skip_on_cran() at the top of
test-lap-solve-batch-coverage.R. Equivalent coverage is
exercised off-CRAN by test-lap-solve-batch-coverage-2.R,
test-lap-solve-batch-coverage-3.R,
test-lap-solve-batch-extended.R,
test-batch-coverage-final.R,
test-batch-processing.R, and
test-batch-kbest-extended.R.rbind(left, right). The pooled
within-group estimator
((n_L-1)*S_L + (n_R-1)*S_R) / (n_L+n_R-2) is the convention
used by optmatch::match_on() and aligns Mahalanobis
behaviour across the matching packages a user is likely to compare
against. Users who relied on the old default can recover it explicitly
with
match_couples(..., sigma = cov(rbind(left[, vars], right[, vars]))).
The previous docstring already documented the default as “pooled
covariance”; this release makes the code match the documentation.full_match() gains
method = "optimal" (new default) using a min-cost
max-flow solver (Dijkstra + Johnson potentials) that finds the globally
optimal group assignment minimizing total distance:
min_controls per groupn_left > n_rightsolve_full_matching.cpp (self-contained
MCMF)method = "greedy" preserved for fast approximate
matchingfull_match() examplefull_match() function assigns
every unit to a matched group with variable ratios (1:k or k:1):
caliper (absolute) or
caliper_sd (SD-based)min_controls,
max_controlsfull_matching_result S3 classcem_match() function implements
coarsened exact matching:
grouping
parametercutpoints
parametercem_result S3 class with matched units and
strata summarysubclass_match() function divides
units into propensity score strata:
subclass_result S3 class with subclass
summarymatch_data() generic converts any
couplr result to analysis-ready format with treatment,
weights, subclass, and distance
columns. Methods for all result types (matching, full, CEM,
subclass).as_matchit() converter creates
matchit-class objects from couplr results, enabling interop
with cobalt, marginaleffects, and other MatchIt ecosystem packages.bal.tab() methods for all
couplr result types. Requires cobalt package (in Suggests).rcond()
instead of fragile det() == 0sigma parameter in
match_couples(), greedy_couples(), and
compute_distance_matrix() for user-supplied covariance
matricesbalance_diagnostics() and join_matched()
are now S3 generics with methods for all result types. Existing code is
100% backward-compatible.full_match() - Variable-ratio full matchingcem_match() - Coarsened exact matchingsubclass_match() - Propensity score
subclassificationmatch_data() - Unified analysis-ready outputas_matchit() - Convert to MatchIt formatratio parameter
in match_couples() and greedy_couples().
Matches k control units to each treated unit by replicating the cost
matrix, then deduplicates assignments.replace
parameter. Each treated unit independently selects its nearest control,
allowing controls to be reused across multiple treated units.ps_match() function wraps
match_couples() with logistic regression:
glm objectcardinality_match() function
maximizes sample size subject to balance constraints:
max_std_diff (default: 0.1 for
excellent balance)batch_fractionsensitivity_analysis() function
implements Rosenbaum bounds:
print(), summary(),
plot()autoplot() methods for ggplot2-based
visualizations (requires ggplot2):
autoplot.matching_result(): histogram, density, or ecdf
of distancesautoplot.balance_diagnostics(): love plot, histogram,
or variance ratio plotautoplot.sensitivity_analysis(): gamma vs p-value
curvesummary.matching_result() now
reports match rate and distance percentilesps_match() - Propensity score matching with logit
calipercardinality_match() - Balance-constrained cardinality
matchingsensitivity_analysis() - Rosenbaum bounds sensitivity
analysisselect() in vignettes by
using explicit dplyr::select() to prevent masking by MASS
or other packagesThe package now includes intelligent preprocessing to improve matching quality:
auto_scale parameter in
match_couples() and greedy_couples() enables
automatic preprocessingpreprocess_matching_vars() function for manual
preprocessing controlComprehensive tools to assess matching quality:
balance_diagnostics() function
computes multiple balance metrics:
balance_table() creates publication-ready formatted
tablesCreate analysis-ready datasets directly from matching results:
join_matched() function automates
data preparation:
left_vars and
right_vars parameters_left,
_right) for overlapping columnspair_id, distance,
block_idaugment() method for
tidymodels integration:
join_matched() parametersinclude_distance - Include/exclude matching
distanceinclude_pair_id - Include/exclude sequential pair
IDsinclude_block_id - Include/exclude block
identifiersleft_id and
right_idPerformance optimization for exploring multiple matching strategies:
compute_distances() function
precomputes and caches distance matrices:
join_matched()distance_object):
match_couples() and
greedy_couples()match_couples(dist_obj, max_distance = 5)update_constraints():
max_distance or calipers without
recomputing distancesmatch_couples(left, right = NULL, vars = NULL, ...)Speed up blocked matching with multi-core processing:
parallel parameter in
match_couples() and greedy_couples():
parallel = TRUE for automatic
configurationparallel = "multisession" or other
future planfuture package:
Like testthat, couplr makes errors light, memorable, and helpful with couple-themed messages:
check_costs parameter (default:
TRUE) in match_couples() and
greedy_couples():
FALSE to skip checks in production codediagnose_distance_matrix():
options(couplr.emoji = FALSE) if preferredpreprocess_matching_vars() - Main preprocessing
orchestratorbalance_diagnostics() - Comprehensive balance
assessmentbalance_table() - Formatted balance tables for
reportingjoin_matched() - Create analysis-ready datasets from
matching resultsaugment.matching_result() - Broom-style interface for
joined datacompute_distances() - Precompute and cache distance
matricesupdate_constraints() - Modify constraints on distance
objectsis_distance_object() - Type checking for distance
objectsdiagnose_distance_matrix() - Comprehensive distance
diagnosticscheck_cost_distribution() - Check for distribution
problemsexamples/auto_scale_demo.R - 5 preprocessing
demonstrationsexamples/balance_diagnostics_demo.R - 6 balance
diagnostic examplesexamples/join_matched_demo.R - 8 joined dataset
demonstrationsexamples/distance_cache_demo.R - Distance caching and
reuse examplesexamples/parallel_matching_demo.R - 7 parallel
processing examplesexamples/error_messages_demo.R - 10 fun error message
demonstrationsThe package has been renamed from lapr to couplr to better reflect its purpose as a general pairing and matching toolkit.
couplr = Optimal pairing and matching via linear assignment
First official stable release with clean, well-organized codebase.
morph_* naming prefixassignment() (low-level) +
lap_solve() (tidy)src/core/ - Utilities and headerssrc/interface/ - Rcpp exportssrc/solvers/ - 14 LAP algorithmssrc/gabow_tarjan/ - Gabow-Tarjan solversrc/morph/ - Image morphingHungarian, Jonker-Volgenant, Auction (3 variants), SAP/SSP, SSAP-Bucket, Cost-scaling, Cycle-cancel, Gabow-Tarjan, Hopcroft-Karp, Line-metric, Brute-force, Auto-select
✅ Tidy tibble interface ✅ Matrix & data frame inputs
✅ Grouped data frames ✅ Batch solving + parallelization ✅ K-best
solutions (Murty, Lawler) ✅ Rectangular matrices ✅ Forbidden
assignments (NA/Inf) ✅ Maximize/minimize ✅ Pixel morphing
visualization
lap_solve() - Main tidy interfacelap_solve_batch() - Batch solvinglap_solve_kbest() - K-best solutionsassignment() - Low-level solverget_total_cost(),
as_assignment_matrix(), etc.pixel_morph(),
pixel_morph_animate()Development history under “lapr” available in git log before v1.0.0.