"""
author: PackSciences
"""
from scipy.special import comb as C
from math import floor
from collections import OrderedDict
import matplotlib.pyplot as plt


def sommande(j, p, n, m):
    """
    Calcule le sommande dans l'expression du modèle
    en fonction des différents paramètres
    Paramètres:
    j : indice muet de sommation
    p : probabilité de l'expérience de Bernoulli
    n : nombre total d'expériences de Bernoulli réalisé
    m : longueur de la chaine maximale
    """
    if j > floor(n/m):
        raise ValueError("Erreur sur la valeur de j: {}".format(j))
    elif j < 1:
        raise ValueError("Erreur sur la valeur de j: {}".format(j))
    if p > 1 or p < 0:
        raise ValueError("Erreur sur la valeur de p: {}".format(p))
    res = (-1)**(j+1)
    res *= (p + ((n-(j*m)+1)/j)*(1-p))
    res *= C(n-j*m, j-1)
    res *= p**(j*m)
    res *= (1-p)**(j-1)
    return res


def de_moivre(p, n, m):
    """
    Réalise la somme du modèle de de Moivre pour la longueur de la plus
    longue chaîne pour n expériences de Bernoulli.
    Correspond à P(l_n >= m) (probabilité cumulative)
    Paramètres:
    p : probabilité de l'expérience de Bernoulli
    n : nombre total d'expériences de Bernoulli réalisé
    m : longueur de la chaine maximale
    """
    if n > 1000:
        raise ValueError("L'algorithme diverge numériquement")
    somme = 0
    for j in range(1, floor(n/m) + 1):
        somme += sommande(j, p, n, m)
    if somme > 1:
        if n > 100:
            # la remise à 1 est nécessaire pour des valeurs de n > 100
            somme = 1
        else:
            raise ValueError("La somme trouvée est supérieure à 1 "
                             "pour n < 100")
    return somme


def graphique_cumulative(p, n_min, n_max, n_pas, m):
    """
    Renvoie la probabilité cumulative en fonction du nombre d'expériences n
    de n_min à n_max avec un pas n_pas, d'avoir une longueur >= m, à p fixé
    """
    plt.close()
    dic_resu = OrderedDict()
    for n in range(n_min, n_max, n_pas):
        dic_resu[n] = de_moivre(p, n, m)
    plt.scatter(list(dic_resu.keys()),
                [dic_resu[cle] for cle in dic_resu.keys()],
                label="$P(X_{n} ≥ m)$ as a function of n")
    plt.title("Probability to obtain a chain of length ≥ {0:}"
              "\n as a function of the number of Bernoulli trials realized"
              " of probability {1:.2%}".format(m, p))
    plt.ylabel("Cumulative probability")
    plt.xlabel("Number of trials")
    plt.legend()
    plt.show()
    return dic_resu

p = 0.93
n_min = 2
n_max = 1000
pas = 1
m = 100

graphique_cumulative(p, n_min, n_max, pas, m)