#OK to post homework
#Yuxuan Yang, April 3rd, Assignment 19
with(combinat):


#1
UmbT:=proc(P,var,N) local n:
[seq(Umb(P^n,var),n=1..N)]:
end:
#UmbT(x + y, [x, y], 20);
#A000984
#UmbT(x+y+z,[x,y,z],20);
#A002893
#UmbT(x1+x2+x3+x4,[x1,x2,x3,x4],20);
#A002895
#UmbT(1+2*x+3*y,[x,y],20);
#Not Found
#UmbT((1+x)*(1+y)*(1+z),[x,y,z],20);
#Not Found


#Maple code for Lecture 19

#C19.txt: March 31, 2021. Maple code for Lecture 19, Using the umbral approach to count Condorcet vote-profile scenarios
Help19:=proc(): print(`Umb1(M,var), Umb(P,var), Zs(N) `): end:
#Umb1(M,var): inputs a MONOMIAL M in the list of variables var and applies the "umbra"
#c*var[1]^a1*...*var[k]^ak->c*(a1+...+ak)!/(a1!*..*ak!)
Umb1:=proc(M,var) local f,d,i,c:
for i from 1 to nops(var) do
d[i]:=degree(M,var[i]):
od:
c:=normal(M/mul(var[i]^d[i],i=1..nops(var))):

if {seq(normal(diff(c,var[i])),i=1..nops(var))}<>{0} then
  RETURN(FAIL):
fi:

c*add(d[i],i=1..nops(var))!/mul(d[i]!,i=1..nops(var)):

end:


#Umb(P,var): applying the above umbra to the polynomial P
#Debugged after class
Umb:=proc(P,var) local P1,i:
if P=0 then
 RETURN(0):
fi:

P1:=expand(P):
if not type(P1,`+`) then
RETURN(Umb1(P,var)):
fi:


add(Umb1(op(i,P1),var),i=1..nops(P1)):
end:

#Zs(N): The first N terms of the sequence: number of vote-profiles with the Condorcet scenario with v voters
Zs:=proc(N) local f,var,i,L:
f:=(1+t^3*x123*x231*x312)*t^3*x312*x123*x231/
((1-t^2*x123*x321)*(1-t^2*x213*x312)*(1-t^2*x312*x123)*(1-t^2*x132*x231)*(1-t^2*x123*x231)*(1-t^2*x231*x312));
f:=taylor(f,t,N+4):
var:=[x123,x132,x213,x231,x312,x321]:

L:=[seq(expand(coeff(f,t,i)),i=1..N)]:

[seq(2*Umb(L[i],var),i=1..nops(L))]:
end:



#C18.txt: Maple code for Lecture 18, original and Generaized Condorcet scenario
Help18:=proc(): print(`IsCond3G(C,G), Cond3G(v,G) , SeqC(N,G), SeqCe(N,G), `):
print(`SeqCo(N,G), NuC(v,G) `):end:

with(combinat):
#i beats j if the number of voters that prefer i to j is larger than
#the number of voters that prefer j to i
#
#S[i,j]:=Number of voters that prefer i to j 
#i beats j if S[i,j]-S[j,i]>=1, 1 is the cut-off
#
#Let G be a 2 by 3 cut-off matrix, [[0,G12,G13],[0,0,G23]]

#IsCond3G(C,G): inputs a vote-count profile outputs true iff 1->2->3->1,
IsCond3G:=proc(C,G) local i,P,T:
P:=permute(3):

for i from 1 to nops(C) do
 T[P[i]]:=C[i]:
od:

#1->2->3->1
if (T[[1,2,3]]+T[[1,3,2]]+T[[3,1,2]])- (T[[2,1,3]]+T[[2,3,1]]+T[[3,2,1]])>=G[1][2] and
(T[[1,2,3]]+T[[2,1,3]]+T[[2,3,1]])-(T[[1,3,2]]+T[[3,1,2]]+T[[3,2,1]])>=G[2][3] and
-(T[[2,3,1]]+T[[3,2,1]]+T[[3,1,2]])+(T[[2,1,3]]+T[[1,2,3]]+T[[1,3,2]])<G[1][3] then
 true:
 else
 false:
fi:
end:


Cond3G:=proc(v,G) local A,a,C:
A:=Comps(v,6):
C:={}:
#C is the SET of vote-count profiles with v voters that lead to 1->2->3->1
for a in A do
if IsCond3G(a,G) then
 C:=C union {a}:
fi:
od:
C:
end:

#Nuc(v,G): The total number of voter-profiles (not to be confused with vote-count-profiles) with 1->2->3->1 with v voters
NuC:=proc(v,G) local A,su,a,i:
A:=Comps(v,6):
su:=0:
for a in A do
 if IsCond3G(a,G) then
  su:=su+add(a)!/mul(a[i]!,i=1..nops(a)):
 fi:
od:
su:
end:





#SeqCo(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with 2*n-1 voters
SeqCo:=proc(N,G) local n:
[ seq(nops(Cond3G(2*n-1,G)),n=1..N)]:
end:

#SeqCe(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with 2*n voters
SeqCe:=proc(N,G) local n:
[ seq(nops(Cond3G(2*n,G)),n=1..N)]:
end:

#SeqC(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with v voters
SeqC:=proc(N,G) local n:
[ seq(nops(Cond3G(n,G)),n=1..N)]:
end:

Help17:=proc(): print(`Comps(n,k), GP1(L,n,d), GP(L,n) , IsCond3(C) `):end:



#Comps(n,k): The set of all compositions of n into exactly k parts (with non-negative entrties)
Comps:=proc(n,k) local S,a,S1,s:
option remember:
if k=0 then
 if n=0 then
   RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:
S:={}:

for a from 0 to n do
 S1:=Comps(n-a,k-1):
 S:=S union {seq([op(s),a],s in S1)}:
 od:
 S:
end:



#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1:=proc(L,n,d) local a,i,P,eq,var:
if nops(L)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,P),i=1..d+2)}:
var:=solve(eq,var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d+3..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:



#GP(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP:=proc(L,n) local d,P:
for d from 0 to nops(L)-2 do
 P:=GP1(L,n,d):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:


#IsCond3(C): inputs a vote-count profile outputs true iff 1->2->3->1,
IsCond3:=proc(C) local i,v,P,T:
P:=permute(3):
v:=add(C):
for i from 1 to nops(C) do
 T[P[i]]:=C[i]:
od:

#1->2->3->1
if T[[1,2,3]]+T[[1,3,2]]+T[[3,1,2]]>v/2 and 
T[[1,2,3]]+T[[2,1,3]]+T[[2,3,1]]>v/2 and
T[[2,3,1]]+T[[3,2,1]]+T[[3,1,2]]>v/2 then
 true:
 else
 false:
fi:
end: