scripts/builtin/stableMarriage.dml

#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

# This script computes a solution for the stable marriage problem.
#
# result description:
#
# If cell [i,j] is non-zero, it means that acceptor i has matched with
# proposer j. Further, if cell [i,j] is non-zero, it holds the preference
# value that led to the match.
# Proposers.mtx:
# 2.0,1.0,3.0
# 1.0,2.0,3.0
# 1.0,3.0,2.0
# 
# Since ordered=TRUE, this means that proposer 1 (row 1) likes acceptor 2
# the most, followed by acceptor 1 and acceptor 3.
# If ordered=FALSE, this would mean that proposer 1 (row 1) likes acceptor 3
# the most (since the value at [1,3] is the row max),
# followed by acceptor 1 (2.0 preference value) and acceptor 2 (1.0 preference value).
# 
# Acceptors.mtx:
# 3.0,1.0,2.0
# 2.0,1.0,3.0
# 3.0,2.0,1.0
# 
# Since ordered=TRUE, this means that acceptor 1 (row 1) likes proposer 3
# the most, followed by proposer 1 and proposer 2.
# If ordered=FALSE, this would mean that acceptor 1 (row 1) likes proposer 1
# the most (since the value at [1,1] is the row max),
# followed by proposer 3 (2.0 preference value) and proposer 2
# (1.0 preference value).
# 
# Output.mtx (assuming ordered=TRUE):
# 0.0,0.0,3.0
# 0.0,3.0,0.0
# 1.0,0.0,0.0
# 
# Acceptor 1 has matched with proposer 3 (since [1,3] is non-zero) at a
# preference level of 3.0.
# Acceptor 2 has matched with proposer 2 (since [2,2] is non-zero) at a
# preference level of 3.0.
# Acceptor 3 has matched with proposer 1 (since [3,1] is non-zero) at a
# preference level of 1.0.
#
# INPUT:
# ----------------------------------------------------------------------------------
# P        proposer matrix P.
#          It must be a square matrix with no zeros.
# A        acceptor matrix A.
#          It must be a square matrix with no zeros.
# ordered  If true, P and A are assumed to be ordered,
#          i.e. the leftmost value in a row is the most preferred partner's index.
#          i.e. the leftmost value in a row in P is the preference value for the acceptor with
#          index 1 and vice-versa (higher is better).
# verbose  if the algorithm should print verbosely
# ----------------------------------------------------------------------------------
#
# OUTPUT:
# ----------------------------------------------------------------------------------------
# result_matrix  Result Matrix
# ----------------------------------------------------------------------------------------

m_stableMarriage = function(Matrix[Double] P, Matrix[Double] A, Boolean ordered = TRUE, Boolean verbose = FALSE)
  return (Matrix[Double] result_matrix)
{
  # variable names follow publication

  print("STARTING STABLE MARRIAGE");
  assert(nrow(A) == nrow(P));
  assert(ncol(A) == ncol(P));

  if(nrow(P) != ncol(P) | nrow(A) != ncol(A))
    stop("StableMarriage error: Wrong Input! Both P and A must be square.")

  n = nrow(P)
  # Let S be the identity matrix
  S = diag(matrix(1.0, rows=n, cols=1))
  result_matrix = matrix(0.0, rows=n, cols=n)
  # Pre-processing
  if(!ordered) {
    # If unordered, we need to order P
    ordered_P = matrix(0.0, rows=n, cols=n)
    transposed_P = t(P)

    parfor(i in 1:n)
      ordered_P[,i] = order(target=transposed_P, by=i, decreasing=TRUE, index.return=TRUE)
    P = t(ordered_P)
  }
  else {
    # If ordered, we need to unorder A
    unordered_A = matrix(0.0, rows=n, cols=n)
    # Since cells can be converted to unordered indices independently, we can nest parfor loops.
    parfor(i in 1:n) {
      parfor(j in 1:n, check=0) 
        unordered_A[i, as.scalar(A[i, j])] = n - j + 1
    }
    A = unordered_A
  }

  proposer_pointers = matrix(1.0, rows=n, cols=1)

  while(sum(S) > 0) {
    stripped_preferences = S %*% P
    mask_matrix = matrix(0.0, rows=n, cols=n)
    for(i in 1:n) {
      max_proposal = as.scalar(stripped_preferences[i, as.scalar(proposer_pointers[i])])
      if(max_proposal != 0) {
        proposer_pointers[i] = as.scalar(proposer_pointers[i]) + 1
        mask_matrix[max_proposal, i] = 1
      }
    }
    # make Hadamard Product
    Propose_round_results = mask_matrix * A
    best_proposers_vector = rowIndexMax(Propose_round_results)
    prev_best_proposers = rowIndexMax(result_matrix)

    for(i in 1:n, check=0) {
      new_best_proposer_index = as.scalar(best_proposers_vector[i])
      new_best = as.scalar(Propose_round_results[i, new_best_proposer_index])

      if(new_best > 0) {
        prev_best_proposer_index = as.scalar(prev_best_proposers[i])
        prev_best = as.scalar(result_matrix[i, prev_best_proposer_index])

        if (new_best > prev_best)
        {
          # Proposal is better than current fiance
          result_matrix[i, prev_best_proposer_index] = 0
          result_matrix[i, new_best_proposer_index] = new_best

          #Disable freshly married man to search for a new woman in the next round
          S[new_best_proposer_index, new_best_proposer_index] = 0

          # If a fiance existed, dump him/her
          if(prev_best > 0)
            S[prev_best_proposer_index, prev_best_proposer_index] = 1
        }
      }
    }
  }
  
  if(verbose)
    print("Result: \n"+toString(result_matrix))
}