Project Euler # 645 - ускорить моделирование методом Монте-Карло на Python

Question

Я пытаюсь решить Q645. Хотя логика, использованная для моего кода, кажется подходящей, сам код слишком медленный для большого числа, требуемого в этом вопросе. Могу я запросить предложения по улучшению производительности моего кода?

Вопрос как по ссылке: https://projecteuler.net/problem=645

Мой код Python выглядит следующим образом:

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)

И моя логика следующая:

Для day_list внутри функции Exp (D) , где D - количество дней в году, нули означают отсутствие праздников, а единицы - праздники. Первоначально day_list состоит из нулей, так как для начала нет выходных .

Правила определения случайного дня ( d ) как праздника следующие:

1. В начале правления нынешнего императора его день рождения объявлен праздником с этого года.

2. Если и день до и после дня d - праздничные дни, то d также становится праздником.

Затем я применяю правила, указанные для вопроса, чтобы постепенно добавлять праздники (праздники) в day_list . После num_emperor числа императоров все дни ( d ) в day_list станут 1, т.е. все дни станут праздничными . Здесь нужно выйти из цикла while_loop в функции Exp (D) и подсчитать необходимое количество императоров. Чтобы получить среднее количество императоров, необходимое для того, чтобы все дни стали праздниками ( avg_n_emperor ), я затем применяю метод Монте-Карло.

Для моего текущего кода время занимает следующее:

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

при котором время увеличивается прибл. линейно с итерами .

Тем не мение,

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

уже занимает около 68 секунд, а вопрос задает D = 10000. Не говоря уже о том, что итоги, необходимые для того, чтобы ответ был точным в пределах 4 цифр после десятичных знаков (как того требует вопрос), тоже будут намного больше, чем 1000000 ...

Мы будем благодарны за любые предложения по ускорению моего кода! :)

Answer 1

Добро пожаловать в Code Review. Хорошая реализация, легко читается и понимается.

Оптимизация

Есть некоторые «дорогие» операции, которые можно упростить. Ниже я прокомментировал соответствующие части:

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

Автор \$O(D)\$Я имею в виду, что в худшем случае такая операция будет повторяться Dраз, где D- количество дней.

Условие в цикле while можно упростить, проверив, составляет ли количество праздников <дней:

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

Вторая оптимизация - избежать внутренних циклов. После того, как новый день рождения рассчитан, достаточно «осмотреться» в этот конкретный день:

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

В %предотвращаете оператора переполнение массива, так что вам не нужны исключения улова.

Среднее значение:

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

Что касается стиля, @Peilonrayz уже предоставил отличный обзор.

Answer 2

Во-первых, давайте сделаем ваш код немного чище:

• Вы можете использовать, statistics.meanа не делать my_mean.

• Вам следует использовать forцикл, а не цикл while monte_carlo.

• Вам вообще не нужно делать assign n_empererв функции.

• Expтак и Dдолжно быть lower_snake_case. Это так как они являются функциями и переменными.

• Вы должны поставить пробелы вокруг всех операторов.

• После запятых должен быть пробел.

• Вы должны иметь несколько лучших имен, day_listможет быть days, Dможет быть что-то вроде days, summможет быть total, itersможет быть amounts.

• Вы можете просто использовать all(day_list)вместо all((d == 1 for d in day_list)).

• Не используйте ==для сравнения с одиночками вроде False. Было бы лучше, если бы вы вместо этого использовали not.

• Это не проверяет, равны ли оба значения 1, он проверяет, является ли первое истинным, а второе - единым. Это означает, что если вы установите day_list[index - 1]два, это все равно будет правдой.

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

  Чтобы проверить, что они оба равны одному, вы должны использовать:

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

  Вместо этого я бы просто проверил, правдивы ли они.

• Вам не нужно, if ind == 0:как если бы было ind0, тогда ind - 1будет -1.

• Вы можете просто использовать, (ind + 1) % len(days)чтобы избавиться от необходимости 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)))

Теперь, когда код красивый и небольшой, мы можем сосредоточиться на том, что вызывает проблемы с производительностью.

1. Следующее anyработает в \$O(n)\$время, где \$n\$это длина days. Это означает, что в худшем случае он будет работать, сколько бы дней ни было каждый раз, когда вы его вызываете.

   not all(days)
   

   Мы можем добиться большего, добавляя переменную с приращением каждый раз, когда мы меняем 0 на 1. Затем мы можем сравнить это с, days_in_yearчтобы увидеть, заполнен ли список. Это будет работать в \$O(1)\$ время, что дает значительную экономию.

2. Если в уже существующий праздник рождается новый император, то дополнительных выходных не будет.

3. Когда рождается новый император, вам не нужно проверять, можно ли изменить каждый ноль, вам нужно только проверить два. Это сократит другой \$O(n)\$операция к \$O(1)\$.
   Скажем, у нас есть следующее days:

   0123456
   1000010
   

   Если новый день рождения:

   • 6 - Поскольку и 5, и 0 уже равны 1, дополнительные выходные не могут быть выполнены.

   • 3 - Поскольку 4 - это 0, а 5 - это 1, 4 может стать 1. Поскольку 2 - это 0, а 1 - это 0, то 3 не может стать 1.

     Это не может распространяться наружу.

Answer 3

На самом деле, мой отзыв должен был бы звучать так: «Это не сработает, вы не получите требуемой точности с таким экспериментом. Вам нужен другой подход» .

Но вот симуляция времени O (D). Вместо того, чтобы потенциально генерировать уже произошедшие дни рождения снова и снова, я сосредотачиваюсь только на новых днях рождения. То есть я вначале перемешиваю все возможные дни рождения, а потом просто их перебираю. Конечно, это означает, что я не могу просто так делать emperors += 1. Вместо этого я добавляю ожидаемое количество новых императоров, необходимое для празднования нового дня рождения.

При 1000 симуляциях моему ноутбуку требуется около 0,6 секунды для D = 365, 1,8 секунды для D = 1000 или 19 секунд для 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)))

Мех. Просто попробовал emperor += 1, это заняло около 1,35, 4,1 и 62 секунды:

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)))