# OK to post homework
# RDB, 2022-03-06, HW13

# 1.

# Let's call SVpF the "subset" approach and SVp the "permutation" approach.

# The permutation approach is the one we discussed in class: Permute the
# players and give each player the "value of the room" once they appear minus
# "the value of the room" without them. Do this for all permutations and divide
# by n!.

# The subset approach is the same, but traverses the permutations in a
# different way. Fix a player i, then pick a set of players S which contains i.
# The set S \ {i} is what comes "before" i. No matter what order those players
# enter, the value of the room before i enters will be the same, and the value
# of the room after i enters will be the same. The order of the players *after*
# i is also irrelevant. If there are n players, then there are (|S| - 1)! * (n
# - |S|)! ways to arrange the players in this setup for a fixed S, which
# assigns i a value of f(S) - f(S \ {i}). The total sum for player i will be
# the same as above.

# In the permutation approach, we look over n! permutations and do some
# constant-time operations, so the runtime is Theta(n!).

# In the subset approach, we look over at most 2^n sets for each of the n
# players and do some constant-time operations, so the runtime is Theta(n *
# 2^n).

# n * 2^n grows *much* slower than n!, so the subset approach is going to win
# out for even moderately small n. (But both will become infeasible before too
# far.)

# 2.

# The axiomatic approach leads to this:

# Let f(S) be our "value" function. Write f(S) = sum(g(T), T <= S) for some g.
# We can determine g via Mobius inversion. Once we do, we assign the atomic
# game corresponding to S the value g(S) / |S|!. The value for player i is then

#       V(i) = sum(g(S) / |S|, S <= [n], i in S).

# We want to show that this equals

#       R(i) = sum((|S| - 1)! (n - |S|)! / n! * (f(S) - f(S \ {i})), S <= [n]).

# Let's toil away on V(i).

#       V(i) = sum(sum(f(T) (-1)^(|S| - |T|) / |S|, T <= S), S <= [n], i in S)
#            = sum(f(T) (-1)^|T| sum((-1)^(|S|) / |S|, S >= T, i in S), T <= [n]).

# We can break this into two sums: one where i is in T (call it A(i)), and one
# where it is not (call it B(i)).

# Then:

#   A(i) = sum(f(T) (-1)^|T| sum((-1)^(|S|) / |S|, S >= T, i in S), T <= [n], i in T).
#        = sum(f(T) (-1)^|T| sum((-1)^(s + |T|) binomial(n - |T|, s) / (s + |T|), s >= 0), T <= [n], i in T)
#        = sum(f(T) sum((-1)^s binomial(n - |T|, s) / (s + |T|), s >= 0), T <= [n], i in T)

#   B(i) = sum(f(T) (-1)^|T| sum((-1)^(|S|) / |S|, S >= T, i in S), T <= [n], i not in T).
#        = sum(f(T) (-1)^|T| sum((-1)^(s + |T|) binomial(n - |T| - 1, s - 1) / (s + |T|), s >= 1), T <= [n], i not in T)
#        = sum(f(T) sum((-1)^(s + 1) binomial(n - |T| - 1, s) / (s + |T| + 1), s >= 0), T <= [n], i not in T).
#        = -sum(f(T) sum((-1)^s binomial(n - |T| - 1, s) / (s + |T| + 1), s >= 0), T <= [n], i not in T).

# We want the two outer sums to have the same range. Fortunately,
# {T: i not in T} = {T \ {i}: i in T}, Thus:

#   B(i) = -sum(f(T \ {i}) sum((-1)^s binomial(n - |T|, s) / (s + |T|), s >= 0), T <= [n], i in T).

# It follows that

#   V(i) = A(i) + B(i)
#        = sum((f(T) - f(T \ {i})) sum((-1)^s binomial(n - |T|, s) / (s + |T|), s >= 0), T <= [n], i in T).

# CLAIM:

