Gale-Shapley Algorithm for R

[VictorDelgado]

Dear R-helpers,

I'm not a speciallist in writing complex functions, and the function still very rusty (any kind of suggestions are very welcome). I want to implement Gale-Shapley algorithm for R Language. It is based on Gale and Shapley (1962), and it has evolved to several applications in many languages, C++, JAVA, Perl, SQL, and so on. I manage to edit one version of it  to R.2.13.1. So, I ask if it was allready implemented (I couldn't find any on R topic), and if there is models and manners to make it more efficiently, add errors check, options, etc.

At Berkeley's MathSite there is a very straighfoward example of the algorithm and its steps.

My implementation follow the principle:

1. All men (or women) seeks for their best partner.

2. If there is no tie in a column (or row), stop.

3. If there is a tie, removes the worst-partners-tied and seek again the second-best (till n-best) alternative.

The function is working right up to 6x6 matrices. But it needs a lot of improvement.

Here it is the "gsa" function:

gsa(m, n, preference.row, preference.col, first)

###
Where:

m: for number of rows
n:  number of columns
preference.row: matrix with preference ordered in its positions by row (see example).
preference.col: matrix with preference ordered in its positions by column (see example).
first: Who is the first to propose (1 to men, 2 to women).
########

gsa <- function(m, n, preference.row, preference.col, first)
{
# m: number of rows (men)
# n: number of columns (women)
# first 1 for row (men); and 2 for column (women)
#
# Two Auxiliary functions:
# 1:
min.n <- function(x,n,value=TRUE){
s <- sort(x, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}

# 2:

max.n <- function(x,n,value=TRUE){
s <- sort(x, decreasing=TRUE, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}
#############################################################

s <- NULL
test_s <-NULL
loop <- 2 # O loop é necessário a partir do 2.
step.1 <- matrix(0,ncol=n, nrow=m)
step.2 <- matrix(0,ncol=n, nrow=m)
store <- NULL
r <- NULL

# Men proposing first:

if (first==1)
        {
step.1 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.1[i,][preference.row[i,]==min.n(preference.row[i,],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.1[,position1[k]]==1)
position3[[k]] <- which(preference.col[,position1[k]]>min(preference.col[position2[[k]],position1[k]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.1[position3[[k]],position1[k]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.1[,t])==0){0}else{which(step.1[,t]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(sort(store))
for (j in position5){step.1[j,][preference.row[j,]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
                                        }else{
step.1 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.1[i,][preference.row[i,]==min(preference.row[i,])] <- 1}
return(step.1)}
loop <- loop + 1
                        } #end of while
        }

# Women proposing first:

if (first==2)
        {
step.2 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.2[,i][preference.col[,i]==min.n(preference.col[,i],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.2[position1[k],]==1)
position3[[k]] <- which(preference.row[position1[k],]>min(preference.row[position1[k],position2[[k]]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.2[position1[k],position3[[k]]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.2[t,])==0){0}else{which(step.2[t,]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(store)
for (j in position5){step.2[,j][preference.col[,j]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
                                        }else{
step.2 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.2[i,][preference.col[,i]==min(preference.col[,i])] <- 1}
step.2}
loop <- loop + 1
                        } # End of 2nd while
        }
if (first==1) {print(step.1)}
if (first==2) {print(step.2)}
}

#####################
# Here it goes one 4x4 example:

m <- seq(1:4)
n <- seq(1:4)
preference.row <- matrix(0,ncol=length(m), nrow=length(m))
preference.col <- matrix(0,ncol=length(n), nrow=length(n))

for (i in 1:length(m))
{
preference.row[i,] <- sample(m, size=length(m), rep=FALSE)
preference.col[,i] <- sample(n, size=length(n), rep=FALSE) # Note a orientação por coluna!
}
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=2)

# The result is a zero-one matrix which indicates blocking pairs.
############################################################

Thank you, and please let me know, any bugs and improvements.

[VictorDelgado]

VictorDelgado wrote

Dear R-helpers,

I'm not a speciallist in writing complex functions, and the function still very rusty (any kind of suggestions are very welcome). I want to implement Gale-Shapley algorithm for R Language. It is based on Gale and Shapley (1962), and it has evolved to several applications in many languages, C++, JAVA, Perl, SQL, and so on. I manage to edit one version of it  to R.2.13.1. So, I ask if it was allready implemented (I couldn't find any on R topic), and if there is models and manners to make it more efficiently, add errors check, options, etc.

At Berkeley's MathSite there is a very straighfoward example of the algorithm and its steps.

My implementation follow the principle:

1. All men (or women) seeks for their best partner.

2. If there is no tie in a column (or row), stop.

3. If there is a tie, removes the worst-partners-tied and seek again the second-best (till n-best) alternative.

The function is working right up to 6x6 matrices. But it needs a lot of improvement.

Here it is the "gsa" function:

gsa(m, n, preference.row, preference.col, first)

###
Where:

m: for number of rows
n:  number of columns
preference.row: matrix with preference ordered in its positions by row (see example).
preference.col: matrix with preference ordered in its positions by column (see example).
first: Who is the first to propose (1 to men, 2 to women).
########

gsa <- function(m, n, preference.row, preference.col, first)
{
# m: number of rows (men)
# n: number of columns (women)
# first 1 for row (men); and 2 for column (women)
#
# Two Auxiliary functions:
# 1:
min.n <- function(x,n,value=TRUE){
s <- sort(x, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}

# 2:

max.n <- function(x,n,value=TRUE){
s <- sort(x, decreasing=TRUE, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}
#############################################################

s <- NULL
test_s <-NULL
loop <- 2 # O loop é necessário a partir do 2.
step.1 <- matrix(0,ncol=n, nrow=m)
step.2 <- matrix(0,ncol=n, nrow=m)
store <- NULL
r <- NULL

# Men proposing first:

if (first==1)
        {
step.1 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.1[i,][preference.row[i,]==min.n(preference.row[i,],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.1[,position1[k]]==1)
position3[[k]] <- which(preference.col[,position1[k]]>min(preference.col[position2[[k]],position1[k]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.1[position3[[k]],position1[k]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.1[,t])==0){0}else{which(step.1[,t]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(sort(store))
for (j in position5){step.1[j,][preference.row[j,]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
                                        }else{
step.1 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.1[i,][preference.row[i,]==min(preference.row[i,])] <- 1}
return(step.1)}
loop <- loop + 1
                        } #end of while
        }

# Women proposing first:

if (first==2)
        {
step.2 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.2[,i][preference.col[,i]==min.n(preference.col[,i],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.2[position1[k],]==1)
position3[[k]] <- which(preference.row[position1[k],]>min(preference.row[position1[k],position2[[k]]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.2[position1[k],position3[[k]]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.2[t,])==0){0}else{which(step.2[t,]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(store)
for (j in position5){step.2[,j][preference.col[,j]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
                                        }else{
step.2 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.2[i,][preference.col[,i]==min(preference.col[,i])] <- 1}
step.2}
loop <- loop + 1
                        } # End of 2nd while
        }
if (first==1) {print(step.1)}
if (first==2) {print(step.2)}
}

#####################
# Here it goes one 4x4 example:

m <- seq(1:4)
n <- seq(1:4)
preference.row <- matrix(0,ncol=length(m), nrow=length(m))
preference.col <- matrix(0,ncol=length(n), nrow=length(n))

for (i in 1:length(m))
{
preference.row[i,] <- sample(m, size=length(m), rep=FALSE)
preference.col[,i] <- sample(n, size=length(n), rep=FALSE) # Note a orientação por coluna!
}
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=2)

# The result is a zero-one matrix which indicates blocking pairs.
############################################################

Thank you, and please let me know, any bugs and improvements.

I have implemented some changes to see "loop" iterations:
loop <- 1

if (first==1) {print(step.1)}
if (first==2) {print(step.2)}
print(loop)
}

And just added some Examples from Gale and Shapley (1962) College Admissions And the Stability of Marriage:

# 1:

m1 <- c(1,2,3); m2 <- c(3,1,2); m3 <- c(2,3,1)
n1 <- c(3,1,2) ;n2 <- c(2,3,1); n3 <- c(1,2,3)
preference.row <- matrix(c(m1, m2, m3), ncol=3, byrow=TRUE)
preference.col <- matrix(c(n1, n2, n3), ncol=3)
gsa(m = 3, n = 3, preference.row = preference.row, preference.col = preference.col, first=1)
gsa(m = 3, n = 3, preference.row = preference.row, preference.col = preference.col, first=2)

# 2 :

m1 <- c(1,2,3,4) ; m2 <- c(1,4,3,2); m3 <- c(2,1,3,4); m4 <- c(4,2,3,1)
n1 <- c(3,4,2,1); n2 <- c(3,1,4,2); n3 <- c(2,3,4,1); n4 <- c(3,2,1,4)
preference.row <- matrix(c(m1, m2, m3, m4), ncol=4, byrow=TRUE)
preference.col <- matrix(c(n1, n2, n3, n4), ncol=4)
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=1)
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=2)

#3:

m1 <- c(1,2,3,4); m2 <- c(1,2,3,4); m3 <- c(3,1,2,4); m4 <- c(2,3,1,4)
n1 <- c(3,4,1,2); n2 <- c(2,3,4,1); n3 <- c(1,2,3,4); n4 <- c(3,4,2,1)
preference.row <- matrix(c(m1, m2, m3, m4), ncol=4, byrow=TRUE)
preference.col <- matrix(c(n1, n2, n3, n4), ncol=4)
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=1)
gsa(m = 4, n = 4, preference.row = preference.row, preference.col = preference.col, first=2)

[VictorDelgado]

VictorDelgado wrote

gsa <- function(m, n, preference.row, preference.col, first)
{
# m: number of rows (men)
# n: number of columns (women)
# first 1 for row (men); and 2 for column (women)
#
# Two Auxiliary functions:
# 1:
min.n <- function(x,n,value=TRUE){
s <- sort(x, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}

# 2:

max.n <- function(x,n,value=TRUE){
s <- sort(x, decreasing=TRUE, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}
#############################################################

s <- NULL
test_s <-NULL
loop <- 2 # O loop é necessário a partir do 2.
step.1 <- matrix(0,ncol=n, nrow=m)
step.2 <- matrix(0,ncol=n, nrow=m)
store <- NULL
r <- NULL

# Men proposing first:

if (first==1)
        {
step.1 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.1[i,][preference.row[i,]==min.n(preference.row[i,],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.1[,position1[k]]==1)
position3[[k]] <- which(preference.col[,position1[k]]>min(preference.col[position2[[k]],position1[k]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.1[position3[[k]],position1[k]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.1[,t])==0){0}else{which(step.1[,t]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(sort(store))
for (j in position5){step.1[j,][preference.row[j,]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
                                        }else{
step.1 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.1[i,][preference.row[i,]==min(preference.row[i,])] <- 1}
return(step.1)}
loop <- loop + 1
                        } #end of while
        }

# Women proposing first:

if (first==2)
        {
step.2 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.2[,i][preference.col[,i]==min.n(preference.col[,i],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.2[position1[k],]==1)
position3[[k]] <- which(preference.row[position1[k],]>min(preference.row[position1[k],position2[[k]]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.2[position1[k],position3[[k]]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.2[t,])==0){0}else{which(step.2[t,]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(store)
for (j in position5){step.2[,j][preference.col[,j]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
                                        }else{
step.2 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.2[i,][preference.col[,i]==min(preference.col[,i])] <- 1}
step.2}
loop <- loop + 1
                        } # End of 2nd while
        }
if (first==1) {print(step.1)}
if (first==2) {print(step.2)}
}

I Just have fixed some problems with the first function. Now it's running with 100x100 (random preferences) matrices. The function still needing some simplification.

gsa <- function(m, n, preference.row, preference.col, first)
{
#
########### TWO VERY USEFUL AUXILIARITY FUNCTIONS:
#
# Returns the n-esim minimun
# If value=TRUE it gives you the value, otherwise it returns the position.

min.n <- function(x,n,value=TRUE){
s <- sort(x, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}

# Same Function for max:

max.n <- function(x,n,value=TRUE){
s <- sort(x, decreasing=TRUE, index.return=TRUE)
if(value==TRUE){s$x[n]} else{s$ix[n]}}
#############################################################

# 1 for men proposing; 2 for women.
s <- NULL
test_s <-NULL
loop <- 1 # Contagem das iterações.
step.1 <- matrix(0,ncol=n, nrow=m)
step.2 <- matrix(0,ncol=n, nrow=m)
store <- NULL
r <- NULL

# Men proposing:

if (first==1)
        {
step.1 <- matrix(0,ncol=n, nrow=m)
for (i in 1:m)
                {
step.1[i,][preference.row[i,]==min.n(preference.row[i,],n=1)] <- 1
                }
for (i in 1:m){s[i] <- sum(step.1[,i])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:n
for (k in 1:length(position1)){position2[[k]] <- which(step.1[,position1[k]]==1)
position3[[k]] <- which(preference.col[,position1[k]]>min(preference.col[position2[[k]],position1[k]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.1[position3[[k]],position1[k]] <- 0}
for (t in 1:n){position4[t] <- if(sum(step.1[,t])==0){0}else{which(step.1[,t]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(sort(store))
for (j in position5){step.1[j,][preference.row[j,]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.1[,i])}
test_s <- s>1
                                        }else{
step.1 <- matrix(0,ncol=m, nrow=n)
for (i in 1:n){step.1[i,][preference.row[i,]==min(preference.row[i,])] <- 1}
return(step.1)}
loop <- loop + 1
                        } #end of while
        }

# Women proposing:

if (first==2)
        {
step.2 <- matrix(0,ncol=n, nrow=m)
for (i in 1:n)
                {
step.2[,i][preference.col[,i]==min.n(preference.col[,i],n=1)] <- 1
                }
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
while (any(test_s==TRUE)==TRUE)
                        {
if (any(test_s==TRUE)==TRUE) {
position1 <- which(s>1)
position2 <- vector('list')
position3 <- vector('list')
position4 <- NULL
position5 <- 1:m
for (k in 1:length(position1)){position2[[k]] <- which(step.2[position1[k],]==1)
position3[[k]] <- which(preference.row[position1[k],]>min(preference.row[position1[k],position2[[k]]]))
x <- which(position3[[k]]%in%position2[[k]])
position3[[k]] <- position3[[k]][x]
step.2[position1[k],position3[[k]]] <- 0}
for (t in 1:m){position4[t] <- if(sum(step.2[t,])==0){0}else{which(step.2[t,]==1)}}
position4 <- position4[position4 >0]
position5 <- position5[-position4]
store <- append(position5, store)
r <- rle(sort(store))
for (j in position5){step.2[,j][preference.col[,j]==r$lengths[r$values==j]+1] <- 1}
for (i in 1:n){s[i] <- sum(step.2[i,])}
test_s <- s>1
                                        }else{
step.2 <- matrix(0,ncol=m, nrow=n)
for (i in 1:m){step.2[i,][preference.col[,i]==min(preference.col[,i])] <- 1}
return(step.2)}
loop <- loop + 1
                        } # End of 2nd while
        }
if (first==1) {print(step.1)}
if (first==2) {print(step.2)}
print(loop)
}