# Ok to post homework
# Lucy Martinez, 03-06-22, Assignment 13
read `C13.txt`:
#----------------------Problem 3-----------------------#
# Implement the more general formula given in the the wikipedia article to get what they call ÆC(v),
# i.e. the Generalized Shapley value. Call it SVC(f,n,C). 
# Make sure that SVC(f,n,{i}) is the same as SVi(f,n,i). Use the first formula given there (analogous to SVi)

SVC:=proc(f,n,C) local N,i,diff,PSd,PSc,T,S:
N:={seq(i,i=1..n)}:
 
diff:=N minus C:
PSd:=powerset(diff):
PSc:=powerset(C):
 
add( add( (((n-nops(T)-nops(C))!*(nops(T))!)/((n-nops(C)+1)!))*(((-1)^(nops(C)-nops(S)))*(f[S union T])) , S in PSc),T in PSd ):
end:

# Checking the previous procedure:
f:=RG(6,10):
SVi(f,6,5);
# 11/6

SVC(f,6,{5});
# 11/6

f1:=RG(4,10):
SVi(f1,4,2);
# 1/3

SVC(f1,4,{2});
# 1/3


#----------------------Problem 4-----------------------#
# Same as above, but using the second formula there. 
# Call it SVC1(f,n,C). Make sure that SVC1(f,n,{i}) is the same as SVi(f,n,i). It is analogous to SV(f,n)
# Hint: what the article calls "synnergy" is our own BM(f,n). What they call v(S) is our f[S]. 
#       What they call Æi(v) is our SVi(f,n,i).

SVC1:=proc(f,n,C) local N,i, PSn,sum, T:
N:={seq(i,i=1..n)}:
PSn:=powerset(N) minus {{}}:
sum:=0:
 
for T in PSn do
 if C intersect T = C then
  sum:=sum+(BM(f,n)[T])/(nops(T)-nops(C)+1): 
 fi:
od:
sum:
end:

# Checking SVC1(f,n,C) works:
f1:=RG(5,7):
SVC1(f1,5,{5});
# 33/20

SVi(f1,5,5);
# 33/20


#----------------------Problem 5-----------------------#

# Implement the formulas for Æi,j(v) Given there, using both formulas, call them 
# SVij1(n,f,i,j) and SVij2(n,f,i,j). Make sure that add(SVij1(n,f,i,j),j=1..n)=SVi(f,n,i),
# add(SVij2(n,f,i,j),j=1..n)=SVi(f,n,i), for randomly chosen functions f, gotten from RSG(8,1000), 
# and all i from 1 to 8.

SVij1:=proc(n,f,i,j) local N, k, PSn, S,t :
N:={seq(k,k=1..n)}:
PSn:=powerset(N) minus {{}}:
 
add( add( ((nops(S)-1)!*(n-nops(S))!)*(f[S]-f[S minus {i}]-f[S minus {j}]+f[S minus {i,j}])*(1/t),t=nops(S)..n), S in PSn )/n!:
end:


SVij2:=proc(n,f,i,j) local N, k, PSn, sum,S, T:
N:={seq(k,k=1..n)}:
PSn:=powerset(N) minus {{}}:
sum:=0:

for S in PSn do
 if {i,j} intersect S ={i,j} then
  sum:=sum+ BM(f,n)[S]/(nops(S))^2: 
 fi:
od:
end:

# Checking SVij1(n,f,i,j) works:
f2:=RSG(8,1000):
for i from 1 to 8 do
 add(SVij1(6,f2,i,j),j=1..6),SVi(f2,6,i);
od;
#             90604/15,  90604/15                  
#            163237/30,  163237/30
#              15341/3,  15341/3
#              22093/4,  22093/4
#            105873/20,  105873/20
#             87149/15,  87149/15
#                    0, 0
#                    0, 0

# Checking SVij2(n,f,i,j) works:
f3:=RSG(8,1000):
for i from 1 to 8 do
 add(SVij2(6,f3,3,j),j=1..6),SVi(f3,6,3);
od;

#         15067/3,  15067/3
#       359059/60,  359059/60
#       105449/20,  105449/20
#       290327/60,  290327/60
#        71603/12,  71603/12
#         24046/5,  24046/5
#               0,  0
#               0,  0