#       sum((-1)^s binomial(n - |T|, s) / (s + |T|), s >= 0) = (|T| - 1)! (n - |T|)! / n! for 1 <= |T| <= n.

# PROOF:

#   Let k = |T|, and F(n, k, s) the summand after dividing through by the
#   right-hand side. That is,

(*
                                s
                            (-1)  binomial(n - k, s) n!
                       F := ---------------------------
                             (s + k) (k - 1)! (n - k)!

*)

#   The sum we want to evaluate is Z(n, k) = sum(F(n, k, s), s=0..n-k).

#   It is routine to verify that

#     (1)          F(n + 1, k, s) - F(n, k, s) = G(n, k, s + 1) - G(n, k, s),

#   where

#           G(n, k, s) = (s(s + k) / (k + s - 1 - n) / (n + 1 - k)) * F(n, k, s).

# If we sum Eq. (1) from s = 0 to s = n - k, we get

#       Z(n + 1, k) - Z(n, k) = sum(G(n, k, s + 1) - G(n, k, s), s=0..n-k)
#                             = G(n, k, n - k + 1) - G(n, k, 0)
#                             = 0.

# This implies Z(n, k) = Z(k, k) for all n >= k:

#   Z(n, k) = Z(k, k) = sum(F(k, k, s), s=0..0)
#                     = F(k, k, 0)
#                     = 1.

# Therefore Z(n, k) = 1 for all n >= k, which implies the identity we want.

# Q.E.D.

# (There's definitely a better way to do this.)

# 3.

read "C13.txt":

with(combinat):

SVC := proc(f, n, coalition)
    g := BM(f, n):

    antiC := {seq(1..n)} minus coalition:
    sets := {seq(coalition union S, S in powerset(antiC))}:
    add(g[T] / (nops(T) - nops(coalition) + 1), T in sets):
end:

# Here's a test.
f := RSG(4, 10):
SVi(f, 4, 1);
SVC(f, 4, {1});

# 4.

# Thanks to Natasha for pointing out the "divide by 0" error in the following
# Wikipedia formula:

#   phi_{ij}(v) = sum((|S| - 1)! (n - |S|)! / n! * (v(S) - v(S \ {i}) - v(S \ {j}) + v(S \ {i, j})) * sum(1 / t, t=|S|..n), S <= [n]).

# Natasha and I checked the original paper and determined that the outer sum
# ought to have the condition "i, j in S," which fixes the divide by zero
# problem. (The authors are probably using a convention which sets the summand
# to zero whenever i and j are not in S.)

SVij1 := proc(f, n, i, j)
    H := S -> add(1 / k, k=nops(S)..n):
    V := S -> f[S] - f[S minus {i}] - f[S minus {j}] + f[S minus {i, j}]:
    add(ifelse({i, j} subset S, (nops(S) - 1)! * (n - nops(S))! * V(S) * H(S), 0), S in powerset(n)) / n!:
end:

# Here's code to test.
# Looks good!
print(`TESTING SVij1`);
n0 := 8:
f := RSG(n0, 1000):
for i from 1 to n0 do
    add(SVij1(f, n0, i, j), j=1..n0) = SVi(f, n0, i);
od;

# 5.

# The formula which comes *after* the previous formula, namely

#   phi_{ij}(v) = sum(w(S) / |S|^2, {i, j} <= S <= [n])     (w = BM(v)),

# seems to have the correct restriction on S. I didn't figure out the correct
# manipulations to prove that this is equivalent to the previous one (involves
# *four* sums!), but it seems to work.

SVij2 := proc(f, n, i, j)
    g := BM(f, n):

    add(ifelse({i, j} subset S, g[S] / nops(S)^2, 0), S in powerset(n)):
end:

print(`TESTING SVij2`);
for i from 1 to n0 do
    add(SVij2(f, n0, i, j), j=1..n0) = SVi(f, n0, i);
od;