Code projet numérique M2
Article
Séance 09/11/2021
Première ébauche bactérie libre, sans interaction. C’est un mouvement Brownien avec une longueur moyenne de persistence.
#%% Modalités
'''
Oral : 15 minutes (10+5 qustions), faire intro puis montrer résultats en disant méthodes utilisées, problèmes rencontrés... Pas grave si pas abouti, ce qui importe est le raisonnement. Il faudra aussi donner le code, donc le commenter et le fournir avec un fichier mode d'emploi, citer les sources si on prend un code déjà fait.
'''
#%% Modalités
'''
Oral : 15 minutes (10+5 qustions), faire intro puis montrer résultats en disant méthodes utilisées, problèmes rencontrés... Pas grave si pas abouti, ce qui importe est le raisonnement. Il faudra aussi donner le code, donc le commenter et le fournir avec un fichier mode d'emploi, citer les sources si on prend un code déjà fait.
'''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
import scipy.optimize as spo
plt.rcParams['axes.grid'] = True #Paramètre esthétique, peut être retiré
plt.rcParams['savefig.format'] = 'png' #On peut sélectionner pdf si on aime les images vectorielles
#%% Définition des paramètres
p=0.5 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
N=100000# Nombre des pas de la simulation
x0,y0=0.5,0.5# Conditions initiales
l0=0.001#Pas de longueurs l0 (pourra être aléatoire par la suite)
#%% Calcul trajectoires sans gradient, une bactérie
def Trajectoire(N,p):
angles=np.random.uniform(0,2*np.pi,N) #Angles à chaque pas de temps
changes=np.random.binomial(1,p,N) # Définit si on chnge de direction pour ce pas de temps, p est la proba de 1, donc de tourner.
changes[0]=1 # Artificiel mais permet de faire marcher le code, ça revient à définir deux conditions initiales
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
x,y=np.zeros(N),np.zeros(N)
x[0]=x0
y[0]=y0
for i in range(0,N-1): # Calcul des trajectoires, avec pour l'instant, pas de CAL
if changes[i]==1 : #On change de direction
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changes[i]==0 : # On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
return x,y
#%%Tracey
x,y=Trajectoire(N,p)[0],Trajectoire(N,p)[1]
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot([x0],[y0],color='r',label='Position initiale',marker='o')
plt.plot([x[N-1]],[y[N-1]],color='g',label='Position finale',marker='o')
plt.xlim([0.2,0.8])
plt.ylim([0.2,0.8])
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une probabilité de changer de direction p="+str(p))
plt.show()
#%% MSD : vérification du mouvement diffusif, construction des MSD
n=200# Nombre de mesures pour moyenner
MSD=np.zeros((N,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
for j in range(n):
x,y=Trajectoire(N,p)[0]-x0,Trajectoire(N,p)[1]-y0
MSD[:,j]=x[:]**2+y[:]**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
#%% Traitement statistique pour tracer MSD
moy=np.zeros(N)
erreurs=np.zeros(N)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(N):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot
temps=np.array([i for i in range(N)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences, p="+str(p))
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
#%% Ajustation du MSD à grand t : on part de N//10 on en aura beaucoup, l'erreur est sur le log, on naura pas log(0) en toute logique.
absc=temps[N//100:N]
ordo=moy[N//100:N]
yerreur=erreurs[N//100:N]/ordo# Erreur du log, sera utile.
xerreur=np.zeros(N-N//100)
'''Copie d'un code déjà utilisé l'année dernière. Ici on trace logMSD en fonction de log(t) pour avoir les pentes et ordonnées origine '''
x = np.log(absc)
y = np.log(ordo)
ux= xerreur
uy= yerreur
#Ici on définit ce qui va nous permettre de faire notre modélisation.
#f est une fonction qui va prendre en entrée le tableau de données en abscisse (ici écrit par x), et les paramètres
#du modèle (ici écrit p), et va ressortir la courbe modèle. Ici on travaille sur un modèle linéaire
#Dx_f est la dérivée du modèle par rapport à x. Ça va servir à calculer la part des incertitudes sur x pour
#l'estimation des paramètres. Ici, on se place toujours dans un cas linéaire
#On définit ensuite un résidu. C'est l'écart au modèle pondéré par l'inverse de cet écart. Les alogrithmes qui
#suivent sur scipy vont minimiser la somme des carrés de ces résidus
def f(x,p):
a,b = p
return a*x+b
def Dx_f(x,p):
a,b = p
return a
def residual(p,y,x):
return (y-f(x,p))/np.sqrt(uy**2 +(Dx_f(x,p)*ux)**2)
#On commence par donner une valeur initiale pour les paramètres, ici p0. Généralement on s'en fout, mais des fois ça
#peut être important pour éviter un minimum local qui gènerait l'optimisation
p0 = np.array([1.5,10])
#Ici on fait tourner la fonction scipy de minimisation des moindres carrées. On lui fait affecter ces résultats
#dans deux variables.
#La première est un n-uplet (ici un couple pour la régression linéaire) qui correspond aux paramètres optimaux.
#La seconde variable donne la matrice de covariance de ces paramètres (ici une matrice 2x2). Il suffit de prendre ses
#coefficients diagonaux pour obtenir les incertitudes sur les paramètres
#(pour les détails techniques, je vous laisse regarder dans la doc de scipy.optimize.leastsq, en gros on inverse
#une approximation de Hessienne obtenue par des matrice de permutations.).
# result = spo.leastsq(residual,p0,args=(y,x), full_output=True)
def f0(x,a,b):
return a*x+b
popt,pcov=spo.curve_fit(f0,x,y,sigma=uy,absolute_sigma=True)
# popt= result[0]
# pcov = result[1]
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
#Ici j'affiche jsute les coeff. Il faut les arrondir en fonction des valeurs d'incertitudes
print('pente = ' + str(round(popt[0],3)) + ' ± ' + str(round(upopt[0],3)))
print('ordonnée à l\'origine = ' + str(round(popt[1],9)) + ' ± ' + str(round(upopt[1],9)))
#Ici on fait juste le tracé. Rien à signaler de particulier. On crée une liste d'abscisse et d'ordonnée pour le
#modèle.Sur cette manip, j'ai eu le problème que les graduations font n'importe quoi, d'où les 10**6 qui se balade
#partout, c'est juste là pour obtenirun graphe élégant.
#Ensuite, on trace en errorbar les points de mesure et le point expérimental (avec une couleur différente)
#Le reste c'est que de la déco pour obtenir un joli graphe (pensez bien à arrondir dans les fonction round() dont
#le dernier argument est la décimale à laquelle il faut arrondir, pour éviter d'avoir 14 chiffres significatifs)
plt.figure(figsize=(12,9))
x_mod = np.linspace(np.amin(x),np.amax(x),1000)
y_mod = f(x_mod,popt)
plt.plot(x_mod,y_mod,label='Modèle affine $y = a x + b$', color = 'blue', linestyle= '-')
plt.plot(x,y)
plt.errorbar(x,y,xerr=ux,yerr=uy,marker='+', color = 'red', linestyle= '',label='Simulation')
#plt.errorbar(x_exp,y_exp,xerr=ux_exp,yerr=uy_exp,marker='+', color = 'green', linestyle= '',label='Point(s) pris en direct')
plt.xlabel('ln(t)',fontsize=25)
plt.ylabel('ln(MSD))',fontsize=25)
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)
#plt.xlim(0.9*np.amin(x), 1.01*np.amax(x))
#plt.ylim(0.9*np.amin(y), 1.01*np.amax(y))
plt.legend(loc='upper left',fontsize=15)
plt.title("Ajustement du MSD pour une probabilité de p="+str(p))
# plt.text((np.amin(x)+np.amax(x))/2.5, np.amin(y) + (np.amin(y) + np.amax(y))/10 , 'Pente a = ' + str(round(popt[0],15)) + ' ± ' + str(round(upopt[0],10)) + ' unité',fontsize=16)
# plt.text((np.amin(x)+np.amax(x))/2.5, np.amin(y) + (np.amin(y) + np.amax(y))/20, 'Ordonnée à l\'origine = ' + str(round(popt[1],10)) + ' ± ' + str(round(upopt[1],10)) + ' unité',fontsize=16)
plt.show()
Séance 18/11
'''Ici, la bactérie décroit la fréquence de rotation quand la oncentration augmente. La réorientation est toujours aléatoire, mais elle change le temps de parcours en pratique. Nous, on a accès à la proba de tourner. Si elle s'approche de nutriments, elle allonge son temps de parcours, donc plutot ligne droite, on va diminuer .' '''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
#%%Définition des paramètres
p=0.1 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=10000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace
x0,y0=l/2,l/2# Conditions initiales
l0=0.001#Pas de longueurs l0
sigma=l/10 # Pour la gausienne, on normalise en l
#%% Profil de concentration, ici gaussien
def C(x,y): #pos est un tableau de coordonnées (x,y)
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
#%% Représentation
x,y=np.linspace(0,l,Ne),np.linspace(0,l,Ne)
M=np.zeros((Ne,Ne))# Position Ne : abscisse l.
for j in range(Ne):
M[:,j]=C(x[:],y[j])
#%% Tracé
plt.figure()
plt.imshow(M)
plt.colorbar(label='Concentration')
plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(x[::Ne//10],2))#On ne met que 10 ticks
plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(y[::-Ne//10],2))# - pour avoir le 0 au bon endroit
plt.title('Profil de concentration dans la boîte : gausienne de largeur l/10',fontsize=18)
plt.xlabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
plt.ylabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
# plt.legend(fontsize=15)
plt.grid(False)
plt.show()
#%% Marche avec concentration
def marche_grad(N):
x,y=np.zeros(N),np.zeros(N)
x[0]=0.4
y[0]=0.4
# P=C(x[0],y[0])#proba d'aller en ligne droite
P=0
angles=np.random.uniform(0,2*np.pi,N)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
for i in range(N-1):
changement=np.random.binomial(1,P)# Si 1, même direction
if changement==0:
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
P=C(x[i+1],y[i+1])# si C est max, on va en ligne droite c'est bon, normalisée à 1
return x,y
#%%
M=marche_grad(Nt)
x,y=M[0],M[1]
#%% Tracé
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot(x[0],y[0],color='r',label='Position initiale',marker='o')
plt.plot([x[Nt-1]],[y[Nt-1]],color='g',label='Position finale',marker='o')
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une concentration gaussienne centrée en l/2,l/2")
plt.show()
#%% MSD
n=50# Nombre de mesures pour moyenner
MSD=np.zeros((Nt,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
for j in range(n):
x,y=marche_grad(Nt)
x,y=x-x0,y-y0
MSD[:,j]=x[:]**2+y[:]**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
#%% Traitement statistique pour tracer MSD
moy=np.zeros(Nt)
erreurs=np.zeros(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot des situations finales
plt.figure()
for i in range(n):
plt.plot(MSD[Nt,i])
#%% Plot
temps=np.array([i for i in range(Nt)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences, p="+str(p))
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
23/11
'''Ici, la bactérie décroit la fréquence de rotation quand la oncentration augmente. La réorientation est toujours aléatoire, mais elle change le temps de parcours en pratique. Nous, on a accès à la proba de tourner. Si elle s'approche de nutriments, elle allonge son temps de parcours, donc plutot ligne droite, on va diminuer .' '''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
#%%Définition des paramètres
p=0.8 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=10000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace
x0,y0=l/2,l/2# Conditions initiales
l0=0.001#Pas de longueurs l0
sigma=l/2 # Pour la gausienne, on normalise en l
#%% Profil de concentration, ici gaussien
def C(x,y): #pos est un tableau de coordonnées (x,y)
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
#%% Représentation
xC,yC=np.linspace(0,l,Ne),np.linspace(0,l,Ne)
MC=np.zeros((Ne,Ne))# Position Ne : abscisse l.
for j in range(Ne):
MC[:,j]=C(xC[:],yC[j])
#%% Tracé
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(xC[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(yC[::-Ne//10],2))# - pour avoir le 0 au bon endroit
plt.title('Profil de concentration dans la boîte : gausienne de largeur l/10',fontsize=18)
plt.xlabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
plt.ylabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
# plt.legend(fontsize=15)
plt.grid(False)
plt.show()
#%% Marche avec concentration
def marche_grad(N):
x,y=np.zeros(N),np.zeros(N)
x[0]=l/2#np.random.uniform(0.2,0.8)
y[0]=l/2#np.random.uniform(0.2,0.8)
# P=C(x[0],y[0])#proba d'aller en ligne droite
P=0
angles=np.random.uniform(0,2*np.pi,N)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
for i in range(N-1):
changement=np.random.binomial(1,P)# Si 1, même direction
if changement==0:
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
# P=0.3
P=C(x[i+1],y[i+1])**0.4# si C est max, on va en ligne droite c'est bon, normalisée à 1
return x,y
#%%
M=marche_grad(Nt)
x,y=M[0],M[1]
#%% Tracé
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot(x[0],y[0],color='r',label='Position initiale',marker='o')
plt.plot([x[Nt-1]],[y[Nt-1]],color='g',label='Position finale',marker='o')
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une concentration gaussienne centrée en l/2,l/2")
plt.show()
#%% Positions finales des bactéries
n=1000# Nombre de mesures pour moyenner
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
X=np.zeros((Nt,n))#n colonnes, Nt lignes. Une colonne correspond à une expérience
Y=np.zeros((Nt,n))
for j in range(n):
x,y=marche_grad(Nt)
X[:,j]=x[:]
Y[:,j]=y[:]
plt.scatter(x[Nt-1],y[Nt-1])
plt.title("Positions finales des bactéries")
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(x[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(y[::-Ne//10],2))# - pour avoir le 0 au bon endroit
# plt.xlim(0,l)
# plt.ylim(0,l)
# plt.aspect("equal")
plt.show()
#%% Histogrammes 2D en positions
posx=X.reshape(1,n*Nt)[0]
posy=Y.reshape(1,n*Nt)[0]
plt.figure()
H=np.histogram2d(posx,posy,bins=100)
plt.imshow(H[0])
plt.show()
#%% MSD
n=1000# Nombre de mesures pour moyenner
MSD=np.zeros((Nt,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
for j in range(n):
x,y=marche_grad(Nt)
plt.scatter(x[Nt-1],y[Nt-1])
x,y=x-x0,y-y0
MSD[:,j]=x[:]**2+y[:]**2
plt.title("Positions finales des bactéries")
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(x[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(y[::-Ne//10],2))# - pour avoir le 0 au bon endroit
# plt.xlim(0,l)
# plt.ylim(0,l)
plt.aspect("equal")
plt.show()
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
#%% Traitement statistique pour tracer MSD
moy=np.zeros(Nt)
erreurs=np.zeros(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot
temps=np.array([i for i in range(Nt-1)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences, p="+str(p))
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
25/11
'''Ici, la bactérie décroit la fréquence de rotation quand la oncentration augmente. La réorientation est toujours aléatoire, mais elle change le temps de parcours en pratique. Nous, on a accès à la proba de tourner. Si elle s'approche de nutriments, elle allonge son temps de parcours, donc plutot ligne droite, on va diminuer .' '''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
import scipy.optimize as spo
plt.rcParams['axes.grid'] = True #Paramètre esthétique, peut être retiré
plt.rcParams['savefig.format'] = 'png' #On peut sélectionner pdf si on aime les images vectorielles
#%%Définition des paramètres
p=0.8 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=500000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace
x0,y0=l/2,l/2# Conditions initiales
l0=0.001#Pas de longueurs l0
sigma=l/2 # Pour la gausienne, on normalise en l
#%% Profil de concentration, ici gaussien
def C(x,y): #pos est un tableau de coordonnées (x,y)
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
#%% Représentation
xC,yC=np.linspace(0,l,Ne),np.linspace(0,l,Ne)
MC=np.zeros((Ne,Ne))# Position Ne : abscisse l.
for j in range(Ne):
MC[:,j]=C(xC[:],yC[j])
#%% Tracé
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(xC[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(yC[::-Ne//10],2))# - pour avoir le 0 au bon endroit
plt.title('Profil de concentration dans la boîte : gausienne de largeur l/10',fontsize=18)
plt.xlabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
plt.ylabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
# plt.legend(fontsize=15)
plt.grid(False)
plt.show()
#%% Marche avec concentration
def marche_grad(N):
x,y=np.zeros(N),np.zeros(N)
x[0]=np.random.uniform(0.2,0.8)
y[0]=np.random.uniform(0.2,0.8)
# P=C(x[0],y[0])#proba d'aller en ligne droite
P=0
angles=np.random.uniform(0,2*np.pi,N)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
for i in range(N-1):
changement=np.random.binomial(1,P)# Si 1, même direction
if changement==0:
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
# P=0.3
P=C(x[i+1],y[i+1])**0.4# si C est max, on va en ligne droite c'est bon, normalisée à 1
return x,y
#%%
M=marche_grad(Nt)
x,y=M[0],M[1]
#%% Tracé
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot(x[0],y[0],color='r',label='Position initiale',marker='o')
plt.plot([x[Nt-1]],[y[Nt-1]],color='g',label='Position finale',marker='o')
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une concentration gaussienne centrée en l/2,l/2")
plt.show()
#%% Positions finales des bactéries
n=100# Nombre de mesures pour moyenner
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
X=np.empty((Nt,n))#n colonnes, Nt lignes. Une colonne correspond à une expérience
Y=np.empty((Nt,n))
for j in range(n):
x,y=marche_grad(Nt)
X[:,j]=x[:]
Y[:,j]=y[:]
plt.plot(x[Nt-1],y[Nt-1],'.')
plt.title("Positions finales des bactéries")
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(x[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(y[::-Ne//10],2))# - pour avoir le 0 au bon endroit
# plt.xlim(0,l)
# plt.ylim(0,l)
# plt.aspect("equal")
plt.show()
#%% Histogrammes 2D en positions
d=2# Plus c'est proche d e1, plus on enlève de positions initiales
Mx=X[Nt//d:Nt-1,:]# On enlève les premières positions de chaque expérience parce qu'on veut le régime stationnaire
My=Y[Nt//d:Nt-1,:]
posx=Mx.flatten()# Concatène les lignes ensemble
posy=My.flatten()
plt.figure()
H=np.histogram2d(posx,posy,bins=100)
plt.pcolormesh(H[1],H[2],H[0])
plt.colorbar()
plt.xlabel("X",fontsize=14)
plt.ylabel('Y',fontsize=14)
plt.title("Histogramme 2D sur les "+str(Nt-Nt//d)+" dernières positions de "+str(n)+" bactéries")
plt.show()
#%% MSD
MSD=np.empty((Nt,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
'''On va reprendre les X,Y calculés juste avant, pas la peine d'en prendre d'autres.' '''
for j in range(n):
MSD[:,j]=(X[:,j]-x0)**2+(Y[:,j]-y0)**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
moy=np.empty(Nt)
erreurs=np.empty(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot
temps=np.array([i for i in range(Nt)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences")
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
''' Ca donne un MSD bizarre, normal. Il faut voir si on peut ajuster la fin par une droite. On réutilise l'ajustement du programme de la particule libre, voir les commentaires sur ce dernier'''
#%% Ajustation du MSD à grand t : on part de N//10 on en aura beaucoup, l'erreur est sur le log, on n'aura pas log(0) en toute logique.
dec=1# combien de décades on prend sur les points
skip=5000# On ne lit qu'un point sur skip pour éviter les biais de l'erreur.
absc=temps[Nt//10**dec:Nt]
ordo=moy[Nt//10**dec:Nt]
yerreur=erreurs[Nt//10**dec:Nt]
'''Copie d'un code déjà utilisé l'année dernière. Ici on trace logMSD en fonction de log(t) pour avoir les pentes et ordonnées origine '''
uy= yerreur/ordo
x = np.log(absc)
y = np.log(ordo)
''' On va enlever quelques termes pour éviter de biaiser l'ajustement'''
y=y[::skip]
x=x[::skip]
uy=uy[::skip]
ux=np.zeros(len(uy))
def f(x,p):
a,b = p
return a*x+b
def Dx_f(x,p):
a,b = p
return a
def residual(p,y,x):
return (y-f(x,p))/np.sqrt(uy**2 +(Dx_f(x,p)*ux)**2)
p0 = np.array([0.1,1])
result = spo.leastsq(residual,p0,args=(y,x), full_output=True)
popt = result[0]
pcov = result[1]
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
x_mod = np.linspace(np.amin(x),np.amax(x),(Nt-Nt//10**dec)//skip)
y_mod = f(x_mod,popt)
#Ici j'affiche jsute les coeff. Il faut les arrondir en fonction des valeurs d'incertitudes
print('pente = ' + str(round(popt[0],8)) + ' ± ' + str(round(upopt[0],8)))
print('ordonnée à l\'origine = ' + str(round(popt[1],4)) + ' ± ' + str(round(upopt[1],4)))
print('Coefficient de corrélation R='+str(np.corrcoef(y,y_mod)[0,1]))
plt.figure(figsize=(12,9))
plt.plot(x_mod,y_mod,label='Modèle affine $y = a x + b$', color = 'blue', linestyle= '-')
plt.plot(x,y,"o")
plt.errorbar(x,y,xerr=ux,yerr=uy,marker='+', color = 'red', linestyle= '',label='Simulation')
plt.xlabel('ln(t)',fontsize=25)
plt.ylabel('ln(MSD))',fontsize=25)
plt.legend(loc='upper left',fontsize=15)
plt.title("Ajustement du MSD pour une probabilité de p="+str(p))
plt.show()
07/12
'''Ici, la bactérie décroit la fréquence de rotation quand la oncentration augmente. La réorientation est toujours aléatoire, mais elle change le temps de parcours en pratique. Nous, on a accès à la proba de tourner. Si elle s'approche de nutriments, elle allonge son temps de parcours, donc plutot ligne droite, on va diminuer .' '''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
import scipy.optimize as spo
plt.rcParams['axes.grid'] = True #Paramètre esthétique, peut être retiré
plt.rcParams['savefig.format'] = 'png' #On peut sélectionner pdf si on aime les images vectorielles
#%%Définition des paramètres
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=50000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace
x0,y0=l/2,l/2# Conditions initiales
l0=l/100#Pas de longueurs l0
sigma=l/2 # Pour la gausienne, on normalise en l
n=500# Nombre de mesures pour moyenner
#%% Profil de concentration, ici gaussien
def C(x,y): #pos est un tableau de coordonnées (x,y)
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
#%% Représentation
xC,yC=np.linspace(0,l,Ne),np.linspace(0,l,Ne)
MC=np.zeros((Ne,Ne))# Position Ne : abscisse l.
for j in range(Ne):
MC[:,j]=C(xC[:],yC[j])
#%% Tracé
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(xC[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(yC[::-Ne//10],2))# - pour avoir le 0 au bon endroit
plt.title('Profil de concentration dans la boîte : gausienne de largeur l/10',fontsize=18)
plt.xlabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
plt.ylabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
# plt.legend(fontsize=15)
plt.grid(False)
plt.show()
#%% Marche avec concentration
def marche_grad(N):
x,y=np.zeros(N),np.zeros(N)
x[0]=np.random.uniform(0.2,0.8)
y[0]=np.random.uniform(0.2,0.8)
# P=C(x[0],y[0])#proba d'aller en ligne droite
P=0.5
angles=np.random.uniform(0,2*np.pi,N)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
for i in range(N-1):
changement=np.random.binomial(1,P)# Si 1, même direction
if changement==0:
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
# P=0.3
P=C(x[i+1],y[i+1])**0.3# si C est max, on va en ligne droite c'est bon, normalisée à 1
return x,y
#%%
M=marche_grad(Nt)
x,y=M[0],M[1]
#%% Tracé
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot(x[0],y[0],color='r',label='Position initiale',marker='o')
plt.plot([x[Nt-1]],[y[Nt-1]],color='g',label='Position finale',marker='o')
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une concentration gaussienne centrée en l/2,l/2")
plt.show()
#%% Création des données de 1000 expériences pour le MSD
X=np.empty((Nt,n))#n colonnes, Nt lignes. Une colonne correspond à une expérience
Y=np.empty((Nt,n))
for j in range(n):
x,y=marche_grad(Nt)
X[:,j]=x[:]
Y[:,j]=y[:]
#%% Histogrammes 2D en positions
d=2# Plus c'est proche d e1, plus on enlève de positions initiales
Mx=X[Nt//d:Nt-1,:]# On enlève les premières positions de chaque expérience parce qu'on veut le régime stationnaire
My=Y[Nt//d:Nt-1,:]
posx=Mx.flatten()# Concatène les lignes ensemble
posy=My.flatten()
plt.figure()
H=np.histogram2d(posx,posy,bins=100)
plt.pcolormesh(H[1],H[2],H[0])
plt.colorbar()
plt.xlabel("X",fontsize=14)
plt.ylabel('Y',fontsize=14)
plt.title("Histogramme 2D sur les "+str(Nt-Nt//d)+" dernières positions de "+str(n)+" bactéries")
plt.show()
#%% MSD
MSD=np.empty((Nt,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
'''On va reprendre les X,Y calculés juste avant, pas la peine d'en prendre d'autres.' '''
for j in range(n):
MSD[:,j]=(X[:,j]-x0)**2+(Y[:,j]-y0)**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
moy=np.empty(Nt)
erreurs=np.empty(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot
temps=np.array([i for i in range(Nt)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences")
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
''' Ca donne un MSD bizarre, normal. Il faut voir si on peut ajuster la fin par une droite. On réutilise l'ajustement du programme de la particule libre, voir les commentaires sur ce dernier'''
#%% Ajustation du MSD à grand t : on part de N//10 on en aura beaucoup, l'erreur est sur le log, on n'aura pas log(0) en toute logique.
dec=1# combien de décades on prend sur les points
skip=1000# On ne lit qu'un point sur skip pour éviter les biais de l'erreur.
absc=temps[Nt//10**dec:Nt]
ordo=moy[Nt//10**dec:Nt]
yerreur=erreurs[Nt//10**dec:Nt]
'''Copie d'un code déjà utilisé l'année dernière. Ici on trace logMSD en fonction de log(t) pour avoir les pentes et ordonnées origine '''
uy= yerreur/ordo
x = np.log(absc)
y = np.log(ordo)
''' On va enlever quelques termes pour éviter de biaiser l'ajustement'''
y=y[::skip]
x=x[::skip]
uy=uy[::skip]
ux=np.zeros(len(uy))
def f(x,p):
a,b = p
return a*x+b
def Dx_f(x,p):
a,b = p
return a
def residual(p,y,x):
return (y-f(x,p))/np.sqrt(uy**2 +(Dx_f(x,p)*ux)**2)
p0 = np.array([0.1,1])
result = spo.leastsq(residual,p0,args=(y,x), full_output=True)
popt = result[0]
pcov = result[1]
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
x_mod = np.linspace(np.amin(x),np.amax(x),(Nt-Nt//10**dec)//skip)
y_mod = f(x_mod,popt)
#Ici j'affiche jsute les coeff. Il faut les arrondir en fonction des valeurs d'incertitudes
print('pente = ' + str(round(popt[0],8)) + ' ± ' + str(round(upopt[0],8)))
print('ordonnée à l\'origine = ' + str(round(popt[1],4)) + ' ± ' + str(round(upopt[1],4)))
print('Coefficient de corrélation R='+str(np.corrcoef(y,y_mod)[0,1]))
plt.figure(figsize=(12,9))
plt.plot(x_mod,y_mod,label='Modèle affine $y = a x + b$', color = 'blue', linestyle= '-')
plt.plot(x,y,"o")
plt.errorbar(x,y,xerr=ux,yerr=uy,marker='+', color = 'red', linestyle= '',label='Simulation')
plt.xlabel('ln(t)',fontsize=25)
plt.ylabel('ln(MSD))',fontsize=25)
plt.legend(loc='upper left',fontsize=15)
plt.title("Ajustement du MSD pour t grand")
plt.show()
#%% Récupération des données : un fichier numpy pour l'utiliser, et un txt pour l'afficher, mais c'est lourd pour les applications.
np.save("positionsX_03_"+str(n)+"_.npy",X)
np.save("positionsY_03_"+str(n)+"_.npy",Y)
# np.savetxt('positionsX.txt', X)
# np.savetxt('positionsY.txt', Y)
#%% Positions radiale en fonction du temps
exp=10
plt.figure()
plt.plot(np.arange(Nt),np.sqrt((X[:,exp]-x0)**2+(Y[:,exp]-y0)**2))
plt.show()
Code pour la conso :
''' Une nouvelle version avec consommation de nutriment prise en compte'''
#%%Modules
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
# import scipy.optimize as spo
#%% Paramètres
p=0.8 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=50000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace : va impacter la zone autour de laquelle la bactérie mange
x0,y0=l/2,l/2# Conditions initiales
l0=l/1000#Pas de longueurs l0
sigma=l/10 # Pour la gausienne, on normalise en l
d=l/Ne
x0,y0=l/2,l/2
gnomgnom=1/100
n=100
extent=0,l,0,l
#%% Fonction mapping et définition de l'espace
def C(x,y):
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
Position=np.linspace(0,l,Ne)
xx,yy=np.meshgrid(Position,Position)
z=C(xx,yy)
plt.figure()
plt.imshow(z,extent=extent)
plt.colorbar(label='concentration')
plt.title("Consommation des bactéries pour n=0 bactéries")
plt.legend()
plt.grid(False)
plt.show()
#%% Fonction d'avancée des bactéries
def marche_conso(N,n): # On met maintenant le nombre de bactéries dans cette fonction, puisqu'elles vont se voir
X=np.empty((N,n))#Les positions de toutes les bactéries en fonction du temps, une colonne correspond à une bactérie. Cette fois, on va modifier ligne par ligne simultannément
Y=np.empty((N,n))
angles=np.random.uniform(0,2*np.pi,(N,n))
X[0,:],Y[0,:]=np.random.uniform(x0-sigma,x0+sigma,n),np.random.uniform(y0-sigma,y0+sigma,n)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
nutriments=z
P=np.zeros(n)
for i in range(N-1):
for j in range(n):
changement=np.random.binomial(1,P[j])# Si 1, même direction
if changement==0:
X[i+1,j]=X[i,j]+pas_x[i,j]
Y[i+1,j]=Y[i,j]+pas_y[i,j]
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i,j]=pas_x[i-1,j]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i,j]=pas_y[i-1,j]
X[i+1,j]=X[i,j]+pas_x[i,j]
Y[i+1,j]=Y[i,j]+pas_y[i,j]
ix=np.int(X[i+1,j]*Ne/l)# Indentations de l'espace, on convertit la position en indice sur la carte de l'espace
iy=np.int(Y[i+1,j]*Ne/l)
if ix<Ne and iy<Ne: # Elle sent le gradient, et on n'aura pas de soucis de définition d'indices
nutriments[ix,iy]=np.abs(nutriments[ix,iy]-gnomgnom)
P[j]=nutriments[ix,iy]**0.5# La probabilité est toujours définie par la concentration, mais mise à jour à cause de la consommation
else : # Sinon, elle sort de la boîte, on la laisse diffuser et si elle revient tant mieux.
P[j]=0
return X,Y,nutriments,pas_x
#%% Test
x,y,nut,px=marche_conso(Nt,n)
#%%
plt.figure()
plt.plot(x[:,2],y[:,2])
M=nut-z
plt.figure()
plt.imshow(nut,extent=extent)
plt.colorbar(label='concentration')
plt.title("Consommation des bactéries pour n="+str(n)+" bactéries")
plt.legend()
plt.grid(False)
plt.show()
09/12
Marche gradient
'''Ici, la bactérie décroit la fréquence de rotation quand la oncentration augmente. La réorientation est toujours aléatoire, mais elle change le temps de parcours en pratique. Nous, on a accès à la proba de tourner. Si elle s'approche de nutriments, elle allonge son temps de parcours, donc plutot ligne droite, on va diminuer .' '''
#%% Def bibliotheques
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
import scipy.optimize as spo
plt.rcParams['axes.grid'] = True #Paramètre esthétique, peut être retiré
plt.rcParams['savefig.format'] = 'png' #On peut sélectionner pdf si on aime les images vectorielles
#%%Définition des paramètres
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=50000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace
x0,y0=l/2,l/2# Conditions initiales
l0=l/100#Pas de longueurs l0
sigma=l/2 # Pour la gausienne, on normalise en l
n=500# Nombre de mesures pour moyenner
#%% Profil de concentration, ici gaussien
def C(x,y): #pos est un tableau de coordonnées (x,y)
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
#%% Représentation
xC,yC=np.linspace(0,l,Ne),np.linspace(0,l,Ne)
MC=np.zeros((Ne,Ne))# Position Ne : abscisse l.
for j in range(Ne):
MC[:,j]=C(xC[:],yC[j])
#%% Tracé
plt.figure()
plt.imshow(MC,extent=[0,l,0,l])
plt.colorbar(label='Concentration')
# plt.xticks([i*Ne//10 for i in range(10)],labels=np.round(xC[::Ne//10],2))#On ne met que 10 ticks
# plt.yticks([i*Ne//10 for i in range(10)],labels=np.round(yC[::-Ne//10],2))# - pour avoir le 0 au bon endroit
plt.title('Profil de concentration dans la boîte : gausienne de largeur l/10',fontsize=18)
plt.xlabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
plt.ylabel('Position dans la boîte (en longueur de boîte)',fontsize=15)
# plt.legend(fontsize=15)
plt.grid(False)
plt.show()
#%% Marche avec concentration
def marche_grad(N):
x,y=np.zeros(N),np.zeros(N)
x[0]=np.random.uniform(0.2,0.8)
y[0]=np.random.uniform(0.2,0.8)
# P=C(x[0],y[0])#proba d'aller en ligne droite
P=0.5
angles=np.random.uniform(0,2*np.pi,N)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
for i in range(N-1):
changement=np.random.binomial(1,P)# Si 1, même direction
if changement==0:
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i]=pas_x[i-1]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i]=pas_y[i-1]
x[i+1]=(x[i]+pas_x[i])
y[i+1]=(y[i]+pas_y[i])
# P=0.3
P=C(x[i+1],y[i+1])**0.5# si C est max, on va en ligne droite c'est bon, normalisée à 1
return x,y
def exp_grad(n,N):# fonction qui crée les données de plusieurs expériences
X=np.empty((N,n))#n colonnes, N lignes. Une colonne correspond à une expérience
Y=np.empty((N,n))
for j in range(n):
x,y=marche_grad(N)
X[:,j]=x[:]
Y[:,j]=y[:]
return X,Y
#%%
M=marche_grad(Nt)
x,y=M[0],M[1]
#%% Tracé
plt.figure()
plt.plot(x,y,label='Trajectoire')
plt.plot(x[0],y[0],color='r',label='Position initiale',marker='o')
plt.plot([x[Nt-1]],[y[Nt-1]],color='g',label='Position finale',marker='o')
plt.legend(loc='best')
plt.title("Trajectoire de la bactérie pour une concentration gaussienne centrée en l/2,l/2")
plt.show()
#%% Création des données de 1000 expériences pour le MSD
X,Y=exp_grad(n,Nt)
#%% Calcul du MSD
MSD=np.empty((Nt,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
'''On va reprendre les X,Y calculés juste avant, pas la peine d'en prendre d'autres.' '''
for j in range(n):
MSD[:,j]=(X[:,j]-x0)**2+(Y[:,j]-y0)**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
moy=np.empty(Nt)
erreurs=np.empty(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot du MSD
temps=np.array([i for i in range(Nt)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences")
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
''' Ca donne un MSD bizarre, normal. Il faut voir si on peut ajuster la fin par une droite. On réutilise l'ajustement du programme de la particule libre, voir les commentaires sur ce dernier'''
#%% Ajustement du MSD à grand t : on part de N//10 on en aura beaucoup, l'erreur est sur le log, on n'aura pas log(0) en toute logique.
dec=1# combien de décades on prend sur les points
skip=1000# On ne lit qu'un point sur skip pour éviter les biais de l'erreur.
absc=temps[Nt//10**dec:Nt]
ordo=moy[Nt//10**dec:Nt]
yerreur=erreurs[Nt//10**dec:Nt]
'''Copie d'un code déjà utilisé l'année dernière. Ici on trace logMSD en fonction de log(t) pour avoir les pentes et ordonnées origine '''
uy= yerreur/ordo
x = np.log(absc)
y = np.log(ordo)
''' On va enlever quelques termes pour éviter de biaiser l'ajustement'''
y=y[::skip]
x=x[::skip]
uy=uy[::skip]
ux=np.zeros(len(uy))
def f(x,p):
a,b = p
return a*x+b
def Dx_f(x,p):
a,b = p
return a
def residual(p,y,x):
return (y-f(x,p))/np.sqrt(uy**2 +(Dx_f(x,p)*ux)**2)
p0 = np.array([0.1,1])
result = spo.leastsq(residual,p0,args=(y,x), full_output=True)
popt = result[0]
pcov = result[1]
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
x_mod = np.linspace(np.amin(x),np.amax(x),(Nt-Nt//10**dec)//skip)
y_mod = f(x_mod,popt)
#Ici j'affiche jsute les coeff. Il faut les arrondir en fonction des valeurs d'incertitudes
print('pente = ' + str(round(popt[0],8)) + ' ± ' + str(round(upopt[0],8)))
print('ordonnée à l\'origine = ' + str(round(popt[1],4)) + ' ± ' + str(round(upopt[1],4)))
print('Coefficient de corrélation R='+str(np.corrcoef(y,y_mod)[0,1]))
plt.figure(figsize=(12,9))
plt.plot(x_mod,y_mod,label='Modèle affine $y = a x + b$', color = 'blue', linestyle= '-')
plt.plot(x,y,"o")
plt.errorbar(x,y,xerr=ux,yerr=uy,marker='+', color = 'red', linestyle= '',label='Simulation')
plt.xlabel('ln(t)',fontsize=25)
plt.ylabel('ln(MSD))',fontsize=25)
plt.legend(loc='upper left',fontsize=15)
plt.title("Ajustement du MSD pour t grand")
plt.show()
#%% Récupération des données : un fichier numpy pour l'utiliser, et un txt pour l'afficher, mais c'est lourd pour les applications.
# np.save("positionsX_test_p05_"+str(n)+".npy",X)
# np.save("positionsY_test_p05_"+str(n)+".npy",Y)
# np.savetxt('positionsX.txt', X)
# np.savetxt('positionsY.txt', Y)
#%% Positions radiale en fonction du temps
exp=10
plt.figure()
plt.plot(np.arange(Nt),np.sqrt((X[:,exp]-x0)**2+(Y[:,exp]-y0)**2))
plt.show()
#%% Histogrammes 2D en positions
d=2# Plus c'est proche d e1, plus on enlève de positions initiales
Mx=X[Nt//d:Nt-1,:]# On enlève les premières positions de chaque expérience parce qu'on veut le régime stationnaire
My=Y[Nt//d:Nt-1,:]
posx=Mx.flatten()# Concatène les lignes ensemble
posy=My.flatten()
plt.figure()
H=np.histogram2d(posx,posy,bins=100)
plt.pcolormesh(H[1],H[2],H[0])
plt.colorbar()
plt.xlabel("X",fontsize=14)
plt.ylabel('Y',fontsize=14)
plt.title("Histogramme 2D sur les "+str(Nt-Nt//d)+" dernières positions de "+str(n)+" bactéries")
plt.show()
Marche gnomgnom
''' Une nouvelle version avec consommation de nutriment prise en compte'''
#%%Modules
import numpy as np
from matplotlib import pyplot as plt
# from tqdm import tqdm
# import scipy.optimize as spo
#%% Paramètres
p=0.8 # Proba de tourner à chaque pas de temps : raideur de la trajectoire (0 rectiligne)
l=1# longueur du côté de la boîte carrée, utile plus tard pour le gradient, va en dépendre.
Nt=50000# Nombre des pas de la simulation temporelle
Ne=1000# Discrétisation espace : va impacter la zone autour de laquelle la bactérie mange
x0,y0=l/2,l/2# Conditions initiales
l0=l/1000#Pas de longueurs l0
sigma=l/10 # Pour la gausienne, on normalise en l
x0,y0=l/2,l/2
gnomgnom=1/100
n=100
extent=0,l,0,l
#%% Fonction mapping et définition de l'espace
def C(x,y):
return np.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
Position=np.linspace(0,l,Ne)
xx,yy=np.meshgrid(Position,Position)
z=C(xx,yy)
plt.figure()
plt.imshow(z,extent=extent)
plt.colorbar(label='concentration')
plt.title("Consommation des bactéries pour n=0 bactéries")
plt.legend()
plt.grid(False)
plt.show()
#%% Fonction d'avancée des bactéries
def marche_conso(N,n,Ne): # On met maintenant le nombre de bactéries dans cette fonction, puisqu'elles vont se voir
X=np.empty((N,n))#Les positions de toutes les bactéries en fonction du temps, une colonne correspond à une bactérie. Cette fois, on va modifier ligne par ligne simultannément
Y=np.empty((N,n))
angles=np.random.uniform(0,2*np.pi,(N,n))
X[0,:],Y[0,:]=np.random.uniform(x0-sigma,x0+sigma,n),np.random.uniform(y0-sigma,y0+sigma,n)
pas_x=np.cos(angles)*l0# en i, Pas que la bactérie va faire au prochain temps i+1
pas_y=np.sin(angles)*l0
Position=np.linspace(x0-l/2,x0+l/2,Ne)
xx,yy=np.meshgrid(Position,Position)
z=C(xx,yy)
nutriments=z
P=np.zeros(n)
for i in range(N-1):
for j in range(n):
changement=np.random.binomial(1,P[j])# Si 1, même direction
if changement==0:
X[i+1,j]=X[i,j]+pas_x[i,j]
Y[i+1,j]=Y[i,j]+pas_y[i,j]
elif changement==1 :# On garde la direction du pas d'avant, mais on calcule à partir de la position précédente quand même
pas_x[i,j]=pas_x[i-1,j]# Très important pour le tour suivant pour garder en mémoire la direction
pas_y[i,j]=pas_y[i-1,j]
X[i+1,j]=X[i,j]+pas_x[i,j]
Y[i+1,j]=Y[i,j]+pas_y[i,j]
ix=np.int(X[i+1,j]*Ne/l)# Indentations de l'espace, on convertit la position en indice sur la carte de l'espace
iy=np.int(Y[i+1,j]*Ne/l)
if ix<Ne and iy<Ne: # Elle sent le gradient, et on n'aura pas de soucis de définition d'indices
nutriments[ix,iy]=np.abs(nutriments[ix,iy]-gnomgnom)
P[j]=nutriments[ix,iy]**0.5# La probabilité est toujours définie par la concentration, mais mise à jour à cause de la consommation
else : # Sinon, elle sort de la boîte, on la laisse diffuser et si elle revient tant mieux.
P[j]=0
return X,Y,nutriments
#%% Test
x,y,nut=marche_conso(Nt,n)
#%%
plt.figure()
plt.plot(x[:,2],y[:,2])
M=nut-z
plt.figure()
plt.imshow(nut,extent=extent)
plt.colorbar(label='concentration')
plt.title("Consommation des bactéries pour n="+str(n)+" bactéries")
plt.legend()
plt.grid(False)
plt.show()
Marche aléatoire : MSD
MSD=np.empty((N,n))# Le tableau avec dans chaque colone l'expérience de diffusion à des t-t0 croissants selon les lignes. On va ainsi pouvoir moyenner et faire des stats pour tracer ln(MSD) en fonction de ln(t-t0)
for j in range(n):
x,y=marche_bactérie(Nt,p)
x,y=x-x0,y-y0
MSD[:,j]=x[:]**2+y[:]**2
''' Ici, on n'a pas encore calculé le MSD car on n'a pas moyenné'''
#%% Traitement statistique pour tracer MSD
moy=np.empty(Nt)
erreurs=np.empty(Nt)
t=1.66# Coeff de student à 95% pour 100 valeurs
for i in range(Nt):
moy[i]=np.mean(MSD[i])# MSD[i] est une ligne contenant les n expériences pour ce temps. On moyenne
erreurs[i]=np.std(MSD[i],ddof=1)
erreurs=t*erreurs/np.sqrt(n)# Pour ne le faire qu'une fois
#%% Plot
temps=np.array([i for i in range(Nt)])
plt.figure()
plt.loglog(temps,moy)
plt.errorbar(temps, moy, yerr=erreurs)
plt.grid()
plt.title("MSD en fonction du temps sur "+str(n)+" expériences, p="+str(p))
plt.xlabel("Temps (en nombre de pas)")
plt.ylabel("MSD")
plt.show()
#%% Ajustement du MSD à grand t : on part de N//10 on en aura beaucoup, l'erreur est sur le log, on naura pas log(0) en toute logique.
dec=2# combien de décades on prend sur les poitns
skip=1000# On ne lit qu'un point sur skip pour éviter les biais de l'erreur.
absc=temps[Nt//10**dec:N:skip]
ordo=moy[Nt//10**dec:N:skip]
yerreur=erreurs[Nt//10**dec:N:skip]# Erreur du log, sera utile.
xerreur=np.zeros((Nt-Nt//10**dec)//skip)
'''Copie d'un code déjà utilisé l'année dernière. Ici on trace logMSD en fonction de log(t) pour avoir les pentes et ordonnées origine '''
x = absc
y = ordo
ux= xerreur
uy= yerreur
#Ici on définit ce qui va nous permettre de faire notre modélisation.
#f est une fonction qui va prendre en entrée le tableau de données en abscisse (ici écrit par x), et les paramètres
#du modèle (ici écrit p), et va ressortir la courbe modèle. Ici on travaille sur un modèle linéaire
#Dx_f est la dérivée du modèle par rapport à x. Ça va servir à calculer la part des incertitudes sur x pour
#l'estimation des paramètres. Ici, on se place toujours dans un cas linéaire
#On définit ensuite un résidu. C'est l'écart au modèle pondéré par l'inverse de cet écart. Les alogrithmes qui
#suivent sur scipy vont minimiser la somme des carrés de ces résidus
def f(x,p):
a,b = p
return a*x+b
def Dx_f(x,p):
a,b = p
return a
def residual(p,y,x):
return (y-f(x,p))/np.sqrt(uy**2 +(Dx_f(x,p)*ux)**2)
#On commence par donner une valeur initiale pour les paramètres, ici p0. Généralement on s'en fout, mais des fois ça
#peut être important pour éviter un minimum local qui gènerait l'optimisation
p0 = np.array([0.001,1])
#Ici on fait tourner la fonction scipy de minimisation des moindres carrées. On lui fait affecter ces résultats
#dans deux variables.
result = spo.leastsq(residual,p0,args=(y,x), full_output=True)
popt = result[0]
pcov = result[1]
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
x_mod = np.linspace(np.amin(x),np.amax(x),(N-N//10**dec)//skip)
y_mod = f(x_mod,popt)
#Ici j'affiche jsute les coeff. Il faut les arrondir en fonction des valeurs d'incertitudes
print('pente = ' + str(round(popt[0],8)) + ' ± ' + str(round(upopt[0],8)))
print('ordonnée à l\'origine = ' + str(round(popt[1],4)) + ' ± ' + str(round(upopt[1],4)))
print('Coefficient de corrélation R='+str(np.corrcoef(y,y_mod)[0,1]))
plt.figure(figsize=(12,9))
plt.plot(x_mod,y_mod,label='Modèle affine $y = a x + b$', color = 'blue', linestyle= '-')
plt.plot(x,y,"o")
plt.errorbar(x,y,xerr=ux,yerr=uy,marker='+', color = 'red', linestyle= '',label='Simulation')
plt.xlabel('ln(t)',fontsize=25)
plt.ylabel('ln(MSD))',fontsize=25)
# plt.xticks(fontsize=25)
# plt.yticks(fontsize=25)
#plt.xlim(0.9*np.amin(x), 1.01*np.amax(x))
#plt.ylim(0.9*np.amin(y), 1.01*np.amax(y))
plt.legend(loc='upper left',fontsize=15)
plt.title("Ajustement du MSD pour une probabilité de p="+str(p))
plt.show()