#OK to post homework
#Blair Seidler, 2021-02-14, Assignment 7

with(combinat): with(plots):

Help:=proc():
print(`NiceValentine()`): 
print(`DrawMandelbrot(R,K)`):
print(`RandSPn(n)`):
print(`NuAdjacentMates(SP,n)`,`EstimateAvAdjacentMates(n,K)`,`CalcAvAdjacentMates(n)`):
end:


# 1.


MyVal:=proc():
plot([16*sin(t)^3, 13*cos(t) - 5*cos(2*t) - 2*cos(3*t) - cos(4*t), t = 0 .. 2*Pi], axes = none, color = red)
end:

NiceValentine:=proc():
implicitplot3d((x^2+9*y^2/4+z^2-1)^3-x^2*z^3-9*y^2*z^3/200=0, x=-2..2, y=-2..2, z=-2..2, grid=[100,100,100],
style=surface,axes=none,color=red,title="I Love Experimental Mathematics",
caption="Happy Valentine's Day, Dr.Z!",scaling=constrained):
end:


# 2.

#DrawMandelbrot(R,K): plots the Mandelbrot set by picking all those complex c
#(of the form x+I*y, with x,y in [-4,4] with resolution 1/R, and deciding whether
#to accept it if iterating z goes to z^2+c (starting at z=0) after K iterations stay bounded.
#Also send the output as .jpeg file, for R as large as possible.
DrawMandelbrot:=proc(R,K) local th,f,xmin,xmax,ymin,ymax,res,pts,x,y,z,c,i:
xmin:=-4: xmax:=4: ymin:=-4: ymax:=4: th:=10^3:
f:=(z,c)->evalf(evalc(z^2+c)):
res:=1/R:

pts:=[]:
x:=xmin:
while x+res<=xmax do
  y:=ymin:
  while y+res<=ymax do
    c:=x+I*y:
    z:=0:
    for i from 1 to K do
      z:=f(z,c):
    od:
    if abs(z)<th then
      pts:=[op(pts),[x,y]]:
    fi:
    y:=y+res:
  od:
  x:=x+res:
od:

pointplot(pts):

end:

# 3.

#RandSPn(n): inputs a positive integer n and outputs a (uniformly) at random set
#partition of {1,...,n} using LD(L) and RandSPnk(n,k)
RandSPn:=proc(n) local c,k,L:
L:=[]:
for k from 1 to n do
  L:=[op(L),NuSPnk(n,k)]:
od:
c:=LD(L):
RandSPnk(n,c):
end:

# 4.

#NuAdjacentMates(SP,n): 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,c,s:
c:=0:
for i from 1 to n-1 do
  if not evalb({seq(evalb({i,i+1} subset s),s in SP)}={false}) then
    c:=c+1:
  fi:
od:
c:
end:

#EstimateAvAdjacentMates(n,K): uses the above procedure, and RandSPn(n), K times, to estimate
#the average of the number of adjacent mates in a random set partition of {1,...n}.
EstimateAvAdjacentMates:=proc(n,K) local i,c:
c:=0:
for i from 1 to K do
  c:=c+NuAdjacentMates(RandSPn(n),n):
od:
c/K:

end:

#What did you get for EstimateAvAdjacentMates(100,1000)?
#(run it several times to make sure that they are close)
# Several runs at (100,1000) and one at (100,10000) seem to indicate that 3.4 is the average


#CalcAvAdjacentMates(n): uses the above procedure, and RandSPn(n), K times, to estimate
#the average of the number of adjacent mates in a random set partition of {1,...n}.
CalcAvAdjacentMates:=proc(n) local i,c,k,SP:
c:=0:
for k from 1 to n-2 do
  c:=c+(1+k)*NuSPnk(n-2,k):
od:

c*(n-1)/NuSPn(n):

end:

#### Running this procedure, I get a value of 3.367670044 for the average number of adjacent
#### mates when n=100.


#### From C7.txt ####
#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:

#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:


#### From C4.txt ####
#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:
#END FROM C4.txt