Project Euler #645 -- speed up Monte-Carlo simulation in Python

Question

I am trying to solve Q645. While the logic used for my code seems to be appropriate, the code itself is way too slow for the large number required in this question. May I ask for suggestion(s) to improve my code's performance?

The question is as in the link: https://projecteuler.net/problem=645

My Python code is as follow:

def Exp(D):
    day_list = [0]*D
    num_emperor = 0
    while all((d == 1 for d in day_list)) == False:
        #the birthday of the emperors are independent and uniformly distributed throughout the D days of the year
        bday = np.random.randint(0,D)
        day_list[bday] = 1
        num_emperor+=1
        #indices of d in day_list where d == 0
        zero_ind = (i for i,v in enumerate(day_list) if v == 0)
        for ind in zero_ind:
            try:
                if day_list[ind-1] and day_list[ind+1] == 1:
                    day_list[ind] = 1
            except IndexError:
                if ind == 0:
                    if day_list[-1] and day_list[1] == 1:
                        day_list[0] = 1
                elif ind == len(day_list)-1:
                    if day_list[len(day_list)-2] and day_list[0] == 1:
                        day_list[len(day_list)-1] = 1
    return num_emperor

def my_mean(values):
    n = 0
    summ = 0.0
    for value in values:
        summ += value
        n += 1
    return summ/n

def monte_carlo(iters, D):
    iter = 0
    n_emperor = 0
    while iter < iters:
        n_emperor = Exp(D)
        yield n_emperor
        iter += 1

avg_n_emperor = my_mean(monte_carlo(iters,D))
print(avg_n_emperor)

And my logic is as follow:

For the day_list inside the Exp(D) function, where D is the number of days in a year, zeros mean no holiday, and ones mean holiday. Initially the day_list is all zeros since there is no holiday to begin with.

The rules of defining a random day (d) as a holiday is as follow:

1. At the beginning of the reign of the current Emperor, his birthday is declared a holiday from that year onwards.

2. If both the day before and after a day d are holidays, then d also becomes a holiday.

I then subsequently implement the rules stated for the question, to gradually add holidays (ones) into the day_list. After num_emperor number of emperors, all the days (d) in day_list will become 1, i.e. all days will become holiday. This is the point to quit the while_loop in Exp(D) function and count the number of emperors required. To get the average number of emperors required for all the days to become holidays (avg_n_emperor), I then apply the monte-carlo method.

For my current code, the time takes is as follow:

avg_n_emperor = my_mean(monte_carlo(iters=100000,D=5)) #6-7 seconds

avg_n_emperor = my_mean(monte_carlo(iters=1000000,D=5)) #about 62 seconds

in which the time takes increase approx. linearly with the iters.

However,

avg_n_emperor = my_mean(monte_carlo(iters=1000,D=365)) #about 68 seconds

already takes about 68 seconds, and the question is asking for D=10000. Not to mention that the iters required for the answer to be accurate within 4 digits after the decimal points (as required by the question) would be much larger than 1000000 too...

Any suggestion(s) to speed up my code would be appreciated! :)

Answer 1

Welcome to Code Review. Nice implementation, easy to read and understand.

Optimization

There are some "expensive" operations that can be simplified. Below I commented the relevant parts:

def Exp(D):
    # the method "all" takes O(D)
    while all((d == 1 for d in day_list)) == False:
        # O(D)
        zero_ind = (i for i,v in enumerate(day_list) if v == 0) 
        # O(D)
        for ind in zero_ind:
            # Here there are only O(1) operations
    return num_emperor

By \$O(D)\$ I mean that in the worst case such operation will iterate D times, where D is the number of days.

The condition in the while loop can be simplified by checking if the number of holidays is < days:

def Exp(D):
    holidays = 0
    while holidays < D:
        # increment holidays 
    return num_emperor

The second optimization is to avoid the inner loops. Once the new birthday has been calculated, it's enough to "look around" that specific day:

def Exp(D):
    # ..
    while holidays < D:
        bday = np.random.randint(0,D)
        # Increment holidays only if birthday is not in a holiday
        if day_list[bday] == 0:
            holidays += 1
        day_list[bday] = 1
        num_emperor+=1

        yesterday = (bday - 1) % D
        day_before_yesterday = (bday - 2) % D
        if day_list[day_before_yesterday] == 1 and day_list[yesterday] == 0:
            day_list[yesterday] = 1
            holidays += 1

        tomorrow = (bday + 1) % D
        day_after_tomorrow = (bday + 2) % D
        if day_list[day_after_tomorrow] == 1 and day_list[tomorrow] == 0:
            day_list[tomorrow] = 1
            holidays += 1
    return num_emperor

The % operator prevents to overflow the array, so you don't need to catch exceptions.

Running the average:

avg_n_emperor = my_mean(monte_carlo(iters=1000,D=365))
# Output: 1173.786
# Running time: around 2 seconds

Regarding the style, @Peilonrayz already provided an excellent review.

Answer 2

Firstly lets get your code to be a little cleaner:

• You can use statistics.mean rather than make my_mean.

• You should use a for loop rather than a while loop in monte_carlo.

• You don't need to do assign n_emperer at all in the function.

• Exp and D should be lower_snake_case. This is as they are functions and variables.

• You should put spaces around all operators.

• There should be a space after commas.

• You should have some better names day_list could just be days, D could also be something like days, summ can be total, iters could be amounts.

• You can just use all(day_list) rather than all((d == 1 for d in day_list)).

• Do not use == to compare to singletons like False. It would be better if you instead use not.

• This doesn't check if both of the values are 1 it checks if the first is truthy and the second is one. This means if you set day_list[index - 1] to two it'd still be true.

  day_list[ind - 1] and day_list[ind + 1] == 1
  

  To check they are both equal to one you shoud use:

  day_list[ind - 1] == 1 and day_list[ind + 1] == 1
  

  Here I would instead just check if they are truthy.

• You don't need if ind == 0: as if ind is 0 then ind - 1 will be -1.

• You can just use (ind + 1) % len(days) to remove the need for elif index == len(days)-1:.

import random
import statistics


def simulate(days_in_year):
    days = [0] * days_in_year
    emperors = 0
    while not all(days):
        days[random.randrange(len(days))] = 1
        emperors += 1
        for index, value in enumerate(days):
            if value:
                continue
            if days[index - 1] and days[(index + 1) % len(days)]:
                days[index] = 1
    return emperors


def monte_carlo(amount, days):
    for _ in range(amount):
        yield simulate(days)


print(statistics.mean(monte_carlo(amount, days)))

Now that the code is nice and small we can focus on what is causing performance issues.

1. The following any runs in \$O(n)\$ time, where \$n\$ is the length of days. This means worst case it will run however long days is each time you call it.

   not all(days)
   

   We can do better than that by adding a variable in that increments each time we change a 0 to a 1. We can then compare that to days_in_year to see if the list is full. This will run in \$O(1)\$ time causing a significant saving.

2. If a new emperor is born on an already existing holiday then no extra holidays will be made.

3. When a new emperor is born you don't need to check whether each zero can be changed, you instead only need to check two. This will cut another \$O(n)\$ operation to \$O(1)\$.
   Say we have the following as days:

   0123456
   1000010
   

   If the new birthday is:

   • 6 - Because both 5 and 0 are already 1s no additional holidays can be made.

   • 3 - Because 4 is a 0 and 5 is a 1, 4 can become a 1. Because 2 is a 0 but 1 is a 0 then 3 cannot become a 1.

     This cannot propagate outwards.

Answer 3

Realistically, my review would have to be "That's not gonna work, you won't get the required accuracy with such an experiment. You need a different approach".

But here's an O(D) time simulation. Instead of potentially generating already occurred birthdays over and over again, I focus just on new birthdays. That is, I shuffle all possible birthdays at the start, and then I just go through them. Of course that means I can't just do emperors += 1. Instead, I add the expected number of new emperors needed to come across a new birthday.

With 1000 simulations, it takes my laptop about 0.6 seconds for D=365, 1.8 seconds for D=1000, or 19 seconds for D=10000.

from random import sample
from statistics import mean

def Exp(D):
    emperors = 0
    holidays = set()
    for i, day in enumerate(sample(range(D), D)):
        emperors += D / (D - i)
        holidays.add(day)
        if (day + 2) % D in holidays:
            holidays.add((day + 1) % D)
        if (day - 2) % D in holidays:
            holidays.add((day - 1) % D)
        if len(holidays) == D:
            return emperors

print(mean(Exp(365) for _ in range(1000)))

Meh. Just tried it the emperor += 1 way as well, that took about 1.35, 4.1 and 62 seconds instead:

from random import randrange
from statistics import mean

def Exp(D):
    emperors = 0
    holidays = set()
    while len(holidays) < D:
        emperors += 1
        day = randrange(D)
        if day not in holidays:
            holidays.add(day)
            if (day + 2) % D in holidays:
                holidays.add((day + 1) % D)
            if (day - 2) % D in holidays:
                holidays.add((day - 1) % D)
    return emperors

print(mean(Exp(365) for _ in range(1000)))