#OK to post homework
#Victoria Chayes, 1/24/2021, Assignment 1

with(combinat): 
Help:=proc(): print(`FP(pi),RandAlmostDer(n,k),EstimateAveFP(n,K), EstimateMomFP(n,r,K),ExactMomFP(n,r)`): end:

#############################################

# 3. Write a procedure FP(pi) that inputs a permutation pi of {1,...n}, (where n=nops(pi)) and outputs an integer between 0 and n indicating the number of fixed points, i.e. the number of places i such that pi[i]=i. Note that if pi is a derangement, then FP(pi)=0.

FP:=proc(pi) local count,i,n:
n:=nops(pi):
count:=0:
for i from 1 to n do
  if pi[i]=i then
    count:=count+1:
  fi:
od:
count:
end: 

# 4. Modify RandDer(n) to write a procedure RandAlmostDer(n,k) that inputs a positive integer n, another positive integer k, (between 0 and n) and outputs a random permutation whose number of fixed points is <=k. Note that RandAlmostder(n,0) does the same thing as RandDer(n).

RandAlmostDer:=proc(n,k) local pi:
pi:=randperm(n):
while  FP(pi)>k do
  pi:=randperm(pi):
od:
pi:
end:

# 5. Write a procedure EstimateAveFP(n,K) that inputs a positive integer n, and a large positive integer K, picks randompermutations of {1,...,n} (using randperm(n) in the combinat package of Maple) and outputs the average number of fixed points. By taking K fairly large, say K=10000, and various n, can you guess the exact value?

EstimateAveFP:=proc(n,K) local pi, i, count:
count:=0.0:               #this usually would be 0; I chose 0.0 because it makes count/K a decimal
for i from 1 to K do
  pi:=randperm(n):
  count:=count+FP(pi):
od:
count/K: 
end:

#I ran the following tests:
#EstimateAveFP(2,10000)
#                          1.004400000
#
#EstimateAveFP(5,10000)
#                          0.9895000000
#
#EstimateAveFP(10,10000)
#                          0.9926000000
#
#EstimateAveFP(15,10000)
#                          0.9981000000
#
#EstimateAveFP(20,10000)
#                          1.008700000
#
#My prediction is therefore independent of n, the average number of fixed points is 1.


# 6. Can you find the exact value of the average? Can you prove it mathematically?
# The average number of fixed points is exactly 1: for a random permutation, each entry has probability 1/n of being a fixed point. There are n entries, and sum_{i=1}^n 1/n =1. Alternatively, the program ExactMomFP(n,r) that I wrote for Problem 8 outputs ExactMomFP(n,1)=1. 


# 7. Write a procedure EstimateMomFP(n,r,K) that inputs a positive integer n, another positive integer, r, and a large positive integer K, picks randompermutations of {1,...,n} (using randperm(n) in the combinat package of Maple) and outputs the average number of FP(pi)^r. Note that EstimateMomFP(n,1,K) does the same thing as EstimateAve(n,K). By taking K fairly large, say K=10000, r=2, and various n, can you guess the exact value of the average of FF(pi)^2?

EstimateMomFP:=proc(n,r,K) local pi, i, c:
c:=0.0:
for i from 1 to K do
  pi:=randperm(n):
  c:=c+FP(pi)^r:
od:
c/K: 
end:

#I ran the following tests:
#EstimateMomFP(2,2,10000)
#                          1.941600000
#
#EstimateMomFP(5,2,10000)
#                          2.032800000
#
#EstimateMomFP(10,2,10000)
#                          2.009600000
#
#EstimateMomFP(15,2,10000)
#                          1.976000000
#
#EstimateMomFP(20,2,10000)
#                          1.973800000
#
#EstimateMomFP(100,2,10000)
#                          2.029900000
#
#EstimateMomFP(200,2,10000)
#                          1.963100000
#
#My prediction is therefore independent of n, the average of FP(pi)^2 is 2.


# 8. Can you find the exact value of the average of random variable FP(pi)^2 over all permutations of length n? Can you prove it mathematically?
#
# My instinct as this is experimental math is to first run an experiment and look for a pattern. It is my guess that FP(pi)^k is independent of n at least for large enough n, and is going to be some sort of sequence, and if I can identify that sequence I will better know what I'm looking for. Running 
#
#seq(EstimateMomFP(2,i,10000),i=1..10)
#     1.006800000, 2.013600000, 3.938400000, 7.971200000, 15.70880000, 31.62240000, 64.48640000, 127.9744000, 254.5664000, 510.8736000
#seq(EstimateMomFP(5,i,10000),i=1..10)
#     1.010200000, 1.986600000, 5.139700000, 15.24690000, 49.89550000, 206.4944000, 813.0305000, 3558.047500, 16578.16890, 87267.79750
#
#seq(EstimateMomFP(10,i,10000),i=1..10)
#     1.001300000, 1.932400000, 5.298300000, 14.95400000, 56.34840000, 200.0395000, 849.8835000, 3070.109000, 17112.06940, 94387.61040
#
#seq(EstimateMomFP(20,i,10000),i=1..10)
#     0.9964000000, 2.020000000, 5.122000000, 14.56050000, 57.86080000, 183.3107000, 1090.158100, 4657.022300, 25063.90540, 2.176568281 10 
#
#seems to confirm my hypothesis, and potentially reminds me of Bell's numbers? However, returning to the problem at hand:
#
#The mathematically rigorous way to do this is to come up with and then evaluate the expression of:
# sum_{i=1}^n i^2*[probability there are i fixed points]
#For small n, this is easy to do: ie when n=2, we have [2,1]=0 fixed points and [1,2]=2 fixed points. then 2^2*(1/2)=2. For n=3, we have 6 possible configurations, two with no fixed point, three with one fixed point, and one with three fixed points, so 0^2*(1/3)+1^2*(1/2)+3^2*(1/6)=2. 
#
#To test that this pattern continues, I create the following extremely naive program:

#ExactMomFP(n,r)
#INPUTS: 
#  an integer n, indicating the length of the list of permutations 
#  a number r, for the power FP(pi)^r whose average you want calculated 
#OUTPUTS:
#  The exact value of the average of random variable FP(pi)^r over all permutations of length n.
ExactMomFP:=proc(n,r) local S,P,i,k:
P:=[seq(0,i=1..(n+1))]:
S:=permute(n):
for i from 1 to nops(S) do
 k:=FP(S[i])+1:
 P[k]:=P[k]+1:
od:
add(P[i+1]*i^r,i=0..n)/nops(S):
end:

#Running the tests 
#seq(ExactMomFP(n,2),n=2..9)
#                 2, 2, 2, 2, 2, 2, 2, 2
#confirms that FP(pi)^2=2 at least for sequences up to length 9. As the above program was not written whatsoever for efficiency, testing for n>9 took too long and I terminated the run. Of course, I'm sure there's an easier way to figure out the pattern for probability of there being i fixed points to which we would simply calculate the final line of the program, but I'm lazy and it was easier to make the computer count.
#
#Out of curiosity to see if the above EstimateMomFP estimates for r=3, r=4, etc were correct, I ran the  tests 
#seq(ExactMomFP(n,3),n=2..9)
#                 4, 5, 5, 5, 5, 5, 5, 5
#seq(ExactMomFP(n,4),n=2..9)
#                 8, 14, 15, 15, 15, 15, 15, 15
#seq(ExactMomFP(n,5),n=2..9)
#                 16, 41, 51, 52, 52, 52, 52, 52
#seq(ExactMomFP(n,6),n=2..9)
#                 32, 122, 187, 202, 203, 203, 203, 203
#seq(ExactMomFP(n,7),n=2..9)
#                 64, 365, 715, 855, 876, 877, 877, 877
#
#So for large enough sequences, we seem to recover the pattern of 2,5,15,52,203,877,... confirming the hypothesis that we get Bell's numbers. 

#A comment on the efficiency of ExactMomFP: I actually tried to write a more efficient version, assuming that perhaps storing the list of permutations was taking a lot of time, to try to take advantage of Maple's firstperm and nextperm commands and not have to calculate the permutations myself. However, the code 
#
#ExactMomFP:=proc(n,r) local m,pi,P,i,k,temp:
#P:=[seq(0,i=1..(n+1))]:
#pi:=firstperm(n):
#m:=numbperm(n):
#for i from 1 to m-1 do
# k:=FP(pi)+1:
# P[k]:=P[k]+1:
# pi:=nextperm(pi):
#od:
#k:=FP(pi)+1:
#P[k]:=P[k]+1:
#add(P[i+1]*i^r,i=0..n)/m:
#end:
#
#resulted in runtime
#time(seq(ExactMomFP(n,2),n=2..9))
#                             5.763
#when the original resulted in the runtime 
#time(seq(ExactMomFP(n,2),n=2..9))
#                             3.891
#and hence I have chosen to keep the original rather than get fancy.