# OK to post homework
# Vikrant, Feb 27 2022, Assignment 11

# ================================================================================
# 0. Code that has been given.
# ================================================================================

read("C11.txt"):

# ================================================================================
# 1. Validate GS inputs.
# ================================================================================

IsPerm:=
    proc(pi)
        local i:
        local S:= {seq(i,i=1..nops(pi))}:
        for i in pi do
            S:= S minus {pi[i]}:
        od:
        evalb(S = {}):
    end:

GaleShVal:=
    proc(M,W)
        local X,i:
        for X in [M,W] do
            if not type(X,list) then
                print("M,W must be a list."):
                return(false):
            fi:
            if nops(X) <> nops(M) then
                print("M,W must have the same length."):
                return(false):
            fi:
            for i from 1 to nops(M) do
                if not type(X[i],list) then
                    print("Entries of M,W must be a list."):
                    return(false):
                fi:
                if nops(X[i]) <> nops(M) then
                    print("Entries of M,W must have the same length as W,M."):
                    return(false):
                fi:
                if not IsPerm(X[i]) then
                    print("Entries of M,W must be permutations."):
                    return(false):
                fi:
            od:
        od:
        true:
    end:

GaleShV:=
    proc(M,W)
        if not GaleShVal(M,W) then
            return(FAIL):
        fi:
        GaleSh(M,W):
    end:

GaleShTV:=
    proc(M,W)
        if not GaleShVal(M,W) then
            return(FAIL):
        fi:
        GaleShT(M,W):
    end:

# ================================================================================
# 2. Women proposers.
# ================================================================================

GaleShTW:=
    proc(M,W)
        local O:= GaleShTV(W,M):
        inv(O[1]),O[2]:
    end:

# ================================================================================
# 3. Check proposer's better off.
# ================================================================================

CheckOpt:=
    proc(M,W)
        local i:
        local b:= true:
        local MProposer:= GaleShTV(M,W)[1]:
        local WProposer:= GaleShTW(M,W)[1]:
        for i from 1 to nops(M) do
            b:= b and (M[i][MProposer[i]] <= M[i][WProposer[i]]):
        od:
        b:
    end:

# ================================================================================
# 4. GS with partial preferences.
# ================================================================================

GetMissing:=
    proc(ppi,n)
        local i:
        {seq(i,i=1..n)} minus {op(ppi)}:
    end:

# The problem statement says M[i] and W[j] are not always the same length.
# This suggests the meaning of M[i],W[j] must have changed: It is a list of persons indexed by ranks rather than ranks indexed by persons.
# We will thus need to invert them before feeding them into GaleShT, which would (un)invert them in its implementation.
GaleShP:=
    proc(M,W)
        local i:
        local n:= nops(M):
        local M2:= [seq(inv([op(M[i]),n+1,op(GetMissing(M[i],n))]),i=1..n),[seq(i,i=1..n+1)]]:
        local W2:= [seq(inv([op(W[i]),n+1,op(GetMissing(W[i],n))]),i=1..n),[seq(i,i=1..n+1)]]:
        local Match:= GaleShTV(M2,W2):
        if Match[1][n+1] <> n+1 then
            return(FAIL):
        fi:
        [op(1..n,Match[1])],Match[2]:
    end:

# ================================================================================
# 5. GS vs Naive.
# ================================================================================

(*

FailsTooMuch:=
    proc(n,K)
        local i:
        evalb(nops([seq(`if`(FindStable(RT(n),10^4)<>FAIL,1,NULL),i=1..K)])/K < 1/2):
    end:

randomize(1937681705):
l:= 1:
h:= 30:
while l < h do
    m:= iquo(l+h,2):
    if not FailsTooMuch(m,100) then
        l:= m+1:
    else
        h:= m:
    fi:
od:
# 20 seems to be the threshold. Out of 100 random instances, less than 50% succeed.
l;

*)

(*

randomize(1317627810):
K:= 30:
for n from 1 to 19 do
    gsco:= 0:
    naco:= 0:
    gswins:= 0:
    for i from 1 to K do
        P:= RT(n):
        gs:= GaleShTV(P):
        na:= FindStable(P,10^4):
        gsco:= gsco+gs[2]:
        naco:= naco+10^4:
        if na<>FAIL then
            naco:= naco-10^4+na[2]:
        fi:
        if na=FAIL or gs[2]<=na[2] then
            gswins:= gswins+1:
        fi:
    od:
    print(gsco,naco,gswins/K);
od:

# The Naive local optimization does fare better on occasion.

*)

# ================================================================================
# 6. Estimate number of proposals.
# ================================================================================

EstStat:=
    proc(n,K)
        local i:
        local d:= diff(add([seq(x^(GaleShTV(RT(n))[2]),i=1..K)]),x)/K:
        evalf(subs(x=1,[d,sqrt(diff(d,x)+d-d^2)])):
    end:

(*

randomize(503735861):
for n from 10 by 10 to 150 do
    print(EstStat(n,10^4));
od:

# [23.96310000, 6.100273632]
# [62.45690000, 16.06703278]
# [106.8092000, 26.63632098]
# [156.1588000, 38.05825250]
# [209.1352000, 51.44857161]
# [263.3541000, 64.02588940]
# [319.7320000, 76.16492484]
# [376.2190000, 86.34193558]
# [436.8146000, 100.7316417]
# [494.9200000, 110.8625239]
# [557.2744000, 125.7744477]
# [620.4146000, 138.1635245]
# [685.9910000, 150.3574119]
# [748.2796000, 161.1688581]
# [810.6049000, 175.1450113]

# [Theta(n ln n), Theta(n)]?
# Coupon-collector-like (I've seen this comparison before).

*)

# ================================================================================
# 7. Maybe a not so stupid stable matching algorithm.
# ================================================================================

PNeighbors:=
    proc(pi)
        local pi2,i,t:
        local S:= {}:
        for i from 1 to nops(pi)-1 do
            pi2:= pi:
            t:= pi2[i]:
            pi2[i]:= pi2[i+1]:
            pi2[i+1]:= t:
            S:= S union {pi2}:
        od:
        S:
    end:

MaybeNotSoStupidStable:=
    proc(M,W)
        MNSSS(M,W,randperm(nops(M)),0):
    end:

MNSSS:=
    proc(M,W,pi,co)
        local s,t:
        local n1:= CountInstabilities(M,W,pi):
        local S:= PNeighbors(pi):
        local pi2:
        local n2:= nops(M)^3+1:
        for s in S do
            t:= CountInstabilities(M,W,s):
            if t<n2 then
                pi2:= s:
                n2:= t:
            fi:
        od:
        if n2<n1 then
            return(MNSSS(M,W,pi2,co+1)):
        fi:
        if n1>0 then
            return FAIL:
        fi:
        pi,co:
    end:

# Seems like it's kinda stupid.

# ================================================================================
# 8. Read [Shapley], [Hausken, Mohr].
# ================================================================================

# Done.

# ================================================================================
# oo. Helpers from past homeworks.
# ================================================================================

CountInstabilities:=
    proc(M,W,pi)
        local i1,i2:
        add([seq(seq(`if`(IsStable1(M,W,pi,i1,i2),0,1),i1=1..nops(pi)),i2=1..nops(pi))]):
    end: