#OK to post homework
#AJ Bu, February 11 2021, Assignment 7

with(plots):

######PROBLEM 1######
NiceValentine:=proc():
implicitplot3d((x^2 + 9/4*y^2 + z^2 - 1)^3 - x^2*z^3 - 9/40*y^2*z^3, x = -1.2 .. 1.2, y = -1 .. 1, 
z = -1.1 .. 1.3, numpoints = 100000, color = red, shading = zhue, glossiness = 1, style = surface, 
axes = none, title = "I Love Experimental Math", orientation = [110, 90, 40]):
end:



#######PROBLEM 2#######


DrawMandelbrot:=proc(R,K) local x,y,z,i,j,c,S,i1:
	S:={}:

	for i from -4*R to 4*R do
			x:=i/R: 
		for j from -4*R to 4*R do
			y:=j/R:
			c:=x+y*I: 			
			z:=0:
		
			for i1 from 1 to K do
				z:=evalf(evalc(z^2+c)): 
			od:
	
			if evalf(abs(z))<10 then 
				S:=S union{[x,y]}:
			fi:
			
		od:
	od:
	
	plot(S,style=point):
end:



######PROBLEM 3######
RandSPn:=proc(n) local Die1,k:
	if n=1 then 
		RETURN({{1}}): 
	fi:

	Die1:=[seq(NuSPnk(n,i),i=1..n)]:
	
	k:=LD(Die1):

	RandSPnk(n,k):
end:



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

NuAdjacentMates:=proc(SP,n) local i,j,k,S,s:
	k:=nops(SP):
	i:=0:
	
	for j from 1 to k do
		for s in SP[j] do
			if s+1 in SP[j] then
				i:=i+1:
			fi:
		od:
	od:

	i:
end:



EstimateAvAdjacentMates:=proc(n,K) local m,SP,i:
	m:=0:
	for i from 1 to K do 
		SP:=RandSPn(n):
		m:=m+NuAdjacentMates(SP,n): 
	od:
	m/K:
end:

#Running EstimateAvAdjacentMates(100,1000) multiple times, I 
#consistently got a number between 3.3 and 3.4.
#Some example outputs for evalf(EstimateAvAdjacentMates(100,1000)) were
##		3.468000000
##		3.361000000
##		3.328000000
##		3.416000000
##		3.395000000










######################  FROM C7.txt  ######################
#SPnk(n,k): Inputs pos. integers k and n 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): if 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 Stirling number 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: