Proyecto Euler # 645: acelere la simulación Monte-Carlo en Python

Question

Estoy intentando resolver Q645. Si bien la lógica utilizada para mi código parece ser apropiada, el código en sí es demasiado lento para el gran número requerido en esta pregunta. ¿Puedo pedir sugerencias para mejorar el rendimiento de mi código?

La pregunta es como en el enlace: https://projecteuler.net/problem=645

Mi código de Python es el siguiente:

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)

Y mi lógica es la siguiente:

Para day_list dentro de la función Exp (D) , donde D es el número de días en un año, los ceros significan que no hay vacaciones y los unos significan vacaciones. Inicialmente, day_list es todo ceros ya que no hay vacaciones para empezar .

Las reglas para definir un día aleatorio ( d ) como feriado son las siguientes:

1. Al comienzo del reinado del actual Emperador, su cumpleaños se declara festivo a partir de ese año.

2. Si tanto el día anterior como el posterior al día d son festivos, entonces d también se convierte en festivo.

Posteriormente, implemento las reglas establecidas para la pregunta, para agregar gradualmente días festivos (unos) en la lista de días . Después de num_emperor número de emperadores, todos los días ( d ) en day_list se convertirán en 1, es decir , todos los días se convertirán en festivos . Este es el punto para salir de while_loop en la función Exp (D) y contar el número de emperadores necesarios. Para obtener el número medio de emperadores necesarios para que todos los días se conviertan en festivos ( avg_n_emperor ), aplico el método montecarlo.

Para mi código actual, el tiempo es el siguiente:

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

en el que el tiempo tarda en aumentar aprox. linealmente con los iters .

Sin embargo,

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

ya lleva unos 68 segundos, y la pregunta es D = 10000. Sin mencionar que los iters requeridos para que la respuesta sea precisa dentro de los 4 dígitos después de los puntos decimales (como lo requiere la pregunta) también serían mucho mayores que 1000000 ...

¡Cualquier sugerencia para acelerar mi código será apreciada! :)

Answer 1

Bienvenido a Revisión de código. Buena implementación, fácil de leer y comprender.

Mejoramiento

Hay algunas operaciones "caras" que se pueden simplificar. A continuación les comento las 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)\$Quiero decir que, en el peor de los casos, dicha operación se repetirá Dveces, donde Destá el número de días.

La condición en el ciclo while se puede simplificar comprobando si el número de días festivos es <días:

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

La segunda optimización es evitar los bucles internos. Una vez que se ha calculado el nuevo cumpleaños, es suficiente "mirar alrededor" ese día 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

El %operador evita que se desborde la matriz, por lo que no es necesario detectar excepciones.

Ejecutando el promedio:

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

En cuanto al estilo, @Peilonrayz ya brindó una excelente reseña.

Answer 2

En primer lugar, hagamos que su código sea un poco más limpio:

• Puede usar en statistics.meanlugar de hacer my_mean.

• Debería utilizar un forbucle en lugar de un bucle while monte_carlo.

• No es necesario asignar n_emperernada en la función.

• Expy Ddebería ser lower_snake_case. Esto es porque son funciones y variables.

• Debe poner espacios alrededor de todos los operadores.

• Debe haber un espacio después de las comas.

• Debería tener algunos nombres mejores. day_listPodría ser days, Dtambién podría ser algo como days, summpuede ser total, iterspodría ser amounts.

• Puede usar en all(day_list)lugar de all((d == 1 for d in day_list)).

• No lo utilice ==para comparar con singletons como False. Sería mejor si usaras not.

• Esto no verifica si ambos valores son 1, verifica si el primero es verdadero y el segundo es uno. Esto significa que si se establece day_list[index - 1]en dos, seguirá siendo cierto.

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

  Para comprobar que ambos son iguales al que debe utilizar:

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

  Aquí, en cambio, solo comprobaría si son veraces.

• No necesita if ind == 0:como si fuera ind0, entonces ind - 1será -1.

• Puede usar (ind + 1) % len(days)para eliminar la necesidad 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)))

Ahora que el código es agradable y pequeño, podemos centrarnos en las causas de los problemas de rendimiento.

1. Lo siguiente se anyejecuta en \$O(n)\$tiempo, donde \$n\$es la longitud de days. Esto significa que, en el peor de los casos, se ejecutará sin importar los días que sean cada vez que lo llame.

   not all(days)
   

   Podemos hacerlo mejor agregando una variable que se incremente cada vez que cambiamos un 0 a un 1. Luego podemos comparar eso days_in_yearpara ver si la lista está llena. Esto se ejecutará en \$O(1)\$ tiempo provocando un ahorro significativo.

2. Si un nuevo emperador nace en un feriado ya existente, no se realizarán feriados adicionales.

3. Cuando nace un nuevo emperador, no es necesario verificar si se puede cambiar cada cero, en su lugar, solo necesita verificar dos. Esto cortará otro \$O(n)\$operación a \$O(1)\$.
   Digamos que tenemos lo siguiente como days:

   0123456
   1000010
   

   Si el nuevo cumpleaños es:

   • 6 - Debido a que tanto el 5 como el 0 ya son 1, no se pueden hacer vacaciones adicionales.

   • 3 - Dado que 4 es un 0 y 5 es un 1, 4 puede convertirse en un 1. Como 2 es un 0 pero 1 es un 0, entonces 3 no puede convertirse en un 1.

     Esto no se puede propagar hacia afuera.

Answer 3

Siendo realistas, mi revisión tendría que ser "Eso no va a funcionar, no obtendrás la precisión requerida con un experimento de este tipo. Necesitas un enfoque diferente" .

Pero aquí hay una simulación de tiempo O (D). En lugar de generar potencialmente cumpleaños ya ocurridos una y otra vez, me concentro solo en nuevos cumpleaños. Es decir, mezclo todos los cumpleaños posibles al principio y luego los reviso. Por supuesto que eso significa que no puedo simplemente hacer emperors += 1. En cambio, agrego la cantidad esperada de nuevos emperadores necesarios para encontrarse con un nuevo cumpleaños.

Con 1000 simulaciones, mi computadora portátil tarda aproximadamente 0,6 segundos para D = 365, 1,8 segundos para D = 1000 o 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. También lo probé de la misma emperor += 1manera, eso tomó alrededor de 1.35, 4.1 y 62 segundos en su lugar:

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