---
title: 'Matching Algorithms in R and C++: An Introduction to matchingR'
output: knitr:::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Matching Algorithms in R and C++: An Introduction to matchingR}
---

```{r, results = "hide", echo=FALSE, message = FALSE}
library(matchingR)
```

# Introduction

`matchingR` is an R package which quickly computes the [Gale-Shapley algorithm](https://www.jstor.org/stable/2312726) [@gale1962college], [Irving's algorithm for the stable roommate problem](https://www.sciencedirect.com/science/article/pii/0196677485900331/) [@irving1985roommates], and the [top trading cycle algorithm](https://www.sciencedirect.com/science/article/abs/pii/0304406874900330/) [@shapley1973cores] for large matching markets. The package provides functions to compute the solutions to the
  [stable marriage problem](https://en.wikipedia.org/wiki/Stable_matching), the
  [college admission problem](https://en.wikipedia.org/wiki/Hospital_resident), the
  [stable roommates problem](https://en.wikipedia.org/wiki/Stable_roommates_problem), and the
  [house allocation problem](https://web.stanford.edu/~niederle/HouseAllocation.pdf).

The package may be useful when the number of market participants is large or when many matchings need to be computed (e.g., for simulation or estimation purposes). It has been used in practice to compute the Gale-Shapley stable matching with 30,000 participants on each side of the market.

Matching markets are common in practice and widely studied by economists. Popular examples include

 * the [National Resident Matching Program](https://www.nrmp.org/) that matches graduates from medical school to residency programs at teaching hospitals throughout the United States
 * the matching of students to schools including the [New York City High School Match](https://www.jstor.org/stable/4132848) or the the [Boston Public School Match](https://www.jstor.org/stable/4132849) (and many more)
 * the matching of kidney donors to recipients in [kidney exchanges](https://www.jstor.org/stable/4132851).

# Two-sided Matching Markets: Gale-Shapley Algorithm

Consider the following marriage market: There are `N` men and `N` women. Each man, `m`, receives utility `uM(w, m)` from a match with woman `w`, and similarly each woman receives a payoff of `uW(m, w)` from being matched with a man.

A matching assigns men to women such that each man is assigned to one woman and each woman is assigned to one man. A matching is **stable** if there is no man and woman who would jointly prefer to be matched to each other over their current spouses. In other words, a matching is stable if there are no pairs `(m, w'), (m', w)`, such that `m` is matched with `w'`, `m'` is matched with `w`, and both `uW(m, w) > uW(m', w)` and `uM(m, w) > uM(m, w')`.

For example, we might have preferences for men given by
```{r}
uM = matrix(c(1.0, 0.5, 0.0,
              0.5, 0.0, 0.5,
              0.0, 1.0, 1.0), nrow = 3, ncol = 3, byrow = TRUE)
```
```{r, echo=FALSE}
dimnames(uM) = list(rows = c('Woman 1', 'Woman 2', 'Woman 3'), cols = c('Man 1', 'Man 2', 'Man 3'))
uM
```
In this example, man `1` receives a payoff of `1.0` from being matched to woman `1`, a payoff of `0.5` from being matched to woman `2` and a payoff of `0.0` from being matched to woman 3 (same logic applies to men `2` and `3`). Similarly, we might have preferences for women given by
```{r}
uW = matrix(c(0.0, 1.0, 0.0,
              0.5, 0.0, 0.5,
              1.0, 0.5, 1.0), nrow = 3, ncol = 3, byrow = TRUE)
```
```{r, echo=FALSE}
dimnames(uW) = list(rows = c('Man 1', 'Man 2', 'Man 3'), cols = c('Woman 1', 'Woman 2', 'Woman 3'))
uW
```
Here, columns in the matrix correspond to women, rows to men. In this example, woman `1` receives a payoff of `0.0` from being matched to man `1`, a payoff of `0.5` from being matched to man `2` and a payoff of `1.0` from being matched to man 3 (same logic applies to women `2` and `3`).

Instead of using payoff matrices, we can also represent preferences using preference orderings. The preference ordering that corresponds to `uM` is
```{r}
prefM = matrix(c(1, 3, 3,
                 2, 1, 2,
                 3, 2, 1), nrow = 3, ncol = 3, byrow = TRUE)
```
```{r, echo=FALSE}
dimnames(prefM) = list(rows = c('Rank 1', 'Rank 2', 'Rank 3'), cols = c('Man 1', 'Man 2', 'Man 3'))
prefM
```
`prefM` states that man `1` prefers woman `1` over woman `2` over woman `3`, etc. The preference ordering that corresponds to `uW` is given by
```{r}
prefW = matrix(c(3, 1, 3,
                 2, 3, 2,
                 1, 2, 1), nrow = 3, ncol = 3, byrow = TRUE)
```
```{r, echo=FALSE}
dimnames(prefW) = list(rows = c('Rank 1', 'Rank 2', 'Rank 3'), cols = c('Woman 1', 'Woman 2', 'Woman 3'))
prefW
```
The matching algorithm discussed below can take either payoff matrices of the type `uM` and `uW` or preference orderings of the type `prefM` and `prefW` as arguments.

The Gale-Shapley algorithm works as follows: Single men ("the proposers") sequentially make proposals to each of their most preferred available women ("the reviewers"). A woman can hold on to at most one proposal at a time. A *single* woman will accept any proposal that is made to her. A woman who already has a proposal will reject any proposal she values less than her current proposal in hand. If a woman receives a proposal from a man that she values more than her current proposal, she will accept the proposal and her previous match will rejoin the line of proposers. This process continues until all men are matched to all women.

For the preferences specified in `uM` and `uW`, we can compute the Gale-Shapley Algorithm by hand. Initially, all men are single.

 1. * Man `1` proposes to woman `1`, his most-preferred choice.
    * Unmatched men: `2`, `3`
 2. * Man `2` proposes to woman `3`, his most-preferred choice.
    * Unmatched men: `3`
 3. * Man `3` proposes to woman `3`, his most-preferred choice.
    * Woman `3` now dumps man `2`.
    * Unmatched men: `2`
 4. * Man `2` proposes to woman `1`, his most-preferred available choice.
    * Woman `1` now dumps man `1`.
    * Unmatched men: `1`
 5. * Man `1` proposes to woman `2`, his most-preferred available choice.
    * All men are now matched.

 The man-optimal stable matching is therefore:

 |Man | Woman|
 |-------|---------|
 | 1 |  2 |
 |  2 |  1 |
 |  3 |  3 |


The package computes the Gale-Shapley algorithm using the function `galeShapley.marriageMarket`:
```{r}
matching = galeShapley.marriageMarket(uM, uW)
```
Note that we can obtain equivalent results when we use `prefM` and `prefW` as arguments:
```{r}
matching = galeShapley.marriageMarket(proposerPref = prefM, reviewerPref = prefW)
```
The function `galeShapley.marriageMarket` returns a list `matching` that includes the vectors `proposals` and `engagements` with the final proposals and engagements, respectively. These two vectors contain the same information (i.e. they tell us who is matched with whom). For the example above, the vector of proposals contains
```{r, echo=FALSE}
dimnames(matching$proposals) = list(rows = c("Man 1", "Man 2", "Man 3"),
                                    cols = c("Proposed to Woman"))
```
```{r}
matching$proposals
```
The first element in the vector tells us that man `1` is matched with woman `2`. Man `2` is matched to woman `1`, and man `3` is matched to woman `3`. The vector of engagement contains
```{r, echo=FALSE}
dimnames(matching$engagements) = list(rows = c("Woman 1", "Woman 2", "Woman 3"),
                                      cols = c("Engaged to Man"))
```
```{r}
matching$engagements
```
The first element in the vector tells us that woman `1` is matched to man `2`, woman `2` will be matched to man `1`, and woman `3` will be matched to man `3`.

We can then check if the computed matching is stable using the function `checkStability`. To check if a matching is stable, we check for each assignment `(m,w)` if there is some other woman `w'` that man `m` would rather be matched with and who would rather be matched to man `m`. This function will return `true` if the matching is stable and `false` otherwise.
```{r}
galeShapley.checkStability(uM, uW, matching$proposals, matching$engagements)
```
For the simple 3-by-3 example, we can perform this check by hand:

 * Man `1` is matched to woman `2`, his second-most preferred choice. His most preferred choice is woman `1`. Woman `1` is matched with man `2` who she prefers over man `1`. Thus man `1` cannot do better than woman `2`.
 * Man `2` is matched to woman `1`, his second-most preferred choice. His most preferred woman is woman `3`, who is     matched with man `3`. Since man `3` is her most-preferred choice, man `2` cannot do better than woman `1`.
 * Man `3` is matched to women `3`, his most preferred choice, so he cannot do better.

Thus, this matching is stable.

The following examples illustrate some additional features of this package.

### Example: Marriage Market
The following is an example of `galeShapley.marriageMarket` with different numbers of participants on each side of the market. There are 2,500 women and 2,000 men. By construction, 500 men will remain unmatched. We randomly generate payoff matrices `uM` and `uW` which are drawn from a uniform distribution (`runif`). We then compute the male-optimal (i.e. men are proposing) and the female-optimal (i.e. woman are proposing) matching.
```{r}
# set seed
set.seed(1)
# set number of men
nmen = 2500
# set number of women
nwomen = 2000
# generate preferences
uM = matrix(runif(nmen*nwomen), nrow = nwomen, ncol = nmen)
uW = matrix(runif(nmen*nwomen), nrow = nmen, ncol = nwomen)
# male-optimal matching
resultsM = galeShapley.marriageMarket(uM, uW)
str(resultsM)
galeShapley.checkStability(uM, uW, resultsM$proposals, resultsM$engagements)
# female-optimal matching
resultsW = galeShapley.marriageMarket(uW, uM)
str(resultsW)
galeShapley.checkStability(uW, uM, resultsW$proposals, resultsW$engagements)
```

### Example: College Admissions Problem
The following is an example of `galeShapley.collegeAdmissions` where 1000 students get matched to 400 colleges, where each college has two slots. By construction, 200 students will remain unmatched. We draw students' and colleges' preferences, `uStudents` and `uColleges`, respectively, by from a uniform distribution.
```{r}
# set seed
set.seed(1)
# set number of students
nstudents = 1000
# set number of colleges
ncolleges = 400
# generate preferences
uStudents = matrix(runif(ncolleges*nstudents), nrow = ncolleges, ncol = nstudents)
uColleges = matrix(runif(nstudents*ncolleges), nrow = nstudents, ncol = ncolleges)
# student-optimal matching
results = galeShapley.collegeAdmissions(studentUtils =  uStudents, collegeUtils =  uColleges, slots = 2)
str(results)
# check if matching is stable
galeShapley.checkStability(uStudents, uColleges, results$matched.students, results$matched.colleges)
```

# One-sided Matching Markets: Irving's Algorithm

This package implements the algorithm by @irving1985roommates for one-sided matching markets.

Consider the following example: A set of `n` potential roommates, each with ranked preferences over all the other potential roommates, are to be matched to rooms, two roommates per room. A matching is **stable** if there is no roommate `r1` that would rather be matched to some other roommate `d2` than to his current roommate `r2` and the other roommate `d2` would rather be matched to `r1` than to his current roommate `d1`.

Preferences of potential roommates are summarized by an $n-1 \times n$ dimensional matrix, e.g., if $n = 6$,
```{r}
pref = matrix(c(3, 6, 2, 5, 3, 5,
                4, 5, 4, 2, 1, 1,
                2, 4, 5, 3, 2, 3,
                6, 1, 1, 6, 4, 4,
                5, 3, 6, 1, 6, 2), nrow = 5, ncol = 6, byrow = TRUE)
```
Column `i` represents the preferences of the `i`th roommate, and row `j` represents the ranking of the roommate whose index is encoded in that row. For example, in the above preference matrix, roommate `1` most prefers to be matched with roommate `3`, followed by `4`, followed by `2`.

The function `roommate.checkPreferences` checks if a given preference order is complete, i.e. if all preferences are fully specified. If the preference order is complete, it will return the proper preference order using `C++` style indexing (beginning at `0` instead of at `1`). If the preference order is incomplete, the function will return an error.
```{r}
roommate.checkPreferences(pref)
```

The algorithm proceeds in two phases.

In phase 1, potential roommates take turns sequentially proposing to the other roommates. Each roommate who is proposed to can accept or reject the proposal. A roommate accepts if he currently has made no better proposal which was accepted to another roommate. If a roommate has accepted a proposal, and then receives a better proposal, he rejects the old proposal and substitutes in the new proposal.

In the above example,

1. Roommate `1` begins by proposing to roommate `3`, his most preferred roommate. `3`, having no better offers, accepts.
2. `2` proposes to `6`, who accepts.
3. `3` proposes to `2`, who accepts.
4. `4` proposes to `5`, who accepts.
5. `5` proposes to `3`, who accepts. `3` cancels his proposal from `1`.
6. `1`, having no proposal, proposes to `4`, who accepts.
7. `6` proposes to `5`, who rejects, having a better proposal from `4`.
8. `6` proposes to `1`, who accepts.

In phase 2, we begin by eliminating all potential roommate matches which are worse than the current proposals held. For example, in the above example, `3` has a proposal from `5`, and so we eliminate `1` and `6` from `3`'s column, and symmetrically eliminate `3` from `1` and `6`'s column. This results in the following 'reduced' preference listing:
```
   6, 2, 5, 3,
4, 5, 4, 2,    1
2, 4, 5, 3, 2,
6, 1,    6, 4, 4
   3,    1,    2
```
These preferences form what is called a 'stable' table, or, 's-table'. ('Stable' for short.) The defining characteristic of a stable table is that if `i` is the most preferred roommate on `j`s list, then `j` is the least preferred roommate on `i`s list. For example, `1` most prefers `4`, but `4` least prefers `1`.

The algorithm proceeds by finding and eliminating 'rotations'. A rotation is a sequence of pairs of roommates, such that there is a distinct roommate in the first position of each pair, the second roommate in each pair least prefers the roommate he is paired with, the first roommate in each pair most prefers the roommate he is paired with, and finally, the first roommate in each pair ranks the second roommate in the following pair second (modulo the number of pairs, that is, the first roommate in the last pair ranks the second roommate in the first pair second) Once a rotation has been identified, removing it results in another stable table.

For example, `(1, 4), (3, 2)` is a rotation in the above table, because `1` loves `4`, `3` loves `2`, `4` hates `1`, `2` hates `3`, `2` is second on `1`s list, and `4` is second on `3`'s list. Eliminating this rotation involves `2` rejecting `3`, `4` rejecting `1`, and then we remove every successive potential roommate as well to preserve the stable table property, resulting in
```
   6,    5, 3,
   5, 4, 2,    1
2, 4, 5, 3, 2,
6, 1,    6, 4, 4
               2
```
A further rotation is `(1, 2), (2, 6), (4, 5)`. Eliminating it yields
```
            3,
   5, 4, 2,    1
   4, 5, 3, 2,
6,
```
A final rotation is `(2, 5), (3, 4)`. Eliminating it yields
```
            3,
         2,    1
   4, 5,
6,

```
Therefore, a stable matching is for `1` and `6` to match, `2` and `4` to match, and `3` and `5` to match.
```{r}
results = roommate(pref = pref)
results
```
The function `roommate.checkStability` can be used to check if the resulting matching is stable:
```{r}
roommate.checkStability(pref = pref, matching = results)
```

### Example: Roommate problem

The function `roommate` can also be called using a payoff matrix, `u`, to specify preferences. In `u`, the element `[i,j]` refers to the payoff that agent `j` gets from being matched to agent `i`. The main diagonal of this matrix contains no information and will be removed by the algorithm.

```{r}
# generate preferences
N = 10
u = matrix(runif(N^2),  nrow = N, ncol = N)
results = roommate(utils = u)
results
roommate.checkStability(utils = u, matching = results)
```

### Example: Roommate problem when no stable matching exists

Note that in the roommate problem, existence of a stable matching is not guaranteed. When no stable matching can be found, the function `roommate` returns `NULL`.
```{r}
set.seed(1)
N = 512
u = matrix(runif(N^2),  nrow = N, ncol = N)
results = roommate(utils = u)
print(results)
```

# Kidneys and Housing: The Top Trading Cycle Algorithm

This package implements the top trading cycle algorithm.

Consider the following problem: A set of $n$ agents each currently own their own home, and have preferences over the homes of other agents. The problem is to trade the homes between the agents in such a way so that no two agents want to swap homes.

Preferences of agents are summarized by an $n \times n$ dimensional matrix, e.g., if $n = 4$,

```{r}
pref = matrix(c(4, 4, 2, 4,
                2, 1, 1, 1,
                1, 2, 3, 3,
                3, 3, 4, 2), nrow = 4, ncol = 4, byrow = TRUE)
```

Column $i$ represents the preferences of the $i$th agent, and row $j$ represents the ranking of the roommate whose index is encoded in that row. For example, in the above preference matrix, agent `1` most prefers the home of agent `4`, followed by `2`, followed by `1`, followed by `3`.

Roughly speaking, the top trading cycle proceeds by identifying cycles of agents, then eliminating those cycles until no agents remain. A cycle is a sequence of agents such that each agent most prefers the next agent's home (out of the remaining unmatched agents), and the last agent in the sequence most prefers the first agent in the sequence's home.

```
4,  4,  2,  4
2,  1,  1,  1
1,  2,  3,  3
3,  3,  4,  2
```

For example, for the above preference matrix, when all the agents are unmatched, the only rotation is $\{4\}$, representing the fact that agent `4` most prefers his own house. Therefore, the algorithm begins by matching agent `4` to himself, and then removing him from the pool:

```
        2
2,  1,  1
1,  2,  3
3,  3,
```

Now, a rotation is $\{1, 2\}$, because `1` most prefers `2`s house, and `2` most prefers `1`s house. So agents `1` and `2` will swap homes, leaving agent `3` all by himself.

```


        3

```

Therefore, the final matching is that agent `1` swaps with agent `2`, and agents `3` and `4` keep their own homes.

```{r}
results = toptrading(pref = pref)
results
```
The function `toptrading.checkStability` can be used to check if the resulting allocation is in fact stable:
```{r}
toptrading.checkStability(pref = pref, matchings = results)
```

# Literature