#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 30 16:49:56 2025

@author: matthieu
"""
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rd
plt.rcParams['figure.dpi'] = 500

S=30
v0=5
L=[v0]
v=v0
for k in range(20):
    v=v+v*(S-v)/(2*S)
    L.append(v)
plt.plot(L, marker="*",linestyle='none')
plt.xlabel("n")
plt.ylabel("v_n")
plt.show()

## plusieurs graphiques
#%%
S=30
for v0 in range(0,40,5):
    L=[v0]
    v=v0
    for k in range(20):
        v=v+v*(S-v)/(2*S)
        L.append(v)
    plt.plot(L, marker="*",linestyle='none')
    plt.xlabel("n")
    plt.ylabel("v_n")
plt.show()

## méthode euler
#%%
A=5
S=30

def F(y):
    return y*(S-y)*(y-A)/(2*S*S)


def euler(f,y0,t_final,n):
    """méthode d'Euler
        f fonction définissant l' équation diff
         y0 valeur initiale en t= 0
         t_final borne de l'intervalle de calcul 
         n nombre de pas
         renvoie la liste des temps et la liste des valeurs de y
    """
    pas= t_final/n
    resultat=[y0]
    y=y0
    for i in range(1,n):
        y=f(y)*pas +y
        resultat.append(y)
    temps=np.linspace(0,t_final,n)
    return temps, resultat      
##
#%%
for y0 in range(0,35):
    t,r=euler(F, y0, 20, 10**5) 
    plt.plot(t,r)   
plt.show()

##
#%%
def galton_watson(lambda_,n):
    population = np.zeros(n+1)
    population[0] = 1
    Z = 1
    for i in range(1,n+1):
        descendants = 0
        for j in range(Z):
            descendants += rd.poisson(lambda_)
        population[i] = descendants
        Z = descendants
        if descendants == 0:
            return population
    return population

##
#%%
lambda_=1.1
for i in range(10):
    plt.plot(galton_watson(lambda_, 20),marker="+")
plt.show()


##
#%%
def galton_watson2(lambda_,n):
    population = np.zeros(n+1)
    population[0] = 1
    Z = 1
    for i in range(1,n+1):
        descendants = 0
        for j in range(Z):
            descendants += rd.poisson(lambda_)
        population[i] = descendants
        Z = descendants
        if descendants == 0:
            return 1
    return 0


def extinction(lambda_):
    N=5000
    S=0
    n=60
    for i in range(N):
        S+=galton_watson2(lambda_,n)
    return S/N

Llambda_=np.arange(0.8,1.2,0.01)
Lplambda_=[extinction(lambda_) for lambda_ in Llambda_ ]
plt.plot(Llambda_,Lplambda_,color='goldenrod')
plt.show()