#OK to post homework
#Wanying Rao, 02/03/2021, Assignment 6

with(plots):

#1
NiceValentine := proc() local p,t:
p := plot([(1 + cos(t))*sin(t), -(1 + cos(t))*cos(t), t = 0 .. 2*Pi], axes=none, thickness=5,colorscheme=["ygradient",["Red","Pink"]]):
t := textplot([0,-0.8,"I love Experimental Mathematics"]):
display({p,t}):
end:

#2
DrawMandelbrot := proc(R,K) local z,i,c,x,y,s:
s := {}:
y := 4:
while y>=-4 do
x:=4:
while x>=-4 do
 c := x+I*y:
 z := 0:
 for i from 1 to K do
  z := evalf(evalc(z^2+c))
 od:
  if abs(z)<10 then
   s := s union {[x,y]}:
  fi:
 x := x-1/R:
od:
y := y-1/R
od:
pointplot(s):
end:

#3
#By using LD(L) and RandSPnk(n,k) write a procedure
#RandSPn(n)
#that inputs a positive integer n and outputs a (uniformly) at random set partition of #{1,...,n}
RandSPn := proc(n) local i,L,r:
L := [seq(NuSPnk(n,i),i=1..n)]:

r := LD(L):
RandSPnk(n,r):
end:


#4
#(i)Write a procedure NuAdjacentMates(SP,n)
#that inputs a set partition, SP, of {1, ...n}, and outputs the number of i between 1 and #n such that i and i+1 are part-mates

NuAdjacentMates := proc(SP,n) local i,s,m:
m := 0:

for i from 1 to n do
 for s in SP do
  if {i, i+1} subset s then
   m := m+1:
  fi:
 od:
od:

m:
end:

#(ii)Write a procedure EstimateAvAdjacentMates(n,K)
#that uses the above procedure, and your RandSPn(n), K times, to estimate the average of #the number of adjacent mates in a random set partition of {1,...n}. What did you get for
#EstimateAvAdjacentMates(100,1000)?

EstimateAvAdjacentMates := proc(n,K) local i,s,a:
a := 0:
for i from 1 to K do
 a := a + NuAdjacentMates(RandSPn(n),n):
od:

a/K:

end:


#I got 1749/500 for EstimateAvAdjcentMates(100,1000). (Also 67/20,167/50)









#From C7.txt
Help7:=proc(): print(` SPnk(n,k), SPn(n) , NuSPnk(n,k), NuSPn(n)`):
print(`RandSPnk(n,k)`):
end:

#SPnk(n,k): Inputs pos. integers k and n and 1<=k<=n
#outputs the SET of set-partitions of {1,...,n} with exactly
#k parts
SPnk:=proc(n,k) local SP,SP1,sp,i:
option remember:
if n=1 then
  if k=1 then 
   RETURN({{{1}}}):
  else
    RETURN({}):
  fi:
fi:

#Case (i): n is a singleton
SP1:=SPnk(n-1,k-1):

SP:={seq(sp union {{n}},sp in SP1)}:

#Case (ii) n has friends
SP1:=SPnk(n-1,k):

for sp in SP1 do
 for i from 1 to nops(sp) do
  SP:=SP union {{op(1..i-1,sp), sp[i] union {n},op(i+1..nops(sp),sp)}}:
 od:
 od: 
SP:
end:  

#SPn(n): The set of all set-partitions of {1,...,n}
#(same as combinat[setpartition](n))
SPn:=proc(n) local k:  { seq(  op(SPnk(n,k)),k=1..n)}:end:

#NuSPnk(n,k): The number of set-partitions of {1,...,n}
#into exactly k parts (aka Strirling numbers of the second kind)
NuSPnk:=proc(n,k) option remember:

if n=1 then
  if k=1 then
   RETURN(1):
  else
   RETURN(0):
 fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end:

NuSPn:=proc(n) local k:
option remember:
add(NuSPnk(n,k),k=1..n):
end:

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:

#RandSPnk(n,k): a random set partition of {1,..,n} with k parts
RandSPnk:=proc(n,k) local c,SP1,k1:
if n=1 then
 if k=1 then
   RETURN({{1}}):
  else
   RETURN(FAIL):
 fi:
fi:
c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]):

if c=1 then
 SP1:=RandSPnk(n-1,k-1):
 RETURN(SP1 union {{n}}):
 else
 SP1:=RandSPnk(n-1,k):
  k1:=rand(1..k)():
  RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}):
  
fi:
end: