#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  9 14:07:57 2023

@author: matthieu

"""



#justifier la fonctionG de transfert







# idées pour un TP sur simulation de vad
import numpy.random  as rd
# nouvelle version numpy.random.Generator
import numpy as np
import matplotlib.pyplot as plt



def Graphe(L):
    try:
        unique, nombre = np.unique(L, return_counts=True)
        plt.stem(unique,nombre/len(L),linefmt='purple',basefmt='purple')
        plt.show()
    except:
        raise ValueError

def Repetition(N,f,*args):
    ''' renvoie un tableau de N valeurs en utilisant le générateur f'''
    L=np.zeros(N)
    for i in range(N):
        L[i]=f(*args)
    return L

    
def Moyenne(S):
    '''renvoie la moyenne de la série 
    cf np.mean et np.sum
    '''
    somme=0
    for x in S:
        somme+=x
    return somme/len(S)


def EcartType(S):
    '''
    cf np.stdv
    en utilisant les operations "elementwise"
    '''     
    return np.sqrt(Moyenne(S**2)-Moyenne(S)**2)

#%%
def Ber(p):
    ''' cf rd.binomial'''
    return rd.random()<p


def Bin(n,p):
    ''' cf rd.binomial'''

    return np.sum(Repetition(n, Ber,p))


def Bin(n,p):
    ''' cf rd.binomial'''

    S=0
    for i in range(n):
        S+=Ber(p)
    return S

def Geometrique(p):
    '''cf rd.geometric'''
    n=0
    S=0 
    while S==0:
        S+=Ber(p)
        n+=1
    return n 




#loi de Poisson a voir lors du cours VA à densité







#%% essai

R=Repetition(100, Bin,10,0.75)
Graphe(R)




def EstValide(P):
    S=0
    Test=True
    i=0
    n=len(P)
    while Test and i<n:
        S+=P[i]
        Test =(0<P[i]<=1)
        i+=1
    if Test==False or abs(S-1)>10**-5:
        return False
    else:
        return True


#%% 

def Repartition(P):
    assert EstValide(P), "tableau non valide"
    R=np.zeros(len(P),float)
    R[0]=P[0]
    
    for i in range(1,len(P)):
        R[i]=R[i-1]+P[i]
    
    return R
    

def Simulation(R):
    i=0
    U=rd.random()
    while i<len(R) and (U >R[i]) :
        i+=1
    return i








def Choisir(S,P):
    assert len(S)==len(P), "les tableaux ne sont pas de même longueur"
    return S[Simulation(Repartition(P))]




def Riemann():
    C=6/(np.pi)**2
    S=0
    i=0
    u=rd.random()
    while u>S:
        i+=1
        S+=C/i**2
       
    return i-1


#%%
N=10**6
R=Repetition(N,Riemann)
# on limite les valeurs extrèmes pour avoir un graphe lisible
Rp=np.zeros(len(R))
for i in range(len(R)):
    Rp[i]=min(R[i],10)

Graphe(Rp)
#calcul de la moyenne sur les données initiales
#si on augmente la taille de l'échantillon,  la somme  semble tendre vers + infty
# les résultats vairient beaucoup
Moyenne(R)







#%%
P=[0.1,0.3,0.5,0.1]
R=Repartition(P)
Graphe(Repetition(10000, Simulation,R))


#%%
N=10
P=np.ones(N)/N
R=Repartition(P)
Graphe(Repetition(10000, Simulation,R))







#%% Autres lois classiques hors programme
#Loi de Pascal

def Pascal(n,p):
    rang=0
    succes=0
    while succes<n:
        succes+=Ber(p)
        rang+=1
    return rang
    


# binome négatif 

def BinNegative(n,p):
    succes=0
    echecs=0
    while succes<n:
        if Ber(p)==1:
            succes+=1
        else:
            echecs+=1
    return echecs

#Loi hypergéométrique

def Hypergeometrique(n,a,b):
    assert n<=a+b, "Problème pas assez de boules"
    NbNoires=a
    NbBlanches=b
    Succes=0
    for i in range(n):
        if Ber(NbNoires/(NbNoires+NbBlanches))==1:
            NbNoires-=1
        else:
            NbBlanches-=1
            Succes+=1
    return Succes
        
        
        
    
    
    
#Urne de Polya
def Polya(n,a,b,h):
    NbNoires=a
    NbBlanches=b
    Succes=0
    for i in range(n):
        if Ber(NbNoires/(NbNoires+NbBlanches))==1:
            NbNoires+=h
        else:
            NbBlanches+=h
            Succes+=1
    return Succes

#%%
Mystere1=Repetition(1000,Bin,20,0.75)
Mystere2=Repetition(1000,Geometrique,0.1)
Mystere3=Repetition(1000,rd.poisson,4)
Graphe(Mystere3)