Projeto Euler # 645 - acelere a simulação de Monte-Carlo em Python

Question

Estou tentando resolver o Q645. Embora a lógica usada para meu código pareça apropriada, o código em si é muito lento para o grande número exigido nesta questão. Posso pedir sugestões para melhorar o desempenho do meu código?

A pergunta é como no link: https://projecteuler.net/problem=645

Meu código Python é o seguinte:

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)

E minha lógica é a seguinte:

Para a day_list dentro da função Exp (D) , onde D é o número de dias em um ano, zeros significam nenhum feriado e uns significam feriado. Inicialmente, a day_list é composta apenas por zeros, já que não há feriado para começar .

As regras para definir um dia aleatório ( d ) como feriado são as seguintes:

1. No início do reinado do atual imperador, seu aniversário é declarado feriado daquele ano em diante.

2. Se o dia anterior e o posterior ao dia d forem feriados, d também se tornará feriado.

Em seguida, implemento posteriormente as regras estabelecidas para a pergunta, para adicionar gradualmente feriados (uns) à day_list . Após num_emperor número de imperadores, todos os dias ( d ) em day_list se tornarão 1, ou seja , todos os dias se tornarão feriados . Este é o ponto para sair do while_loop na função Exp (D) e contar o número de imperadores necessários. Para obter o número médio de imperadores necessários para que todos os dias se tornem feriados ( avg_n_emperor ), aplico o método do monte-carlo.

Para meu código atual, o tempo gasto é o seguinte:

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

em que o tempo leva aumentar aprox. linearmente com os iters .

Contudo,

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

já leva cerca de 68 segundos, e a pergunta está pedindo D = 10000. Sem mencionar que os iters necessários para que a resposta seja precisa dentro de 4 dígitos após as casas decimais (conforme exigido pela pergunta) também seriam muito maiores do que 1000000 ...

Qualquer sugestão (ões) para acelerar meu código seria apreciada! :)

Answer 1

Bem-vindo à revisão do código. Boa implementação, fácil de ler e entender.

Otimização

Existem algumas operações "caras" que podem ser simplificadas. Abaixo comentei as partes relevantes:

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

Por \$O(D)\$Quero dizer que, no pior caso, essa operação irá iterar Dvezes, onde Dé o número de dias.

A condição no loop while pode ser simplificada verificando se o número de feriados é <dias:

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

A segunda otimização é evitar os loops internos. Uma vez que o novo aniversário foi calculado, basta "dar uma olhada" naquele dia específico:

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

O %operador evita o estouro do array, então você não precisa capturar exceções.

Executando a média:

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

Em relação ao estilo, @Peilonrayz já fez uma excelente crítica.

Answer 2

Em primeiro lugar, vamos deixar seu código um pouco mais limpo:

• Você pode usar em statistics.meanvez de fazer my_mean.

• Você deve usar um forloop em vez de um loop while monte_carlo.

• Você não precisa atribuir n_emperernada na função.

• Expe Ddeveria ser lower_snake_case. É assim que são funções e variáveis.

• Você deve colocar espaços ao redor de todos os operadores.

• Deve haver um espaço após as vírgulas.

• Você deve ter alguns nomes melhores, day_listpoderia ser days, Dtambém poderia ser algo como days, summpode ser total, iterspoderia ser amounts.

• Você pode apenas usar em all(day_list)vez de all((d == 1 for d in day_list)).

• Não use ==para comparar com singletons como False. Seria melhor se você usar not.

• Isso não verifica se ambos os valores são 1, ele verifica se o primeiro é verdadeiro e o segundo é um. Isso significa que se você definir day_list[index - 1]como dois, ainda será verdade.

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

  Para verificar se ambos são iguais a um que você deve usar:

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

  Em vez disso, gostaria de apenas verificar se eles são verdadeiros.

• Você não precisa if ind == 0:como se indfosse 0, então ind - 1será -1.

• Você pode apenas usar (ind + 1) % len(days)para remover a necessidade de 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)))

Agora que o código é bom e pequeno, podemos nos concentrar no que está causando problemas de desempenho.

1. O seguinte anyé executado em \$O(n)\$tempo, onde \$n\$é o comprimento de days. Isso significa que, na pior das hipóteses, ele será executado, independentemente dos longos dias de cada vez que você ligar.

   not all(days)
   

   Podemos fazer melhor do que isso adicionando uma variável a esses incrementos cada vez que mudamos de 0 para 1. Podemos então comparar days_in_yearpara ver se a lista está cheia. Isso será executado em \$O(1)\$ tempo causando uma economia significativa.

2. Se um novo imperador nascer em um feriado já existente, nenhum feriado extra será concedido.

3. Quando um novo imperador nasce, você não precisa verificar se cada zero pode ser alterado, ao invés disso, você só precisa verificar dois. Isso vai cortar outro \$O(n)\$operação para \$O(1)\$.
   Digamos que temos o seguinte como days:

   0123456
   1000010
   

   Se o novo aniversário for:

   • 6 - Como 5 e 0 já são 1s, nenhum feriado adicional pode ser concedido.

   • 3 - Como 4 é 0 e 5 é 1, 4 pode se tornar 1. Como 2 é 0, mas 1 é 0, então 3 não pode se tornar 1.

     Isso não pode se propagar para fora.

Answer 3

Realisticamente, minha revisão teria que ser "Isso não vai funcionar, você não obterá a precisão necessária com tal experimento. Você precisa de uma abordagem diferente" .

Mas aqui está uma simulação de tempo O (D). Em vez de potencialmente gerar aniversários já ocorridos repetidamente, concentro-me apenas em novos aniversários. Ou seja, embaralhei todos os aniversários possíveis no início e depois simplesmente os examino. Claro que isso significa que não posso simplesmente fazer emperors += 1. Em vez disso, adiciono o número esperado de novos imperadores necessários para fazer um novo aniversário.

Com 1000 simulações, meu laptop leva cerca de 0,6 segundos para D = 365, 1,8 segundos para D = 1000 ou 19 segundos para 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. Também tentei da mesma emperor += 1forma, que demorou cerca de 1,35, 4,1 e 62 segundos:

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