From 619a9f1f2fda3ab457980e38534bc799738a34e6 Mon Sep 17 00:00:00 2001 From: philcarmona <47579183+philcarmona@users.noreply.github.com> Date: Wed, 27 May 2020 23:08:06 +0200 Subject: [PATCH] Add files via upload --- bdp.py | 3877 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 3877 insertions(+) create mode 100644 bdp.py diff --git a/bdp.py b/bdp.py new file mode 100644 index 0000000..e43460f --- /dev/null +++ b/bdp.py @@ -0,0 +1,3877 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import pandas as pd +from scipy.integrate import odeint,quad +from scipy.stats import kde,beta +#import seaborn +#import mpmath +#from mpmath import mp +#from numba import jit +from importlib import reload +pi=np.pi +from functools import reduce +from scipy.interpolate import CubicSpline +from matplotlib.patches import Rectangle +font = {'family' : 'normal', + 'weight' : 'bold', + 'size' : 22} + +plt.rc('font', **font) +#plt.rc('font', family='serif',size='16') +plt.rc('text', usetex=True) +plt.rc('xtick',labelsize=22) +plt.rc('ytick',labelsize=22) +#plt.rc('axes', labelsize=22) +#plt.rc('axes', titlesize=22) +#couleurcontrole='cyan' +couleurcontrole='#00AFF0' +cgris='lightgrey' + + +def bandegris(a,b,couleur=cgris,alpha=1.0): + plt.axvspan(a, b, color=couleur, alpha=alpha, lw=0) +def axbandegris(ax,a,b,couleur=cgris,alpha=1.0): + ax.axvspan(a, b, color=couleur, alpha=alpha, lw=0) + + +#pour trouver le max d'une fonction à 2 paramètres +#@jit +def tm(f,T,N): + N=int(N) + x=0 + y=1/N + m=f(x,y) + for i in range(N): + for j in range(i+1,N+1): + z=f(i/N,j/N) + if (z>m): + x=i/N + y=j/N + m=z + return((x,y,m)) + +#tm(lambda x,y : y*y-x*x,1.0,5000) + +def rzero(la,mu,T): + lb=quad(la,0,T)[0] + mb=quad(mu,0,T)[0] + return(lb/mb) + + +def lam(x): + return(1.5*(1 +np.sin(2*pi*x))) +def mu(x): + return(1.0) + +rzero(lam,mu,1) + + +def optimrzero(la,mu,rhom,T,N=100): + lb=quad(la,0,T)[0] + mb=quad(mu,0,T)[0] + K= (lb-mb)/rhom + print("K=",K) + m=0.0 + tun,tdeux=0,0 + c=T + h=T/N #le pas + for i in range(N): + for j in range(i+1,N+1): + z=quad(la,i*h,j*h)[0] + if (z>K): + ci=(j-i)*h + print("i,j,ci",i,j,ci) + if(ci < c): + tun,tdeux=i*h,j*h + c=ci + print("tun,tdeux,ci",tun,tdeux,ci) + break + return(tun,tdeux,c) + + +#optimrzero(lam,mu,T=1.0,rhom=0.5) + + + +#la fonction psi +def psi(al): + return(al + (1/pi)*np.sin(pi*al)) + + + +#allons y pour optimiser la proba d'émergence + + +#mardi 27 novembre. Exemple du creneau. Comparaison de deux stratégies. + +def scaling(f,T): + def g(t): + return(f(t/T)) + return(g) +def Scaling(f,T): + def g(t): + return (T*f(t/T)) + return(g) + +def calculepe(lap,Lap,mup,Mup,T): + n=1- np.exp(-T*(Lap(1)-Mup(1))) + def f(t): + a=quad(lambda x : lap(x)*np.exp(-T*(Lap(x)-Mup(x))),0,t)[0] + b=quad(lambda x : lap(x)*np.exp(-T*(Lap(x)-Mup(x))),t,1)[0] + d=np.exp(T*(Lap(t) -Mup(t)))*(b +(1-n)*a)*T + return (n/d) + return(f) #on retourne une fonction t -> pe(tT,T) pour t dans (0,1) +def newcalculepe(lap,Lap,mup,Mup,T,limit=50): + phiun=Lap(1)-Mup(1) + n=1- np.exp(-T*phiun) + def integrand(x): + if (x<1): + return(lap(x)*np.exp(-T*(Lap(x)-Mup(x)))) + else: + return(lap(x-1)*np.exp(-T*(Lap(x-1)-Mup(x-1)+phiun))) + def f(t): + a=quad(integrand,t,t+1,limit=limit)[0] + d=np.exp(T*(Lap(t) -Mup(t)))*a*T + return (n/d) + return(f) #on retourne une fonction t -> pe(tT,T) pour t dans (0,1) + +def nncalculepe(lap,Lap,mup,Mup,T,limit=50): + phiun=Lap(1)-Mup(1) + n=1- np.exp(-T*phiun) + def integrand(x): + if (x<1): + return(mup(x)*np.exp(-T*(Lap(x)-Mup(x)))) + else: + return(mup(x-1)*np.exp(-T*(Lap(x-1)-Mup(x-1)+phiun))) + def f(t): + a=quad(integrand,t,t+1,limit=limit)[0] + d=np.exp(T*(Lap(t) -Mup(t)))*a*T + return (n/(n+d)) + return(f) #on retourne une fonction t -> pe(tT,T) pour t dans (0,1) + +def testcalculepe(gamma=0.3,lzero=1.5,alpha=0.5,T=100): + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(lzero*alpha) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma) + (t-(1-gamma))*lzero*alpha) + + def lam(t,lzero=1.5): + return(lzero*(1+ np.sin(2*pi*t))) + def Lam(t,lzero=1.5): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + + npem=nncalculepe(la,La,mu,Mu,T) + #voirtauxetpe(la,La,mu,Mu,npem,T) + voirtauxetpe(la,La,mu,Mu,npem,T) + +def voirtauxetpe(la,La,mu,Mu,pem,T): + t=np.linspace(0,1,100) + plt.rcParams['figure.figsize'] = [16, 10] + plt.subplot(2,1,1) + plt.title("taux de naissance et de mort") + plt.plot(t,[la(s) for s in t],label=r"$\lambda(t)$") + plt.plot(t,[mu(s) for s in t],label=r"$\mu(t)$") + plt.legend() + + plt.subplot(2,1,2) + phiun=La(1)-Mu(1) + plt.title("probabilite d'extinction, T="+str(T)) + plt.xlabel("instant relatif d'arrivee de l'infecte dans la periode") + plt.plot(t,[pem(s) for s in t],label=r"$p_e(tT,T)$, T="+str(T)) + plt.plot(t,[(La(s)-Mu(s))/phiun for s in t],label=r"$\varphi(t)/\varphi(1)$") + plt.legend() + +#Regardons maintenant ce que donnent les deux strategies de controle +#totu d'abord le cas ou le taux de croissance est au dessus puis en dessous +# du taux de mort. Ici c'est le cas si lzero*alpha < muzero=1 +def prestratb(c,rhom=0.5,gamma=0.3,alpha=0.5,lzero=1.5,T=100): + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(lzero*alpha) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma) + (t-(1-gamma))*lzero*alpha) + + def frho(t): + if (tMu(1)),r"$R_0$ should be greater than $1$" + print("Rzero=",La(1)/Mu(1),"Rzero(rhom)=",rzero) + perho=calculepe(larho,Larho,mu,Mu,T) + t=np.linspace(0,1,100) + plt.title(r"probabilite d'extinction : stratégie $\rho_1$") + plt.xlabel("instant relatif d'arrivee de l'infecte dans la periode") + plt.plot(t,[perho(s) for s in t],label=r"$p_e(tT,T)$") + plt.plot(t,[Larho(s) -Mu(s) for s in t],label=r"$\varphi(\rho)$") + plt.legend() + pemoyennerho=(quad(perho,0,1)[0]) + print("pemoyenne=",pemoyennerho) + tstar=(La(1)-Mu(1))/(lzero-1) + rhomb=rhom*(1-gamma)/tstar + def frhob(t): + if (tMu(1)),r"$R_0$ should be greater than $1$" + print("Rzero=",La(1)/Mu(1),"Rzero(rhom)=",rzero) + voirprobext(larho,Larho,mu,Mu,T,nom=r"stratégie $\rho_1$") + + tstar=(La(1)-Mu(1))/(lzero-1) + # rhomb=rhom*(1-gamma)/tstar + # def frhob(t): + # if (t0) and (cout < phiun) + assert (tun <= tdeux) and (tdeux <= 1-gamma) + assert ((cout/lzero) + deltat <= 1-gamma) + phirhoun=phiun-cout + rhosm=(cout/(lzero*(1-gamma-deltat))) + phirhotun=tun*(lzero*(1-rhosm)-1) + phirhotdeux=phirhotun+deltat*(lzero-1) + print("phiroun,rhosm,phirhotun,phirhotdeux=",phirhoun,rhosm,phirhotun,phirhotdeux) + if (phirhotun<=0): + return(phirhoun/lzero) + elif (phirhoun < phirhotun): + return(phurhoun/(lzero*(1-rhosm))) + elif (phirhoun < phirhotdeux): + return((phirhoun+phirhotun*(rhosm/(1-rhosm)))/lzero) + else: + return((phirhoun+(phirhotdeux-phirhotun)/(lzero*(1-rhosm)-1))/(lzero*(1-rhosm))) + + +def illutheoeff(cout=0.2,gamma=0.3,lzero=2,tun=0,tdeux=0.5,T=100): + etheo=theoeff(cout=cout,gamma=gamma,lzero=lzero,tun=tun,tdeux=tdeux) + #c=1-gamma + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + + phiun=lzero*(1-gamma) -1 + deltat=tdeux-tun + assert (phiun>0) and (cout < phiun) + assert (tun <= tdeux) and (tdeux <= 1-gamma) + assert ((cout/lzero) + deltat <= 1-gamma) + phirhoun=phiun-cout + rhosm=(cout/lzero)/(1-gamma-deltat) + + def frho(t): + if ((t< tun) or (tdeux< t< 1-gamma)): + return(1-rhosm) + else: + return(1.0) + def larho(t): + return(la(t)*frho(t)) + def Intlarho(t): + if (t < tun): + return(t*lzero*(rhosm)) + elif (t< tdeux): + return(tun*lzero*(rhosm)) + elif (t< 1-gamma): + return((tun+(t-tdeux))*lzero*(rhosm)) + else: + return((tun+(1-gamma-tdeux))*lzero*(rhosm)) + + def Larho(t): + return(La(t) -Intlarho(t)) + print(Larho(1),Intlarho(1)) + print("phirhotun,phirhotdeux,phirhoun",Larho(tun)-Mu(tun),Larho(tdeux)-Mu(tdeux),Larho(1)-Mu(1)) + print("pemoyennetheorique=",etheo) + voirprobext(larho,Larho,mu,Mu,T=T) + plt.figure() + t=np.arange(0,1,1/100) + plt.plot(t,[larho(s) for s in t],label=r"$\lambda(t)$") + plt.legend() + + +#vendredi 14 decembre 2018. Comparaison de la strategie naive et de la strategie optimale +def compstrat(lzero=2.0,gamma=0.3,beamer=False): + phiun=lzero*(1-gamma)-1 + assert (phiun>0) + t=np.arange(0,phiun,1/100)#les valeurs prises par le cout + #print(phiun,t) + def couopt(c): + return((phiun-c)/lzero) + def counaif(c): + return(couopt(c)/(1- c/(lzero*(1-gamma)))) + if (beamer): + plt.rcParams['figure.figsize'] = [14, 10] + else: + plt.rcParams['figure.figsize'] = [10,14] + + #plt.title("mean emergence probability for the step example :"+ r"$\lambda_0=$"+str(lzero)+r", $\gamma=$"+str(gamma)) + plt.xlabel("cost") + plt.ylabel("mean emergence probability",fontsize=16) + plt.plot(t,[couopt(s) for s in t],dashes=[2,2],color='black', label=r"optimal") + plt.plot(t,[counaif(s) for s in t], color='black', label=r"naive") + plt.legend() + + plt.savefig("compstratstep.pdf",dpi=300) + + +def ncompstrat(lzero=2.0,rzero=3.0,beamer=False): + assert (rzero>1) + cmax=(rzero -1)/lzero + #t=np.arange(0,cmax,1/100)#les valeurs prises par le cout + t=np.linspace(0,cmax,100)#les valeurs prises par le cout + #print(phiun,t) + def couopt(c): + return((rzero -1 - c*lzero)/lzero) + def counaif(c): + return(couopt(c)*rzero/(rzero - c*lzero)) + if (beamer): + plt.rcParams['figure.figsize'] = [16, 10] + else: + plt.rcParams['figure.figsize'] = [10,16] + + #plt.title("mean emergence probability for the step example :"+ r"$\lambda_0=$"+str(lzero)+r", $\gamma=$"+str(gamma)) + plt.xlabel(r"\bf{Cost of control,} ${C}$",fontsize=30) + plt.ylabel(r"\bf{Mean emergence probability} ${}$",fontsize=30) + plt.yticks((0.0, 0.2,0.4,0.6,0.8), (r'\bf{0}', r'\bf{0.2}', r'\bf{0.4}',r'\bf{0.6}',r'\bf{0.8}'), color='k', size=20) + plt.xticks((0.0, 0.2,0.4,0.6,0.8,1.0), (r'\bf{0}', r'\bf{0.2}', r'\bf{0.4}',r'\bf{0.6}',r'\bf{0.8}',r'\bf{1.0}'), color='k', size=20) + plt.plot(t,[couopt(s) for s in t],dashes=[2,2],color='black', label=r"optimal") + plt.plot(t,[counaif(s) for s in t], color='black', label=r"naive") + aprest=np.arange(cmax,1.5*cmax,1/100) + plt.plot(aprest,[0.0 for s in aprest],color='black') + plt.text(cmax+0.01,0.05,r"${\frac{R_0 -1}{\lambda_0}}$",fontsize=30) + + #plt.legend() + + plt.savefig("ncompstratstep.pdf",dpi=300) + + + + +#Lundi 17 decembre essayons de calculer numeriquement + +def casin(lzero=2.0,T=100,C=0.2): + def la(t): + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def La(t): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + + def sinqextmoy(tun,tdeux,C=C,voir=False): + assert (tun >=0) + if (tun >= 0.5): + return(0.0) + if (tdeux <= tun): + return(0.0) + deltat=tdeux -tun + if (deltat < C): + return(0.0) + rhom=C/deltat + def larho(t): + f=(1-rhom) if ((tun=0) + #assert (rhom>0) and (rhom<=1) + deltat=min(C/rhom,1-tun) + tdeux=tun+deltat + def larho(t): + f=(1-rhom) if ((tun$="+'{:.4f}'.format(m) + plt.title(titre) + + + + + +def simusin(lzero=2.0,C=0.2): + def la(t): + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def La(t): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + + def sinqextmoy(tun,rhom,C=C,voir=False): + deltat=min(C/rhom,1-tun) + tdeux=tun+deltat + def larho(t): + f=(1-rhom) if ((tun$=""" +'{:.4f}'.format(m)+r""" for $t_1$=""" + str(xopt)+r""", $\rho_M=$"""+str(yopt) + ax.set_title(titre) + plt.show() + + +def grillerzerosin(nbpt,C=0.2,tunmax=0.5): + f=rzerosin(C=C) + tun=np.arange(0,tunmax,tunmax/nbpt) + rho=(1.0/nbpt)+np.arange(0,1.0,1.0/nbpt) + zs = np.array([[f(i,j) for j in rho] for i in tun]) + return(zs) + +def testrzerosin(nbpt=100,C=0.2): + z=grillerzerosin(nbpt=nbpt,C=C) + np.save("testrzerosinN"+str(nbpt)+"C"+str(C),z) + return(z) + +def rzerosinvoir(z,C=0.2,tunmax=0.5): + N=(z.shape)[0] + tun=np.arange(0.0,tunmax,tunmax/N) + X,Y=np.meshgrid((1.0/N)+np.arange(0.0,1.0,1/N),tun) + fig, ax = plt.subplots() + levels = np.arange(1.6, 1.8, 0.001) + CS = ax.contour(X,Y,z,levels=levels,cmap='flag') + #CS = ax.contour(X,Y,z,cmap='flag') + ax.clabel(CS, inline=1, fontsize=10) + ax.set_ylabel(r"$t_1$") + ax.set_xlabel(r"$\rho_M$") + m=z.min() + y=np.where(z==m) + xopt=(y[0][0]/N)*tunmax + yopt=y[1][0]/N + titre=r"Basic Reproduction Number: Contour Plot. Cost="+str(C)+"\n"+ r""" minimal $R_0$=""" +'{:.4f}'.format(m)+r""" for $t_1$=""" + str(xopt)+r""", $\rho_M=$"""+str(yopt) + ax.set_title(titre) + plt.show() + +def rzerosin(lzero=2.0,C=0.2): + def la(t): + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def La(t): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + + def rzerorho(tun,rhom,C=C): + deltat=min(C/rhom,1-tun) + tdeux=tun+deltat + return((La(1) -rhom*(La(tdeux)-La(tun)))/Mu(1)) + + return(rzerorho) + +########## pour le cas step, on visualise a la fois plusieurs periodes +def stepvoirphietpe(gamma=0.3,lzero=4.0,lT=[20,50,100]): + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + + t=np.linspace(0,1,100) + #plt.rcParams['figure.figsize'] = [16, 10] + plt.subplot(2,1,1) + plt.title("Birth and Death rates") + plt.plot(t,[la(s) for s in t],label=r"$\lambda(t)$") + plt.plot(t,[mu(s) for s in t],label=r"$\mu(t)$") + plt.legend() + + plt.subplot(2,1,2) + plt.xlabel("Introduction time of the infected") + + phiun=La(1)-Mu(1) + plt.title("Extinction probability") + for T in lT: + pem=nncalculepe(la,La,mu,Mu,T) + plt.plot(t,[pem(s) for s in t],label=r"$p_e(tT,T)$, T="+str(T)) + plt.plot(t,[(La(s)-Mu(s))/phiun for s in t],label=r"$\varphi(t)/\varphi(1)$") + plt.legend() + + +################ jeudi 24 janvier 2019 +#from scipy.integrate import odeint +#tout d'abord ecrire la fonction a optimiser +def lafoncao(l12,mu1,l21,mu2,T,nbpts=200): + def msisi(x,t): + return([(l12(t)+mu1(t))*x[0]-(mu1(t)+l12(t)*x[0]*x[1]), + (l21(t)+mu2(t))*x[1]-(mu2(t)+l21(t)*x[0]*x[1])]) + def f(y): + timeint=np.arange(0,T+1/nbpts,T/nbpts) + z=np.array(odeint(msisi,y,timeint)) + fin=z[-1] + print("y1,y2=",y," z=",z) + plt.plot(timeint,z[:,0],label="y1") + plt.plot(timeint,z[:,1],label="y2") + plt.legend() + r=fin-y + print("r=",r,"norme(r)carree=",(r*r).sum()) + return((r*r).sum()) + return(f) + +def l21(x): + return(1/8) +def l12(x): + return(0.22) +def mu1(x): + return(1/5) +def mu2(x): + return(1/8) +f=lafoncao(l12,mu1,l21,mu2,1) + +#la solution stationnaire, quand les fonctions sont constantes +def sstat(l12,mu1, l21,mu2): + x=1 + return( [1- (l12(x)*l21(x)-mu1(x)*mu2(x))/(l21(x)*(l12(x)+mu1(x))), + 1- (l12(x)*l21(x)-mu1(x)*mu2(x))/(l12(x)*(l21(x)+mu2(x)))]) + + +def champvec(l12,mu1, l21,mu2): + def f(x): + t=1 + return(np.array([(l12(t)+mu1(t))*x[0]-(mu1(t)+l12(t)*x[0]*x[1]), + (l21(t)+mu2(t))*x[1]-(mu2(t)+l21(t)*x[0]*x[1])])) + return(f) +y=sstat(l12,mu1, l21,mu2) +F=champvec(l12,mu1, l21,mu2) +F([1,1]),F(y) #les deux doivent donner 0 +#revenons au test +#f([0.2,0.3]) + + +#on observe meme pour des coeff constants un systeme de lotka volterra avec un equilibre instable + +#vendredi 25 janvier. Utilisons l'approche de Bacaer +def oldfbaca(l12,mu1,l21,mu2,T,tau,nbpts=200): + def msisi(x,t): + return([(l12(t)+mu1(t))*x[0]-(mu1(t)+l12(t)*x[0]*x[1]), + (l21(t)+mu2(t))*x[1]-(mu2(t)+l21(t)*x[0]*x[1])]) + def f(y): + timeint=np.arange(0,T+1/nbpts,T/nbpts) + z=np.array(odeint(msisi,y,timeint)) + fin=z[-1] + print("y1,y2=",y," z=",z) + plt.plot(timeint,z[:,0],label="y1") + plt.plot(timeint,z[:,1],label="y2") + plt.legend() + r=fin-y + print("r=",r,"norme(r)carree=",(r*r).sum()) + return((r*r).sum()) + return(f) + + +#lundi 28 janvier, on ecrit carrement les equations verifiees par z +#on resout l'equation (8) D Bacaer +def fbaca(l12,mu1,l21,mu2,tau,T,nbpts=200,zzero=[1.0,1.0],tv=False): + def msisi(x,t): + s=(tau-t)/T + return([-mu1(s)*x[0] +l12(s)*x[1]*(1-x[0]), + -mu2(s)*x[1] +l21(s)*x[0]*(1-x[1])]) + + timeint=np.arange(0,tau+1/nbpts,tau/nbpts) + z=np.array(odeint(msisi,zzero,timeint)) + #print("T=",T,"Z.shape",z.shape,"z=",z) + i=int(T*nbpts/tau) + res=z[-i:] + if tv: + t=np.arange(1/nbpts,T+1/nbpts,T/nbpts) + plt.plot(timeint,z,label="z") + return(res[::-1]) +def fvraipe(l12,mu1,l21,mu2,T,nbpts=200,zzero=[1.0,1.0]): + def msisi(x,t): + s=t/T + return([mu1(s)*x[0] -l12(s)*x[1]*(1-x[0]), + mu2(s)*x[1] -l21(s)*x[0]*(1-x[1])]) + + timeint=np.arange(0,T+1/nbpts,T/nbpts) + z=np.array(odeint(msisi,zzero,timeint)) + return(z) + + #la solution stationnaire, quand les fonctions sont constantes, pour les z +def zstat(l12,mu1, l21,mu2): + x=1 + return( [(l12(x)*l21(x)-mu1(x)*mu2(x))/(l21(x)*(l12(x)+mu1(x))), + (l12(x)*l21(x)-mu1(x)*mu2(x))/(l12(x)*(l21(x)+mu2(x)))]) + +def zchampvec(l12,mu1, l21,mu2): + def f(x): + t=1 + tau=100 + return(np.array([-mu1(tau-t)*x[0] +l12(tau-t)*x[1]*(1-x[0]), + -mu2(tau-t)*x[1] +l21(tau-t)*x[0]*(1-x[1])])) + return(f) + +def testonsfbacaetfvraipe(): + #on a bien convergence vers la solution constante, quand les coeff le sont + z=fbaca(l12,mu1,l21,mu2,tau=1000,T=100,nbpts=10000) + fin=z[-1] + g=zchampvec(l12,mu1, l21,mu2) + vraisol=zstat(l12,mu1, l21,mu2) + fin,g(fin),vraisol,g(vraisol) #presque 0.0, on est proche de la solution + #et si on lance l'ode depuis cette solution approchee cela pose probleme + zz=fvraipe(l12,mu1,l21,mu2,T=100,nbpts=200,zzero=fin) + zzz=fvraipe(l12,mu1,l21,mu2,T=100,nbpts=200,zzero=vraisol) + plt.plot(zz) + plt.plot(zzz) + + #alors que si on lance l'ode deouis un autre point cela explose + z4=fvraipe(l12,mu1,l21,mu2,T=100,nbpts=200,zzero=[0.1,0.2]) + +#c'est normal si on calcule le Jacobien au point d'équilibre on trouve la trac et le determinant >0 donc le systeme est hyperbolique : le point fixe est instable. +#introuduisons maintenant un coeff periodique + +def tperfbaca(ldu=1,ab=2.0,ep=1,muu=1,mud=1,tau=300,T=50,nb=10): + def l21(x): + return(ldu) + def l12(x): + return(ab*(1+ep*np.sin(2*np.pi*x))) + def mu1(x): + return(muu) + def mu2(x): + return(mud) + z=fbaca(l12,mu1,l21,mu2,tau=tau,T=T,nbpts=tau*nb,tv=True) + +#tperfbaca(ab=2.0) + +#comparons avec l'approximation fast/slow +def fslow(l12,mu1,l21,mu2,T,limit=50): + def mu(x): + return(mu1(x)) + def Mu(x): + return(x *mu1(x)) #seule la fonctionl12 n'est pas constante + def la(x): + return(l21(x)*l12(x)/mu2(x)) + def La(x): + return(quad(la,0,x,limit=limit)[0]) + pem=nncalculepe(la,La,mu,Mu,T) + t=np.linspace(0,1,100) + return(np.array([pem(s) for s in t])) + +#introuduisons maintenant un coeff periodique +def pexper(ldu=1/8,ab=0.22,ep=0.3,muu=1/25,mud=1/8,tau=300,T=50,nb=10): + def l21(x): + return(ldu) + def l12(x): + return(ab*(1+ep*np.sin(2*np.pi*x))) + def mu1(x): + return(muu) + def mu2(x): + return(mud) + def cible1(x): + z=(l12(x)*l21(x)-mu1(x)*mu2(x))/(l21(x)*(l12(x)+mu1(x))) + return(z if (z>=0) else 0.0) + z=fbaca(l12,mu1,l21,mu2,tau=tau,T=T,nbpts=tau*nb) + plt.plot(np.linspace(0,1,z.shape[0]),z[:,0],label="one infected human") + plt.plot(np.linspace(0,1,z.shape[0]),z[:,1],label="one infected vector") + #zz=fslow(l12,mu1,l21,mu2,T=T,limit=50) + #plt.plot(np.linspace(0,1,zz.shape[0]),zz,label="fast/slow approximation (human)") + tt=np.linspace(0,1,100) + tz=np.array([cible1(x) for x in tt]) + plt.plot(tt,tz,label="asymptotic profile for one infected human") + plt.title("Emergence probability") + plt.xlabel("time of introduction") + plt.legend() + + + + +#mettons les memes parametres en moyenne excepte pour le parametre oscillant +def testonspexper(): + pexper(muu=1,ldu=1,ab=2.0,ep=1,mud=1,tau=500,T=100) #ca colle rzero=2 + + plt.figure() + pexper(muu=1,ldu=1,ab=1.2,ep=1,mud=1,tau=500,T=100) #ca ne arche pas car le rzero du system BD de dimension 2 est inferieur a 1, alors que le Rzero du systeme 1d BD approximant est superieueur a 1, ici rzero=1.2 (en prenant les valeurs moyennes : je ne sais pas l'approcher comme le fait Bacaer , je prends la formule du rzero pour coeff constants et je remplace les coeff par leur valeur moyenne). + + + +#on s'appercoit que l'on peut avoir des probabilites d'emergence quasi nulles avec ab ldu - mud muu >0 : ce n'est donc pas le bon critère pour avoir Rzero>1 +#on a un decalage certain avec l'approximation fast/slow + +#il reste a verifier que la fonction proba d'emergence periodique donnee par la'lgorithme de Bacaer est effectivement une solution periodique du systeme differentiel calcule pour les proba d'emergence +def checkbaca(ldu=1,ab=2.0,ep=1.0,muu=1,mud=1,tau=300,T=50,nb=10): + def l21(x): + return(ldu) + def l12(x): + return(ab*(1+ep*np.sin(2*np.pi*x))) + def mu1(x): + return(muu) + def mu2(x): + return(mud) + z=fbaca(l12,mu1,l21,mu2,tau=tau,T=T,nbpts=tau*nb) + plt.plot(np.linspace(0,1,z.shape[0]),z,label="2d approximation") + zz=fvraipe(l12,mu1,l21,mu2,T=T,nbpts=T*nb,zzero=z[0]) + plt.plot(np.linspace(0,1,zz.shape[0]),zz,label="solution initialisee") + plt.legend() + return(z,zz) +#z,zz=checkbaca(ldu=1,ab=2.0,ep=1.0,muu=1,mud=1,tau=300,T=50,nb=10) +#a la fin de la periode cela ne correspond plus mais je ne vois pas pourquoi. + + + + +############################################################################# +###### jeudi 31 janvier 2019 +########## generation des figures complexes pour le papier +#import seaborn as sns; +#sns.set() +#sns.set_color_codes() + +def add_arrow(line, position=None, direction='right', size=15, color=None): + """ + add an arrow to a line. + + line: Line2D object + position: x-position of the arrow. If None, mean of xdata is taken + direction: 'left' or 'right' + size: size of the arrow in fontsize points + color: if None, line color is taken. + """ + if color is None: + color = line.get_color() + + xdata = line.get_xdata() + ydata = line.get_ydata() + + if position is None: + position = xdata.mean() + # find closest index + start_ind = np.argmin(np.absolute(xdata - position)) + if direction == 'right': + end_ind = start_ind + 1 + else: + end_ind = start_ind - 1 + + line.axes.annotate('', + xytext=(xdata[start_ind], ydata[start_ind]), + xy=(xdata[end_ind], ydata[end_ind]), + arrowprops=dict(arrowstyle="->", color=color), + size=size + ) + +def trouvelsinmu(lz,of): + if ((of>=1) or (of< 1-2*lz)): + return((False,0.0,0.0)) + else: + a=0.5 - (np.arcsin(((1-of)/lz) -1))/(2*np.pi) + return((True,a,1.5-a)) + +def periodise(f): + #retourne la fonction qui etait definie sur [0,1] periodisee sur [0,2] + def g(t): + return(f(t-np.floor(t))) + return(g) +def periodiseetintegre(f): + + def g(t): + k=np.floor(t) + return(f(t-k)+k*f(1)) + return(g) + +def tauxphietpe(gamma=0.3,lzero=3.0,lT=[0.2,1,5,20,50,100],offsetlam=0.0): + def flechewinter(hauteur,xmin,xmax,couleur=cgris): + hdelta=0.05 + #plt.arrow(xmax-hdelta,hauteur,-(xmax-xmin)+2*hdelta,0,shape='full', length_includes_head=True,color=couleur,alpha=1.0,linewidth=15,head_starts_at_zero=True,head_width=5e-2) + plt.arrow(xmax,hauteur,-(xmax-xmin)+2*hdelta,0,shape='full', length_includes_head=True, + color=couleur,alpha=1.0,linewidth=15,head_starts_at_zero=True,head_width=5e-2) + def grise(): + bandegris(1-gamma, 1.0) + bandegris(1+1-gamma,2.0) + #atracer,asin,bsin=bdp.trouvelsinmu(lzero*(1-gamma),0.0) + def axsingrise(ax): + if (atracer): + axbandegris(ax,asin,bsin) + axbandegris(ax,asin+1,bsin+1) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero+offsetlam) + else: + return(0.0+offsetlam) + def La(t): + if (t< 1-gamma): + return((lzero+offsetlam)*t) + else: + return(lzero*(1-gamma)+offsetlam*t) + def phi(t): + return(La(t) -Mu(t)) + #il faut ajuster le lzerosin pour que les processus aient le meme rzero + + def lasin(t): + return(0.5*lzero*(1+np.sin(2*pi*t))+offsetlam) + def Lasin(t): + return(0.5*lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))+t*offsetlam) + def phisin(t): + return(Lasin(t) -Mu(t)) + + print("phi(1),phisin(1)",phi(1),phisin(1)) + + pla=periodise(la) + pmu=periodise(mu) + plasin=periodise(lasin) + + et=np.linspace(0,2,200) + t=np.linspace(0,1,100) + + pphi=periodiseetintegre(phi) + pphisin=periodiseetintegre(phisin) + vpphi=np.array([pphi(s) for s in et]) + vpphisin=np.array([phisin(s) for s in et]) + vphimax=max(vpphi.max(),vpphisin.max()) + + nblignes=3 + f,ax=plt.subplots(nrows=nblignes, ncols=2, figsize=[18, 16]) #pour les graphiques beamer + #f,ax=plt.subplots(nrows=3, ncols=2, figsize=[11, 16]) #pour les graphiques A4 + #plt.tight_layout(h_pad=3.5) + plt.subplots_adjust(left=0.125,right=0.9,bottom=0.1,top=0.9,hspace=0.5) + + #on trace les taux de naissance et de mort step + #plt.subplot(nblignes,2,1) + axc=ax[0,0] + axabsticks(axc) + mettrelettreax(axc,'A',gauche=0.0) + if (offsetlam<1): + axbandegris(axc,1-gamma, 1.0) + axbandegris(axc,1+1-gamma,2.0) + axc.set_ylabel(r"\bf{Birth and Death rates}") + axc.set_yticks((0,1,2,3)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}",r"\bf{3}")) + vpla=np.array([pla(s) for s in et]) + vplasin=np.array([plasin(s) for s in et]) + vmax=max(vpla.max(),vplasin.max()) + axc.plot(et,vpla,label=r"$\lambda(t)$",color='black') + axc.plot(et,[pmu(s) for s in et],label=r"$\mu(t)$",color='black',dashes=[6,2]) + axc.axis([0.0, 2.0,0.0, 1.3*vmax]) #par defaut dans matplotlib sans seaborn + axc.text(0.05,1.1,r"$\mu(t)$",fontsize=30) + axc.text(0.05,la(0.0)+0.1,r"$\lambda(t)$",fontsize=30) + + #on trace les taux de naissance et de mort sinusoidaux + #plt.subplot(nblignes,2,2) + axc=ax[0,1] + axabsticks(axc) + mettrelettreax(axc,'B') + axc.set_yticks((0,1,2,3)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}",r"\bf{3}")) + atracer,asin,bsin=trouvelsinmu(lzero*(1-gamma),offsetlam) + axsingrise(axc) + axc.plot(et,vplasin,label=r"$\lambda(t)$",color='black') + axc.plot(et,[pmu(s) for s in et],label=r"$\mu(t)$",color='black',dashes=[6,2]) + axc.axis([0, 2.0, 0.0, 1.3*vmax]) #par defaut dans matplotlib sans seaborn + + ########################################## + ##### courbe phi dans le cas step + ################################## + axc=ax[1,0] + axabsticks(axc) + mettrelettreax(axc,'C',gauche=0.0) + + if (offsetlam<1): + axbandegris(axc,1-gamma, 1.0) + axbandegris(axc,1+1-gamma,2.0) + + axc.set_ylabel(r"\bf{Integrated growth rate} $\varphi(t)$") + axc.set_yticks((0.0,1.0,2.0)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}")) + axc.axis([0, 2.0, 0.0, 1.1*vphimax]) #par defaut dans matplotlib sans seaborn + axc.plot(et,vpphi,label=r"$\varphi(t)$",color='black') + + tstar=(phi(1))/(lzero-1) + minloc=phi(1) + print("tstar=",tstar) + axc.text(tstar-0.13,0.02,r"${t^*}$",fontsize=30) + abst=np.linspace(tstar,1,100) + + line=axc.plot(abst,np.full_like(a=abst,fill_value=minloc),dashes=[6,2],color='black')[0] + + #add_arrow(line,position=0,direction='left',size=30,color='red') + h=0.1 + axc.arrow(tstar+h,minloc,-h,0,shape='full', lw=0, length_includes_head=True, head_width=.05,color='black') + + ybst=np.linspace(0,minloc,100) + + axc.plot(np.full_like(a=ybst,fill_value=tstar),ybst,dashes=[6,2],color='black') + axc.arrow(tstar,h,0,-h,shape='full', lw=0, length_includes_head=True, head_width=.05,color='black') + axbandegris(axc,tstar,1-gamma,alpha=0.5) + axbandegris(axc,1+tstar,1+1-gamma,alpha=0.5) + + ################################################################ + # on trace la courbe phi dans le cas sinusoidal + axc=ax[1,1] + axabsticks(axc) + mettrelettreax(axc,'D') + axc.set_yticks((0.0,1.0,2.0)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}")) + axc.axis([0, 2.0, 0.0, 1.1*vphimax]) #par defaut dans matplotlib sans seaborn + + axc.plot(et,vpphisin,label=r"$\varphi(t)$",color='black') + indl=np.linspace(0.5,1,1000) + x=np.array([phisin(s) for s in indl]) + w=np.where(x==x.min()) + tminloc=indl[w[0][0]] + print("tminloc=",tminloc) + minloc=phisin(tminloc) + #maintenant on recherche tstar + indt=np.linspace(0.0,0.5,1000) + x=np.array([phisin(s) for s in indt]) + w=np.where(np.abs(x-minloc)<= 1e-3) + tstarsin=indt[w[0][0]] + print("tstarsin=",tstarsin,"phisin(tstar)",phisin(tstarsin),"phisin(tminloc)",phisin(tminloc)) + axc.text(tstarsin-0.13,0.02,r"${t^*}$",fontsize=30) + abst=np.linspace(tstarsin,tminloc,100) + + axc.plot(abst,np.full_like(a=abst,fill_value=minloc),dashes=[6,2],color='black') + + #la fleche pointillee depuis le minimum local jusque a tstar + h=0.1 + axc.arrow(tstarsin+h,minloc,-h,0,shape='full', lw=0, length_includes_head=True, head_width=.05,color='black') + + ybst=np.linspace(0,minloc,100) + + axc.plot(np.full_like(a=ybst,fill_value=tstarsin),ybst,dashes=[6,2],color='black') + axc.arrow(tstarsin,h,0,-h,shape='full', lw=0, length_includes_head=True, head_width=.05,color='black') + + def axsinwicetwinter(ax): + axsingrise(axc) + axbandegris(axc,tstarsin,asin,alpha=0.5) + axbandegris(axc,1+tstarsin,1+asin,alpha=0.5) + return() + axsinwicetwinter(axc) + + #### proba emergence cas sinusoidal + #plt.subplot(nblignes,2,6) + axc=ax[2,1] + mettrelettreax(axc,'F') + axabsticks(axc) + axc.set_yticks((0.0,0.2,0.4,0.6)) + axc.set_yticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}")) + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + vpemax=0.0 + cs=CubicSpline(np.array([0.0,10,100]),np.array([0.2,0.5,0.9])) + + for T in lT: + pem=nncalculepe(lasin,Lasin,mu,Mu,T) + ppem=periodise(pem) + valppem=np.array([ppem(s) for s in et]) + axc.plot(et,valppem,label=r"T="+str(T),color=plt.cm.YlOrRd(cs(T))) + vpemax=max(vpemax,valppem.max()) + + axc.axis([0, 2.0, 0.0, 1.1*vpemax]) #par defaut dans matplotlib sans seaborn + def singuess(t): + if (lasin(t)>=mu(t)): + return(1-(mu(t)/lasin(t))) + else: + return(0.0) + + psing=periodise(singuess) + axsinwicetwinter(axc) + axc.plot(et,[psing(s) for s in et],label="guess(t)",color='black',dashes=[6,2]) + + + #### proba emergence cas step + #axc.subplot(nblignes,2,5) + axc=ax[2,0] + mettrelettreax(axc,'E',gauche=0.0) + axabsticks(axc) + axc.set_yticks((0.0,0.2,0.4,0.6)) + axc.set_yticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}")) + if (offsetlam<1): + axbandegris(axc,1-gamma, 1.0) + axbandegris(axc,1+1-gamma,2.0) + + axc.set_ylabel(r"\bf{Emergence probability} ${p_e(t_0T,T)}$") + axc.axis([0, 2.0, 0.0, 1.1*vpemax]) #par defaut dans matplotlib sans seaborn + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + #plt.xlabel("Introduction time of the infected") + + for T in lT: + pem=nncalculepe(la,La,mu,Mu,T) + ppem=periodise(pem) + color=plt.cm.YlOrRd(cs(T)) + axc.plot(et,[ppem(s) for s in et],label=r"T="+str(T),color=color) + + def stepguess(t): + if (la(t)>=mu(t)): + return(1-(mu(t)/la(t))) + else: + return(0.0) + + psg=periodise(stepguess) + #plt.plot(et,[psg(s) for s in et],label=r"guess($t_0$)",color='black',dashes=[6,2]) + axc.plot(et,[psg(s) for s in et],color='black',dashes=[6,2]) + #flechewinter(0.6,bdp.trouvepremierzero(psg,0.1,0.9),1.0) + tstarstep=phi(1)/(lzero-1) + astep,bstep=trouvepremieretdernierzero(ppem,0.1,1.0) + print("tstarstep=",tstarstep) + axbandegris(axc,tstarstep,1-gamma,alpha=0.5) + axbandegris(axc,1+tstarstep,1+1-gamma,alpha=0.5) + #bandegris(astep,bstep,alpha=0.5) + #bandegris(1+astep,1+bstep,alpha=0.5) + + legend = axc.legend(fontsize=12,loc='lower left') + frame = legend.get_frame() + frame.set_linewidth(0.0) + + plt.savefig("ratesandemergenceproba.pdf",dpi=300,pad_inches=0.1) + + + +######################################################################### +#### vendredi premier février 2019 : attaquons nous a la seconde figure + +#tout d'abord on reprend notre calcul de la strategie de controle optimale + +def nsimusin(lzero=4.0,C=0.2,voir=False): + #ne pas ajuster le lzero : cela doit etre fait dans la fonction qui appelle celle ci + + def la(t): + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def La(t): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + pla=periodise(la) + pmu=periodise(mu) + pLa=periodiseetintegre(La) + pMu=periodiseetintegre(Mu) + + def sinqextmoy(tun,rhom,C=C): + #on peut prendre tdeux >=1, mais il faut imposer rhom >= C + #et prendre les fonctions periodisees + assert(rhom >= C) + deltat=C/rhom + tdeux=tun+deltat + if (tdeux <=1): + def rho(t): + return(rhom if ((tun$=""" +'{:.4f}'.format(m)+r""" for $t_1$=""" + str(xopt)+r""", $\rho_M=$"""+str(yopt) + ax.set_title(titre) + plt.show() +def nnsinvoir(z,C=0.2,tunmax=0.5): + N=(z.shape)[0] + tun=np.arange(0.0,tunmax,tunmax/N) + rho=(1.0/N)+np.arange(C,1.0,1.0/N) + + X,Y=np.meshgrid(rho,tun) + fig, ax = plt.subplots() + zmin=z.min() + zmax=z.max() + c = ax.pcolormesh(X,Y,z,cmap='RdBu',vmin=zmin,vmax=zmax) + ax.set_ylabel(r"$t_1$") + ax.set_xlabel(r"$\rho_M$") + m=z.min() + y=np.where(z==m) + #xopt=(y[0][0]/N)*tunmax + xopt=tun[y[0][0]] + #yopt=(y[1][0]/N)*(1-C)+C + yopt=rho[y[1][0]] + titre=r"Mean emergence probability: Cost="+str(C)+"\n"+ r""" minimal $$=""" +'{:.4f}'.format(m)+r""" for $t_1$=""" + str(xopt)+r""", $\rho_M=$"""+'{:.4f}'.format(yopt) + ax.set_title(titre) + fig.colorbar(c,ax=ax) + plt.show() + +#jeudi 7 fevrier : il faut ecrier une autre fonction peifmoy, qui tienne compte que l'on peut avoir tdeux plus grand que 1 + + +def npeinfmoy(la,La,mu,Mu,N=200,voir=False): + def phi(t): + return(La(t)-Mu(t)) #les fonctions en argument sont deja periodisees + + lesphi=np.array([phi(i/N) for i in np.arange(2*N+1)]) + def peinf(i): + t=i/N + #print(i,t,la(t),mu(t),La(t),Mu(t)) + if (la(t) <= mu(t)): + return(0.0) + z=lesphi[i+1:i+N]- lesphi[i] #on regarde la courbe de taux sur une periode apres i + if (z.min()<=0): #s 'il y a un piege, la proba d'extinctinest nulle + return(0.0) + else: + return(1-(mu(t)/la(t))) + ab=np.arange(N+1) + val=np.array([peinf(i) for i in ab]) + return(val.mean()) + + + +def bacanpeinfmoy(la,La,mu,Mu,N=200,voir=False,nbtau=10): + T=40 #plus grand cela fait exploser l'integration numerique + pe=nncalculepe(la,La,mu,Mu,T,limit=100) + ppe=periodise(pe) + def phi(t): + return(La(t)-Mu(t)) #les fonctions en argument sont deja periodisees + + lesphi=np.array([phi(i/N) for i in np.arange(2*N+1)]) + def peinf(i): + t=i/N + #print(i,t,la(t),mu(t),La(t),Mu(t)) + if (la(t) <= mu(t)): + return(0.0) + z=lesphi[i+1:i+N]- lesphi[i] #on regarde la courbe de taux sur une periode apres i + #print("z=",z) + if (z.min()<=0): #s 'il y a un piege, la proba d'extinctinest nulle + return(0.0) + else: + return(1-(mu(t)/la(t))) + ab=np.arange(N+1) + val=np.array([peinf(i) for i in ab]) + def pf(t): + "on fabrique la fonction a partir des valeurs entieres" + return(val[int(N*t)]) + ppf=periodise(pf) + + Lambda=[[la]] + Mu=[mu] + z=vfbaca(Lambda,Mu,T=T,tv=False,nbtau=nbtau) + z=z.reshape(len(z),) + def pz(t): + return(z[int(len(z)*t)]) + ppz=periodise(pz) + + + if voir: + abcisse=np.arange(0.0,2.0,1.0/(2*N)) + ordonnee=[ppf(t) for t in abcisse] + #print("ordonnee=",ordonnee) + plt.plot(abcisse,ordonnee,label=r"$p_{e,\infty}$") + plt.plot(abcisse,[ppe(t) for t in abcisse],label=r"$p_{e,T}(tT)$") + plt.plot(abcisse,[ppz(t) for t in abcisse],label=r"$ppz(t)$") + plt.legend() + #on fait une integration par la methode des trapezes + trapint=np.trapz(val,ab/N) + directint=quad(ppf,0,1)[0] + pemoy=(quad(pe,0,1,limit=100)[0]) + pzmoy=(quad(pz,0,1,limit=100)[0]) + #print("trapint=",trapint,"directint=",directint,"moyenne de val",val.mean()) + print("trapint=",trapint,"directint=",directint,"moyenne de val",val.mean(),"pemoy=",pemoy,"pzmoy=",pzmoy,"z.mean",z.mean()) + return(np.trapz(val,ab/N)) + + +#################################################################################### +########### vendredi 8 fevrier 2019 +########### enfin je m'attaque a la figure 2 + +def controltauxphietpe(gamma=0.3,lzero=3.0,tunmax=1.0,T=70,nbpt=100,C=0.2,beamer=False,tunstepopt=0.2,rhomstepopt=0.85): + r"on fixe a priori une grande periode" + #couleurcontrole='lightblue' + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + def phi(t): + return(La(t) -Mu(t)) + #il faut ajuster le lzerosin pour que les processus aient le meme rzero + + def lasin(t): + return(lzero*(1-gamma)*(1+np.sin(2*pi*t))) + def Lasin(t): + return(lzero*(1-gamma)*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + def phisin(t): + return(Lasin(t) -Mu(t)) + + def grise(): + bandegris(1-gamma, 1.0) + bandegris(1+1-gamma,2.0) + def bandegrisbarre(a,b): + plt.axvspan(a, b, color='lightgray', alpha=0.2, lw=0,hatch='///') + def grisebarre(): + bandegrisbarre(1-gamma, 1.0) + bandegrisbarre(1+1-gamma,2.0) + def flechewinter(hauteur,xmin,xmax,couleur='silver'): + hdelta=0.05 + plt.arrow(xmax-hdelta,hauteur,-(xmax-xmin)+2*hdelta,0,shape='full', length_includes_head=True,color=couleur,linewidth=10,head_starts_at_zero=True) + + + plt.rc('lines', linewidth=3.5, color='b') + + atracer,asin,bsin=trouvelsinmu(lzero*(1-gamma),0.0) + def singrise(): + if (atracer): + bandegris(asin,bsin) + bandegris(asin+1,bsin+1) + + #print("phi(1),phisin(1)",phi(1),phisin(1)) + #on recupere les valeurs de tun et rhom qui donnent le minimum + #tun,rhom=minopt(nbpt,C,lzero,tunmax=1.0) + lzeromodif=round(lzero*(1-gamma)*100)/100 + + tun,rhom=minopt(nbpt,C,lzeromodif,tunmax=1.0,prefixe="unif") + + #on definit les fonctions corespondant au controle optimal pour la sinusoide + deltat=C/rhom + tdeux=tun+deltat + print("optimum sinus:tun,tdeux,rhom=",tun,tdeux,rhom) + if (tdeux <=1): + def rhosinopt(t): + return(rhom if ((tun$=""" +'{:.4f}'.format(m)+r""" for $t_1$=""" + str(xopt)+r""", $\rho_M=$"""+'{:.4f}'.format(yopt)+"\n"+r"minimal $R_0$="+'{:.4f}'.format(rzeropt) +r""" for $t_1$=""" + str(tunrzeropt)+r""", $\rho_M=$"""+'{:.4f}'.format(rhomrzeropt) + #ax.set_title(titre) #a remettre lorsque l'on veut lire les valeurs mai 2019) + #fig.colorbar(c,ax=ax) #pas de barre de couleurs le 29 mai 2019 + #ax.plot(rhomrzeropt,tunrzeropt,marker='x',color=Lacolormarker,markersize=20) + ax.plot(rhomrzeropt,tunrzeropt,marker='x',color='black',markersize=30) + ax.plot(yopt,xopt,marker='x',color=couleurcontrole,markersize=30) + #cercle=plt.Circle((yopt,xopt),0.01,color=couleurcontrole,fill=True) + #ax.add_artist(cercle) + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['bottom'].set_visible(False) + ax.spines['left'].set_visible(False) + plt.text(rhomrzeropt-0.05,tunrzeropt+0.05,r"${R_{0,opt}}$",color='black',fontsize=40) + plt.savefig("optimalpeandrzero.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + plt.show() + + + +#on rajoute en parametre la fonction densite d'introduction pour le calcul de la proba d'emergence + +def demicercledensity(t): + assert((t>=0) and (t<=1)) + return((8/np.pi)*np.sqrt(t-t*t)) +def unifdensity(t): + return(1.0) +def betadensity(t): + return(beta.pdf(t,2.0,5.0)) +def sinusdensity(t): + return(1+np.sin(2*np.pi*t)) + +#pour voir la densité du demicercle +def voirdemicercle(): + t=np.linspace(0,1,100) + plt.plot(t,[demicercledensity(s) for s in t]) + print("surface=",(quad(demicercledensity,0,1))[0]) + +#pour voir la densité beta +def voirbeta(a=4,b=5): + t=np.linspace(0,1,100) + plt.plot(t,[beta.pdf(s,a,b) for s in t]) + + +def nnsimusin(intdens,lzero=4.0,C=0.2,voir=False): + #ne pas ajuster le lzero : cela doit etre fait dans la fonction qui appelle celle ci + + def la(t): + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + def La(t): + return(lzero*(t + (1/(2*pi))*(1-np.cos(2*pi*t)))) + pla=periodise(la) + pmu=periodise(mu) + pLa=periodiseetintegre(La) + pMu=periodiseetintegre(Mu) + + def sinqextmoy(tun,rhom,C=C): + #on peut prendre tdeux >=1, mais il faut imposer rhom >= C + #et prendre les fonctions periodisees + assert(rhom >= C) + deltat=C/rhom + tdeux=tun+deltat + if (tdeux <=1): + def rho(t): + return(rhom if ((tun1): + return(1-(1/a)) + else: + return(0.0) + + #nbpts=nb*tau + Lambda=[[lasin]] + VMu=[mu] + z=vfbaca(Lambda,VMu,T=T,tv=False) + z=z.reshape(len(z),) + pem=nncalculepe(lasin,Lasin,mu,Mu,T) + t=np.linspace(0.0,1.0,len(z)) + valpem=np.array([pem(s) for s in t]) + vg=np.array([guess(s) for s in t]) + plt.title("Linear BD,sinusoidal birth rate") + plt.plot(t,z,label="simulation de Bacaer") + plt.plot(t,valpem,label="calcul integral exact") + plt.plot(t,vg,label=r"guess : $(1-\lambda/\mu)^+$") + plt.legend() + print("difference max",np.abs(valpem-z).max()) + #on n'a pas zero car les tableaux n'ont pas meme dimension + return((z,valpem)) + +def produit(liste): + if (len(liste)==0): + return(1) + else: + return(reduce((lambda x,y : x*y),liste)) + +def guess(lam,mu): + r"lam et mu sont un tableau et un vecteur de taux : retourne les probabilites d extinction" + d=len(mu) + numerateur=produit([lam[i%d,(i+1)%d] for i in range(d)]) -produit(mu) + if (numerateur<=0.0): + return(np.array([0.0 for i in range(d)])) + denominateur=np.array([np.array([produit([mu[j%d] for j in range(i,i+k)])*produit([lam[i%d,(i+1)%d] for i in range(i+k,i+d)]) for k in range(d)]).sum() for i in range(d)]) + #print("lam,mu,numerateur,denominateur",lam,mu,numerateur,denominateur) + return(numerateur/denominateur) + +def rzerovar(lam,mu): + d=len(mu) + numerateur=produit([lam[i%d,(i+1)%d] for i in range(d)]) + return(numerateur/produit(mu)) + +def testeguess(): + lam=np.array([[0,1],[2,3]]) + mu=np.array([1,1]) + return(guess(lam,mu)) + +def lambrecht(t): + return(0.001044*t*(t-12.286)*np.sqrt(ppos(32.461-t))) + +def tlbrech(): + t=np.linspace(15,35,100) + plt.plot(t,[lambrecht(s) for s in t]) + plt.title("Transmission probability as a function of temperature :\n the model of Lambrecht et al.") + plt.xlabel("Temperature") + plt.savefig("Lambrecht.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + + +def fsamodel(Tmoy=25,DT=20): + muunc=(1/6.0)+1.0/(75*365) + av=0.25 + phihv=0.5 + nvsnh=5 #le rapport Nv/Nh + def l34(x): + return(1/9.0) + def l12(x): + return(1/6.0) + def l23(x): + return(av*phihv*nvsnh) + def mu1(x): + return(muunc) + def mu2(x): + return(muunc) + def mu3(x): + return((1.0/5) + (1.0/9.0)) + def mu4(x): + return(1.0/9) + def zero(x): + return(0.0) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def l41(t): + return(av*lambrecht(temperature(t))) + + Mu=[mu1,mu2,mu3,mu4] + Lambda=[[zero,l12,zero,zero],[zero,zero,l23,zero],[zero,zero,zero,l34],[l41,zero,zero,zero]] + return(Lambda,Mu) + + +def tfsa(DT=10,T=365,nbtau=5): + Lambda,Mu=fsamodel(DT=DT) + z=vfbaca(Lambda,Mu,T=T,tv=False,nbtau=nbtau) + t=np.linspace(0.0,1.0,len(z)) + vg=np.array([ppos(guess(afvalue2(Lambda,s),afvalue(Mu,s))) for s in t]) + plt.title("FSA process : emergence probability for one exposed human") + plt.plot(t,z[:,0],label=r"$p_e$") + #plt.plot(t,z,label=r"$p_e$") + #plt.plot(t,[Lambda[3][0](s) for s in t],label=r"$\lambda_{I_V,E_H}$") + plt.plot(t,vg[:,0],label=r"guess") + #plt.plot(t,vg,label=r"guess") + plt.xlabel("Introduction time") + plt.legend() + return(z) + + +############ Vendredi 1er Mars 2019 : modele SA amerique du sud de Zhang et al + +def tmosden(Tmoy=25,Topt=28,DT=10,damping=50): + + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def denfunc(t): + return(np.exp(-(t-Topt)**2/damping)) + def mosquitodensity(t): + return(denfunc(temperature(t))) + + f,ax=plt.subplots(nrows=3, ncols=1, figsize=[8, 12]) + t=np.linspace(0,1,100) + + plt.subplot(3,1,1) + plt.plot(t,[temperature(s) for s in t],label="temperature") + plt.title("Temperature") + plt.xlabel("relative time in the year") + plt.legend() + plt.subplot(3,1,2) + plt.plot(t,[10*mosquitodensity(s) for s in t],label="mosquitodensity") + plt.title("\n Mosquito density vs temperature : Zhang et al model") + plt.xlabel("relative time in the year") + plt.legend() + + plt.subplot(3,1,3) + tempabs=np.linspace(Tmoy-0.5*DT,Tmoy+0.5*DT,100) + vdf=np.array([denfunc(s) for s in tempabs]) + print("denfunc, max, min,ratio",vdf.max(),vdf.min(),vdf.max()/vdf.min()) + plt.plot(tempabs,[denfunc(s) for s in tempabs],label="mosquitodensity") + plt.xlabel(r"Temperature $T$") + plt.legend() + +def samodel(Tmoy=25,Topt=25,DT=20,muv=0.4,kc=10): + betatilde=0.3 + epsilonv=1.0/8 + epsilonh=1.0/7 + muh=0.2 + #muv=1.0 + damping=5.3*5.3 + def l34(x): + return(epsilonv) + def l12(x): + return(epsilonh) + def l23(x): + return(betatilde*kc*np.exp(-(temperature(x)-Topt)**2/damping)) + def mu1(x): + return(epsilonh) + def mu2(x): + return(muh) + def mu3(x): + return(epsilonv+muv) + def mu4(x): + return(muv) + def zero(x): + return(0.0) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def l41(t): + return(betatilde) + + Mu=[mu1,mu2,mu3,mu4] + Lambda=[[zero,l12,zero,zero],[zero,zero,l23,zero],[zero,zero,zero,l34],[l41,zero,zero,zero]] + return(Lambda,Mu) + + +def tsa(Tmoy=25,Topt=28,DT=20,T=365,nbtau=5,muv=0.5,kc=10): + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + Lambda,Mu=samodel(Tmoy=Tmoy,Topt=Topt,DT=DT,muv=muv,kc=kc) + z=vfbaca(Lambda,Mu,T=T,tv=False,nbtau=nbtau) + + #f,ax1=plt.subplots(figsize=[8, 12]) + f,ax1=plt.subplots() + + t=np.linspace(0.0,1.0,len(z)) + vg=np.array([ppos(guess(afvalue2(Lambda,s),afvalue(Mu,s))) for s in t]) + lud=Lambda[1][2] + lihev=np.array([lud(s) for s in t]) + + ax1.set_title("Emergence probability for one exposed human") + ax1.plot(t,z[:,0],label=r"$p_e$") + #plt.plot(t,[temperature(s) for s in t],label=r"temperature") + ax1.plot(t,vg[:,0],label=r"guess",color='b') + ax1.set_xlabel("Introduction time") + ax1.set_ylabel(r"$p_e$",color='b') + ax1.tick_params('y',color='b') + ax1.legend() + + x=np.array([rzerovar(afvalue2(Lambda,s),afvalue(Mu,s)) for s in t]) + + #k=kde.gaussian_kde(x) + # plt.plot(x,k(x),label=r"density estimation") + # plt.xlabel(r"$R_0$") + # plt.legend() + + ax2=ax1.twinx() + ax2.set_ylabel(r"time varying $R_0$",color='r') + ax2.plot(t,lihev,label=r"$\lambda_{1,2}$",color='green') + + ax2.plot(t,x,color='r',label=r"$R_0$") + ax1.tick_params('y',color='r') + ax2.legend() + plt.savefig("SA1.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + + plt.show() + return(z) + + +################ mrecredi 6 mars 2019 : emergence for a vector borne disease +########### essayons de mettre en place une strategie de controle +############ pour le cas ou la periodicite est uniquement dans le rapport nh/nv +########### on diminue ce facteur de rhom pendant l'intervalle de temps t1,t2 + +def tsacontrol(Tmoy=25,DT=8,T=365,nbtau=5,muv=0.5,tun=0.4,rhom=0.9,C=0.2,kc=10): + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + Lambda,Mu=samodel(Tmoy=Tmoy,DT=DT,muv=muv,kc=kc) + z=vfbaca(Lambda,Mu,T=T,tv=False,nbtau=nbtau) + + fC=modelsacontrolled(Tmoy=Tmoy,DT=DT,muv=muv,kc=kc,C=C) + LambdaC,MuC=fC(tun,rhom) + zC=vfbaca(LambdaC,MuC,T=T,tv=False,nbtau=nbtau) + + #f,ax1=plt.subplots(figsize=[8, 12]) + f,ax=plt.subplots(nrows=2,ncols=2, figsize=[11, 16]) + + t=np.linspace(0.0,1.0,len(z)) + vg=np.array([ppos(guess(afvalue2(Lambda,s),afvalue(Mu,s))) for s in t]) + vgC=np.array([ppos(guess(afvalue2(LambdaC,s),afvalue(MuC,s))) for s in t]) + lud=Lambda[1][2] + #print("lud",lud) + lihev=np.array([lud(s) for s in t]) + ludC=LambdaC[1][2] + lihevC=np.array([ludC(s) for s in t]) + + plt.subplot(2,2,1) + plt.title(r"Transmission rate $\lambda_{I^H,E^V}$") + plt.plot(t,lihev,color='black') + + + plt.subplot(2,2,3) + + plt.title("Emergence probability for one exposed human") + plt.plot(t,z[:,0],label=r"$p_e$",color='black') + #plt.plot(t,[temperature(s) for s in t],label=r"temperature") + plt.plot(t,vg[:,0],label=r"guess",color='black',dashes=[6,2]) + plt.xlabel("Introduction time") + plt.legend() + + plt.subplot(2,2,2) + plt.title(r"Controlled Transmission rate $\lambda_{I^H,E^V}$") + plt.plot(t,lihevC,color='black') + + + plt.subplot(2,2,4) + + plt.title("Controlled Emergence probability for one exposed human") + plt.plot(t,zC[:,0],label=r"$p_e$",color='black') + #plt.plot(t,[temperature(s) for s in t],label=r"temperature") + plt.plot(t,vgC[:,0],label=r"guess",color='black',dashes=[6,2]) + plt.xlabel("Introduction time") + plt.legend() + + + plt.savefig("SA1control.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + + + + + + +def modelsacontrolled(Tmoy=25,DT=20,muv=0.4,kc=10,C=0.2): + r""" retourne une fonction qui prend pour parametre tun et rhom et retourne les parametres controles, c'est a dire les fonctions lambbda et mu avec la diminution du facteur kc entre tun et tdeux""" + betatilde=0.3 + epsilonv=1.0/8 + epsilonh=1.0/7 + muh=0.2 + #muv=1.0 + damping=5.3*5.3 + def l34(x): + return(epsilonv) + def l12(x): + return(epsilonh) + def mu1(x): + return(epsilonh) + def mu2(x): + return(muh) + def mu3(x): + return(epsilonv+muv) + def mu4(x): + return(muv) + def zero(x): + return(0.0) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def l41(t): + return(betatilde) + + def sinqextmoy(tun,rhom): + #on peut prendre tdeux >=1, mais il faut imposer rhom >= C + #et prendre les fonctions periodisees + assert(rhom >= C) + deltat=C/rhom + tdeux=tun+deltat + #print("tun,tdeux,rhom",tun,tdeux,rhom) + if (tdeux <=1): + def facteur(x): + if ((x<=tun) or (x>=tdeux)): + return(1.0) + else: + return(1-rhom) + else: + def facteur(x): + if ((x<= tdeux -1) or (x>=tun)): + return(1-rhom) + else: + return(1.0) + #print("facteur x=milieu tun,tdeux",facteur(0.5*(tun+tdeux))) + def l23(x): + return(facteur(x)*betatilde*kc*np.exp(-(temperature(x)-Tmoy)**2/damping)) + + + Mu=[mu1,mu2,mu3,mu4] + Lambda=[[zero,l12,zero,zero],[zero,zero,l23,zero],[zero,zero,zero,l34],[l41,zero,zero,zero]] + return(Lambda,Mu) + + return(sinqextmoy) + + + +def samodelbis(Tmoy=25,DT=20,muv=0.4,kc=10,betatilde=0.3): + r""" on remplace le betatilde, proba de transmission fixe, par une fonction de Lambrecht""" + #betatilde=0.3 + epsilonv=1.0/8 + epsilonh=1.0/7 + muh=0.2 + #muv=1.0 + damping=5.3*5.3 + def l34(x): + return(epsilonv) + def l12(x): + return(epsilonh) + def l23(x): + return(2*betatilde*lambrecht(temperature(x))*kc*np.exp(-(temperature(x)-Tmoy)**2/damping)) + def mu1(x): + return(epsilonh) + def mu2(x): + return(muh) + def mu3(x): + return(epsilonv+muv) + def mu4(x): + return(muv) + def zero(x): + return(0.0) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def l41(t): + #modif du 22 mai + #return(2*betatilde*lambrecht(temperature(t))) + return(lambrecht(temperature(t))) + + Mu=[mu1,mu2,mu3,mu4] + Lambda=[[zero,l12,zero,zero],[zero,zero,l23,zero],[zero,zero,zero,l34],[l41,zero,zero,zero]] + return(Lambda,Mu) + + + +def modelsacontrolledbis(Tmoy=25,DT=20,muv=0.4,kc=10,C=0.2,betatilde=0.3): + r""" retourne une fonction qui prend pour parametre tun et rhom et retourne les parametres controles, c'est a dire les fonctions lambbda et mu avec la diminution du facteur kc entre tun et tdeux : pour le modele bis""" + #betatilde=0.3 + epsilonv=1.0/8 + epsilonh=1.0/7 + muh=0.2 + #muv=1.0 + damping=5.3*5.3 + def l34(x): + return(epsilonv) + def l12(x): + return(epsilonh) + def mu1(x): + return(epsilonh) + def mu2(x): + return(muh) + def mu3(x): + return(epsilonv+muv) + def mu4(x): + return(muv) + def zero(x): + return(0.0) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + def l41(t): + return(2*betatilde*lambrecht(temperature(t))) + + def sinqextmoy(tun,rhom): + #on peut prendre tdeux >=1, mais il faut imposer rhom >= C + #et prendre les fonctions periodisees + assert(rhom >= C) + deltat=C/rhom + tdeux=tun+deltat + #print("tun,tdeux,rhom",tun,tdeux,rhom) + if (tdeux <=1): + def facteur(x): + if ((x<=tun) or (x>=tdeux)): + return(1.0) + else: + return(1-rhom) + else: + def facteur(x): + if ((x<= tdeux -1) or (x>=tun)): + return(1-rhom) + else: + return(1.0) + #print("facteur x=milieu tun,tdeux",facteur(0.5*(tun+tdeux))) + def l23(x): + return(facteur(x)*2*betatilde*lambrecht(temperature(x))*kc*np.exp(-(temperature(x)-Tmoy)**2/damping)) + + + Mu=[mu1,mu2,mu3,mu4] + Lambda=[[zero,l12,zero,zero],[zero,zero,l23,zero],[zero,zero,zero,l34],[l41,zero,zero,zero]] + return(Lambda,Mu) + + return(sinqextmoy) + + +def tsacontrolbis(Tmoy=25,DT=8,T=365,nbtau=5,muv=0.5,tun=0.4,rhom=0.9,C=0.2,kc=10,betatilde=0.3): + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + Lambda,Mu=samodelbis(Tmoy=Tmoy,DT=DT,muv=muv,kc=kc,betatilde=betatilde) + z=vfbaca(Lambda,Mu,T=T,tv=False,nbtau=nbtau) + + fC=modelsacontrolledbis(Tmoy=Tmoy,DT=DT,muv=muv,kc=kc,C=C,betatilde=betatilde) + LambdaC,MuC=fC(tun,rhom) + zC=vfbaca(LambdaC,MuC,T=T,tv=False,nbtau=nbtau) + + + t=np.linspace(0.0,1.0,len(z)) + def guessori(s): + return(ppos(guess(afvalue2(Lambda,s),afvalue(Mu,s)))) + def guesscontrol(s): + return(ppos(guess(afvalue2(LambdaC,s),afvalue(MuC,s)))) + + vg=np.array([guessori(s) for s in t]) + vgC=np.array([guesscontrol(s) for s in t]) + lud=Lambda[1][2] + #print("lud",lud) + lihev=np.array([lud(s) for s in t]) + ludC=LambdaC[1][2] + lihevC=np.array([ludC(s) for s in t]) + + + + lquatreun=Lambda[3][0] + liveh=np.array([lquatreun(s) for s in t]) + lquatreunC=LambdaC[3][0] + livehC=np.array([lquatreunC(s) for s in t]) + + #f,ax1=plt.subplots(figsize=[8, 12]) + f,(ax1,ax2,ax3)=plt.subplots(nrows=3,ncols=1, figsize=[16,20]) + f.subplots_adjust(top=0.9,bottom=0.1,left=0.125,right=0.9,hspace=0.5) + + mettrelettreax(ax1,'A') + mettrelettreax(ax2,'B') + mettrelettreax(ax3,'C') + + #plt.subplot(2,1,2) + for a in (ax1,ax2,ax3): + axdemiabsticks(a) + + ax3.axis([0.0,1.0,0.0,0.4]) + ax2.axis([0.0,1.0,0.0,4.0]) + ax1.axis([0.0,1.0,0.6,1.0]) + + ax3.set_ylabel(r"\bf{Emergence probability}",fontsize=22) + ax3.plot(t,z[:,0],label=r"$p_e$",color='black') + #plt.plot(t,vg[:,0],label=r"guess",color='black',dashes=[6,2]) #le 4 juin + ax3.set_xlabel(r"\bf{Time of Pathogen introduction} ${t_0}$",fontsize=22) + + ax3.plot(t,zC[:,0],label=r"$p_e$",color=couleurcontrole,linestyle=':') + #plt.plot(t,vgC[:,0],label=r"guess",dashes=[6,2],color=couleurcontrole) + + #a1,b1=trouvepremieretdernierzero(guesscontrol,0.0,1.0,tolerance=0.01) + def pemcontrol(s): + N=len(zC[:,0])-1 + i=int(s*N) + return(zC[i,0]) + a1,b1=trouvepremieretdernierzero(pemcontrol,0.0,1.0,tolerance=0.03) + tdeux=tun+(C/rhom) + axbandegris(ax3,a1,tun,couleur=couleurcontrole,alpha=0.15) + axbandegris(ax3,tun,tdeux,couleur=couleurcontrole,alpha=0.3) + + #plt.subplot(2,1,1) + ax2.set_ylabel(r"$\lambda_{I^H,E^V}$",fontsize=30) + #ax2.set_ylabel(r"\bf{Transmission rate} $\lambda_{I^H,E^V}$",fontsize=22) + ax2.plot(t,lihev,color='black') + ax2.plot(t,lihevC,color=couleurcontrole,linestyle=':') + + #ax1.set_ylabel(r"\bf{Transmission rate} $\lambda_{I^V,E^H}$",fontsize=22) + ax1.set_ylabel(r" $\lambda_{I^V,E^H}$",fontsize=30) + ax1.plot(t,liveh,color='black') + + + plt.savefig("SA2control.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + + + +####################################"" Jeudi 14 mars +################## trace de la heatmap de pe dans le cas step +#################" pour voir qu'il ya a toute une variete de strategies optimales + + + +#### vendredi premier février 2019 : attaquons nous a la seconde figure + +#tout d'abord on reprend notre calcul de la strategie de controle optimale + +def nsimustep(lzero=4.0,C=0.2,voir=False,gamma=0.3): + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + + return(lzero*(1+np.sin(2*pi*t))) + def mu(t): + return(1.0) + def Mu(t): + return(t) + pla=periodise(la) + pmu=periodise(mu) + pLa=periodiseetintegre(La) + pMu=periodiseetintegre(Mu) + + def sinqextmoy(tun,rhom,C=C): + #on peut prendre tdeux >=1, mais il faut imposer rhom >= C + #et prendre les fonctions periodisees + assert(rhom >= C) + deltat=C/rhom + tdeux=tun+deltat + if (tdeux <=1): + def rho(t): + return(rhom if ((tun$", transform=ax2.transAxes, + fontsize=40) + #ax3.set_title(r"$$",fontsize=30) + #ax3.axis([0.0,0.0,0.0,0.0]) + #enlevecadreax(ax3) + + #ax2.set_ylabel(r"\bf{Start of control} ${t_1}$",fontsize=40) + ax2.set_xlabel(r"\bf{Intensity of control} ${\rho_M}$",fontsize=40,labelpad=10) + ax2.set_yticks((0.0,0.2,0.4,0.6,0.8)) + ax2.set_yticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}",r"\bf{0.8}")) + ax2.set_xticks((0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)) + ax2.set_xticklabels((r"\bf{0.2}",r"\bf{0.3}",r"\bf{0.4}",r"\bf{0.5}",r"\bf{0.6}",r"\bf{0.7}",r"\bf{0.8}",r"\bf{0.9}",r"\bf{1.0}")) + ax2.axis([C,1.1,0,1.0]) + m=zsin.min() + y=np.where(zsin==m) + xopt=tun[y[0][0]] + yopt=rho[y[1][0]] + ax2.plot(rhomrzeropt,tunrzeropt,marker='x',color='red',markersize=30,mew=4) + ax2.plot(yopt,xopt,marker='x',color=couleurcontrole,markersize=30,mew=4) + + + plt.savefig("heatmapstepandsinus.pdf",dpi=300,bbox_inches='tight',pad_inches=0.2) + plt.show() + + + ################# Jeudi 14 mars : maintenant on montre l'influence de la temperature locale et surtout quand elle est diferenete de la temperature otimale de developpement des vecteurs + +def tsageo(Tmoy1=25,Tmoy2=28,Topt=28,DT=8,T=365,nbtau=5,muv=0.4,tun=0.4,rhom=0.9,C=0.2,kc=10,complet=False,beamer=True): + def flechewinter(hauteur,xmin,xmax,couleur=cgris): + hdelta=0.05 + #plt.arrow(xmax-hdelta,hauteur,-(xmax-xmin)+2*hdelta,0,shape='full', length_includes_head=True,color=couleur,alpha=1.0,linewidth=15,head_starts_at_zero=True,head_width=5e-2) + plt.arrow(xmax,hauteur,-(xmax-xmin)+2*hdelta,0,shape='full', length_includes_head=True,color=couleur,alpha=1.0,linewidth=40,head_starts_at_zero=True,head_width=2e-2) + def temperature(t): + return(Tmoy+0.5*DT*np.sin(2*pi*t)) + Lambda1,Mu1=samodel(Tmoy=Tmoy1,Topt=Topt,DT=DT,muv=muv,kc=kc) + z1=vfbaca(Lambda1,Mu1,T=T,tv=False,nbtau=nbtau) + + Lambda2,Mu2=samodel(Tmoy=Tmoy2,Topt=Topt,DT=DT,muv=muv,kc=kc) + z2=vfbaca(Lambda2,Mu2,T=T,tv=False,nbtau=nbtau) + + #f,ax1=plt.subplots(figsize=[8, 12]) + nblignes= 2 + if (beamer): + f,ax=plt.subplots(nrows=nblignes,ncols=2, figsize=[16, 10]) + else: + f,ax=plt.subplots(nrows=nblignes,ncols=2, figsize=[10, 16]) + + #plt.tight_layout() + plt.subplots_adjust(left=0.065,right=0.97,bottom=0.1,top=0.9,hspace=0.5) + plt.rc('text', usetex=True) + + + t=np.linspace(0.0,1.0,len(z1)) + vg1=np.array([ppos(guess(afvalue2(Lambda1,s),afvalue(Mu1,s))) for s in t]) + vg2=np.array([ppos(guess(afvalue2(Lambda2,s),afvalue(Mu2,s))) for s in t]) + lud1=Lambda1[1][2] + #print("lud",lud) + lihev1=np.array([lud1(s) for s in t]) + lud2=Lambda2[1][2] + lihev2=np.array([lud2(s) for s in t]) + vmin=min(lihev1.min(),lihev2.min()) + vmax=max(lihev1.max(),lihev2.max()) + def fguess2(x): + if (x<1): + #print("int(x*len(z1))",int(x*len(z1))) + return(vg2[int(x*len(z1)),0]) + else: + return(vg2[-1,0]) + awinter,bwinter=trouvepremieretdernierzero(fguess2,0.5,1.0) + + #plt.subplot(nblignes,2,1) + axc=ax[0,0] + axdemiabsticks(axc) + mettrelettreax(axc,'A',taille=50,gauche=0) + axc.axis([0, 1.0, 0.0,4.0]) + axc.set_yticks((0.0, 1.0,2.0,3.0,4.0)) + axc.set_yticklabels((r'\bf{0}', r'\bf{1}', r'\bf{2}',r'\bf{3}',r'\bf{4}'), color='k', size=20) + axc.set_ylabel( r"\bf{Transmission rate} $\lambda_{I^V,E^H}$",fontsize=22) + axc.plot(t,lihev1,color='black') + + #plt.subplot(nblignes,2,2) + axc=ax[0,1] + axdemiabsticks(axc) + mettrelettreax(axc,'B',taille=50) + axc.axis([0, 1.0, 0.0,4.0]) + axc.set_yticks((0.0, 1.0,2.0,3.0,4.0)) + axc.set_yticklabels((r'\bf{0}', r'\bf{1}', r'\bf{2}',r'\bf{3}',r'\bf{4}'), color='k', size=20) + axc.plot(t,lihev2,color='black') + axbandegris(axc,awinter,bwinter) + + #plt.subplot(nblignes,2,3) + axc=ax[1,0] + axdemiabsticks(axc) + mettrelettreax(axc,'C',taille=50,gauche=0) + axc.axis([0, 1.0, 0.0,0.4]) + axc.set_yticks((0.0, 0.1,0.2,0.3)) + axc.set_yticklabels((r'\bf{0}', r'\bf{0.1}', r'\bf{0.2}',r'\bf{0.3}'), color='k', size=20) + if (not(complet)): + axc.set_ylabel(r"\bf{Emergence probability}",fontsize=22) + axc.plot(t,z1[:,0],label=r"$p_e$",color='black') + else: + axc.set_ylabel(r"\bf{Emergence probabilities} ${p_{e,i}}$",fontsize=22) + axc.plot(t,z1[:,0],label=r"Exposed Human",color='black') + axc.plot(t,z1[:,1],label=r"Infectious Human",color='black',dashes=[2,2]) + axc.plot(t,z1[:,2],label=r"Exposed Vector",color='red') + axc.plot(t,z1[:,3],label=r"Infectious Vector",color='red',dashes=[2,2]) + axc.plot(t,vg1[:,0],label=r"guess",color='black',dashes=[1,1],alpha=0.5) + + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + + #axc.subplot(nblignes,2,4) + axc=ax[1,1] + axdemiabsticks(axc) + mettrelettreax(axc,'D',taille=50) + axc.axis([0, 1.0, 0.0,0.4]) + axc.set_yticks((0.0, 0.1,0.2,0.3)) + axc.set_yticklabels((r'\bf{0}', r'\bf{0.1}', r'\bf{0.2}',r'\bf{0.3}'), color='k', size=20) + + if complet: + axc.set_yticks((0.0, 0.1,0.2,0.3)) + axc.set_yticklabels((r'\bf{0}', r'\bf{0.1}',r'\bf{0.2}',r'\bf{0.3}'), color='k', size=20) + axc.plot(t,z2[:,0],label=r"Exposed Human",color='black') + axc.plot(t,z2[:,1],label=r"Infectious Human",color='black',dashes=[2,2]) + axc.plot(t,z2[:,2],label=r"Exposed Vector",color='red') + axc.plot(t,z2[:,3],label=r"Infectious Vector",color='red',dashes=[2,2]) + else: + axc.plot(t,z2[:,0],label=r"Exposed Human",color='black') + axc.plot(t,vg2[:,0],label=r"guess",color='black',dashes=[1,1],alpha=0.5) + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + def fpe2(x): + return(z2[int(x*len(z1)),0]) + trueawinter=trouvepremierzero(fpe2,0.0,awinter,tolerance=4e-3) + maxguess2=(vg2[:,0]).max() + print("trueawinter=",trueawinter,"maxguess2",maxguess2) + a1,b1=trouvepremieretdernierzero(fpe2,0.02,0.98,tolerance=5e-3) + axbandegris(axc,awinter,bwinter) + axbandegris(axc,a1,awinter,couleur=cgris,alpha=0.5) + + nomfig="zikaemergence.pdf" + if complet: + nomfig="complet"+nomfig + plt.savefig(nomfig,dpi=300,bbox_inches='tight',pad_inches=0.2) + +Lacolormap='Greys' +Lacolormarker='cyan' + + +################################ Mardi 19 mars 2019 +def stepnaivevsoptimal(gamma=0.3,lzero=3.0/0.7,tunmax=1.0,T=70,nbpt=100,C=0.2,beamer=False): + r"on fixe a priori une grande periode" + #couleurcontrole='lightblue' + def mu(t): + return(1.0) + def Mu(t): + return(t) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + def phi(t): + return(La(t) -Mu(t)) + + pla=periodise(la) + ############################## on definit la fonction optimale pour le cas step + + tdso=1-gamma-C*lzero/(lzero-1) + rhoms=C/(1-gamma -tdso) + print("tdso=",tdso) + def lasopt(t): + if (t <=tdso): + return(lzero) + elif (t<= 1-gamma): + return(1) + else: + return(0.0) + plasopt=periodise(lasopt) + + ###################" puis on definit la fonction naive + rhomnaive=C/(1-gamma) + def lasnaive(t): + if (t<=1-gamma): + return(lzero*(1-rhomnaive)) + else: + return(0.0) + + plasnaive=periodise(lasnaive) + ############# la on va tracer les courbes + et=np.linspace(0,2,200) + t=np.linspace(0,1,100) + #plt.rcParams['figure.figsize'] = [16, 10] + if (beamer): + fig, ax = plt.subplots(figsize=[16, 10]) + else: + fig, ax = plt.subplots(figsize=[10, 16]) + + + ax.set_xlabel('time') + ax.set_ylabel('rates') + ax.plot(et,[plasnaive(s) for s in et],label=r"$\lambda_{naive}$",color='lightgray') + ax.plot(et,[plasopt(s) for s in et],label=r"$\lambda_{opt}$",color=couleurcontrole) + ax.plot(et,[pla(s) for s in et],label=r"$\lambda$",color='black') + rectopt=Rectangle((tdso,lzero*(1-rhoms)),width=1-gamma-tdso,height= lzero*rhoms,color=couleurcontrole,alpha=0.5,fill=True) + ax.add_patch(rectopt) + rectnaive=Rectangle((0.0,lzero*(1-rhomnaive)),width=1-gamma,height= lzero*rhomnaive,color='lightgrey',alpha=0.5,fill=True) + ax.add_patch(rectnaive) + + ax.legend() + ax.plot(et,[1.0 for s in et],label=r"$\mu(t)$",color='black',dashes=[6,2]) + #plt.axis('scaled') + ax.text(0.05,1.1,r"$\mu(t)$",fontsize=22) + ax.text(0.05,la(0.0)+0.1,r"$\lambda(t)$",fontsize=22,color='black') + + + plt.savefig("stepnaivevsoptimal.pdf",dpi=300) + + + +def trouvecontour(z,tun,rho,disc=0.0001): + w=np.where(z<=z.min()+disc) + lesx=w[0] + lesy=w[1] + xc=lesx[0] + yc=lesy[0] + lec=[(rho[yc],tun[xc])] + for i in range(1,len(lesx)): + if (lesx[i]> xc): + xc=lesx[i] + yc=lesy[i] + lec.append((rho[yc],tun[xc])) + + return(lec) + + +########################## Mercredi 20 mars 2020 : c'est le printemps +def explaincontrolstrategy(gamma=0.3,lzero=3.0,beamer=False,tun=0.0,tdeux=0.8,C=0.2): + r"Explication de la strategie de controle avec la courbe rho et la courbe lambda_rho" + #couleurcontrole='lightblue' + def mu(t): + return(1.0) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def rho(t): + if ((t>= tun) and (t<= tdeux)): + return(rhom) + else: + return(0.0) + def larho(t): + return(la(t)*(1-rho(t))) + + #on definit les fonctions corespondant au controle optimal pour la sinusoide + rhom=C/(tdeux -tun) + assert(rhom <=1) + + pla=periodise(la) + pmu=periodise(mu) + plarho=periodise(larho) + prho=periodise(rho) + + + ############# la on va tracer les courbes + et=np.linspace(0,2,200) + t=np.linspace(0,1,100) + #plt.rcParams['figure.figsize'] = [16, 10] + if (beamer): + f,ax=plt.subplots(nrows=2, ncols=1, figsize=[16, 11]) + else: + f,ax=plt.subplots(nrows=2, ncols=1, figsize=[11, 16]) + + + + plt.subplot(2,1,1) + plt.ylabel(r"$\rho(t)$",fontsize=22) + plt.plot(et,[prho(s) for s in et],color='black') + plt.axis([0.0, 2.0,0.0, 1.2*rhom]) + h=0.05 + plt.text(h/2,rhom+h/2,r"$\rho_M$",fontsize=22) + plt.text(tun-2*h,h/2,r"$t_1$",fontsize=22) + plt.text(tdeux+h,h/2,r"$t_2$",fontsize=22) + + + plt.subplot(2,1,2) + plt.ylabel(r"$\lambda_\rho(t)$",fontsize=22) + plt.axis([0.0, 2.0,0.0, 1.2*lzero]) + + plt.plot(et,[pla(s) for s in et],color='black',label=r"$\lambda(t)$") + plt.plot(et,[plarho(s) for s in et],color=couleurcontrole,label=r"$\lambda_\rho(t)$") + plt.legend(fontsize=16) + + + plt.savefig("explainstepcontrolstrategy.pdf",dpi=300) + +################################################################## +########## Mercredi 19 avril 2019 : winter coming pour le cas step + +def winteriscoming(gamma=0.3,lzero=3.0,lT=[0.2,1,5,20,50,100],lzmin=0.5,pgr=True,beamer=True): + #cgris='lightgrey' + def mu(t): + return(1.0) + def Mu(t): + return(t) + def laoff(t): + if (t< 1-gamma): + return(lzero) + else: + return(lzmin) + def Laoff(t): + if (t< 1-gamma): + return((lzero)*t) + else: + return(lzero*(1-gamma)+lzmin*(t-(1-gamma))) + def la(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def La(t): + if (t< 1-gamma): + return((lzero)*t) + else: + return(lzero*(1-gamma)) + def phi(t): + return(La(t) -Mu(t)) + + def phioff(t): + return(Laoff(t) -Mu(t)) + + def grise(): + bandegris(1-gamma, 1.0) + bandegris(1+1-gamma,2.0) + def axgrise(ax): + axbandegris(ax,1-gamma, 1.0) + axbandegris(ax,1+1-gamma, 1+1.0) + + def bandegrisbarre(a,b): + plt.axvspan(a, b, color='lightgray', alpha=0.5, lw=0,hatch='///') + def grisebarre(): + bandegrisbarre(1-gamma, 1.0) + bandegrisbarre(1+1-gamma,2.0) + + plt.rc('lines', linewidth=3.5, color='b') + + pla=periodise(la) + pmu=periodise(mu) + plaoff=periodise(laoff) + + et=np.linspace(0,2,200) + t=np.linspace(0,1,100) + #plt.rcParams['figure.figsize'] = [16, 10] + nblignes=3 if pgr else 2 + if (beamer): + f,ax=plt.subplots(nrows=nblignes, ncols=2, figsize=[16,13]) #pour les graphiques beamer + else: + f,ax=plt.subplots(nrows=nblignes, ncols=2, figsize=[10,16]) #pour les graphiques beamer + #f.tight_layout() + f.subplots_adjust(top=0.9,bottom=0.1,left=0.125,right=0.9,hspace=0.5) + + #plt.subplot(nblignes,2,1) + axc=ax[0,0] + axabsticks(axc) + axgrise(axc) + mettrelettreax(axc,'A',taille=50,gauche=0.0) + axc.set_ylabel(r"\bf{Birth and Death rates}",fontsize=20) + vpla=np.array([pla(s) for s in et]) + vplaoff=np.array([plaoff(s) for s in et]) + vmax=max(vpla.max(),vplaoff.max()) + axc.plot(et,vplaoff,label=r"$\lambda(t)$",color='black') + axc.plot(et,[pmu(s) for s in et],label=r"$\mu(t)$",color='black',dashes=[1,1]) + axc.axis([0.0, 2.0,0.0, 3.5]) #par defaut dans matplotlib sans seaborn + axc.text(0.05,1.1,r"$\mu(t)$",fontsize=22) + axc.text(0.05,laoff(0.0)+0.1,r"$\lambda(t)$",fontsize=22) + axc.set_yticks((0,1,2,3)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}",r"\bf{3}")) + + #plt.subplot(nblignes,2,2) + axc=ax[0,1] + mettrelettreax(axc,'B',taille=50) + axabsticks(axc) + axgrise(axc) + axc.plot(et,vpla,label=r"$\lambda(t)$",color='black') + axc.plot(et,[pmu(s) for s in et],label=r"$\mu(t)$",color='black',dashes=[1,1]) + axc.axis([0, 2.0, 0.0, 3.5]) #par defaut dans matplotlib sans seaborn + axc.set_yticks((0,1,2,3)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}",r"\bf{3}")) + + if (pgr): + ########################################## + ##### courbe phi dans le cas step + ################################## + #plt.subplot(nblignes,2,4) + axc=ax[1,1] + mettrelettreax(axc,'D',taille=50) + axabsticks(axc) + axgrise(axc) + axc.set_yticks((0,1,2)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}")) + pphi=periodiseetintegre(phi) + pphioff=periodiseetintegre(phioff) + vpphi=np.array([pphi(s) for s in et]) + vpphioff=np.array([pphioff(s) for s in et]) + vphimax=max(vpphi.max(),vpphioff.max()) + axc.axis([0, 2.0, 0.0, 1.1*vphimax]) #par defaut dans matplotlib sans seaborn + axc.plot(et,vpphi,label=r"$\varphi(t)$",color='black') + tstar=(phi(1))/(lzero-1) + minloc=phi(1) + #la ligne pointillee horizontale depuis le minimum local + abst=np.linspace(tstar,1,100) + #je la vire le 28 mai 2019 + #line=plt.plot(abst,np.full_like(a=abst,fill_value=minloc),dashes=[6,2],color='black')[0] + + h=0.1 + + def bgclair(ax): + #enfin la bande gris clair au lieu de la fleche winter + axbandegris(ax,tstar,1-gamma,alpha=0.5) + axbandegris(ax,1+tstar,1+1-gamma,alpha=0.5) + bgclair(axc) + hep=0.04 + axc.arrow(1.0-0.02, 0.72, -(1.0 -tstar)+0.05,0.0, color='black', head_length = 0.06,linewidth=7 ,head_width = 0.09, length_includes_head = True,shape='full') + axc.arrow(1.0+1.0-0.02, 1.48, -(1.0 -tstar)+0.05,0.0, color='black', head_length = 0.06,linewidth=7 ,head_width = 0.09, length_includes_head = True,shape='full') + print("tstar=",tstar) + + + ################################################################ + # on trace la courbe phi dans le cas step avec offset + #plt.subplot(nblignes,2,3) + axc=ax[1,0] + mettrelettreax(axc,'C',taille=50,gauche=0.0) + + axabsticks(axc) + axgrise(axc) + axc.set_yticks((0,1,2)) + axc.set_yticklabels((r"\bf{0}",r"\bf{1}",r"\bf{2}")) + + axc.set_ylabel(r"\bf{Integrated growth rate}",fontsize=18) + axc.text(0.2,0.7,r"$\varphi(t)$") + axc.axis([0, 2.0, 0.0, 1.1*vphimax]) #par defaut dans matplotlib sans seaborn + #plt.title(r"\bf{Integrated growth rate}",fontsize=22) + axc.plot(et,vpphioff,label=r"$\varphi(t)$",color='black') + + + #### proba emergence cas step avec offset + if (pgr): + #plt.subplot(nblignes,2,5) + axc=ax[2,0] + else: + #plt.subplot(nblignes,2,3) + axc=ax[1,0] + mettrelettreax(axc,'E',taille=50,gauche=0.0) + axgrise(axc) + axabsticks(axc) + axc.set_yticks((0,0.2,0.4,0.6,0.8)) + axc.set_yticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}",r"\bf{0.8}")) + axc.axis([0, 2.0, 0.0, 0.82]) #par defaut dans matplotlib sans seaborn + + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + axc.set_ylabel(r"\bf{Emergence probability}",fontsize=20) + axc.text(0.1,0.63,r"$p_e(t_0T,T)$",color='k') + vpemax=0.0 + cs=CubicSpline(np.array([0.0,10,100]),np.array([0.2,0.5,0.9])) + + for T in lT: + pem=nncalculepe(laoff,Laoff,mu,Mu,T) + ppem=periodise(pem) + valppem=np.array([ppem(s) for s in et]) + axc.plot(et,valppem,label=r"T="+str(T),color=plt.cm.YlOrRd(cs(T))) + vpemax=max(vpemax,valppem.max()) + + def offguess(t): + if (laoff(t)>=mu(t)): + return(1-(mu(t)/laoff(t))) + else: + return(0.0) + + poffg=periodise(offguess) + axc.plot(et,[poffg(s) for s in et],color='black',dashes=[6,2]) + legend = axc.legend(fontsize=12,loc='lower left') + frame = legend.get_frame() + frame.set_linewidth(0.0) + + #### proba emergence cas step + if (pgr): + #plt.subplot(nblignes,2,6) + axc=ax[2,1] + else: + #plt.subplot(nblignes,2,4) + axc=ax[1,1] + mettrelettreax(axc,'F',taille=50) + axabsticks(axc) + axgrise(axc) + axc.set_yticks((0,0.2,0.4,0.6,0.8)) + axc.set_yticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}",r"\bf{0.8}")) + + axc.axis([0, 2.0, 0.0, 0.82]) #par defaut dans matplotlib sans seaborn + axc.set_xlabel(r"\bf{Time of pathogen introduction} ${t_0}$",fontsize=22) + + for T in lT: + pem=nncalculepe(la,La,mu,Mu,T) + ppem=periodise(pem) + color=plt.cm.YlOrRd(cs(T)) + axc.plot(et,[ppem(s) for s in et],label=r"T="+str(T),color=color) + + bgclair(axc) + def stepguess(t): + if (la(t)>=mu(t)): + return(1-(mu(t)/la(t))) + else: + return(0.0) + + psg=periodise(stepguess) + axc.plot(et,[psg(s) for s in et],label=r"guess($t_0$)",color='black',dashes=[6,2]) + + if (pgr): + plt.savefig("winteriscoming.pdf",dpi=300,bbox_inches='tight') + else: + plt.savefig("winteriscomingsansphi.pdf",dpi=300,bbox_inches='tight') + + +def absticks(): + plt.xticks((0.0,0.5,1.0,1.5,2.0),(r"\bf{0}",r"\bf{0.5}",r"\bf{1.0}",r"\bf{1.5}",r"\bf{2.0}")) +def axabsticks(ax): + ax.set_xticks((0.0,0.5,1.0,1.5,2.0)) + ax.set_xticklabels((r"\bf{0}",r"\bf{0.5}",r"\bf{1.0}",r"\bf{1.5}",r"\bf{2.0}")) + +def demiabsticks(): + plt.xticks((0.0,0.2,0.4,0.6,0.8,1.0),(r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}",r"\bf{0.8}",r"\bf{1}")) +def axdemiabsticks(ax): + ax.set_xticks((0.0,0.2,0.4,0.6,0.8,1.0)) + ax.set_xticklabels((r"\bf{0}",r"\bf{0.2}",r"\bf{0.4}",r"\bf{0.6}",r"\bf{0.8}",r"\bf{1}")) + + + +################ Mardi 4 juin : illustrons sur l'exemple de sylvain laconvergence de l'ode de Bacaer vers la prob d'emergence + +def tbacasyl(T=100,lzero=0.2,muzero=0.1,gamma=0.3,nbperiodes=3): + def mu(t): + return(muzero) + def Mu(t): + return(muzero*t) + #il faut ajuster le lzerosin pour que les processus aient le meme rzero + + def lastep(t): + if (t< 1-gamma): + return(lzero) + else: + return(0.0) + def Lastep(t): + if (t< 1-gamma): + return(lzero*t) + else: + return(lzero*(1-gamma)) + def guess(t): + a=lastep(t)/mu(t) + if (a>1): + return(1-(1/a)) + else: + return(0.0) + plastep=periodise(lastep) + pmu=periodise(mu) + #nbpts=nb*tau + Lambda=[[plastep]] + VMu=[pmu] + #z=vfbaca(Lambda,VMu,T=T,tv=False) + fig, ax = plt.subplots(2,1,figsize=[16, 12]) + plt.tight_layout() + + plt.subplot(2,1,1) + z=nvfbaca(Lambda,VMu,T=T,tv=False,nbtau=nbperiodes)#fait en plus un plot + z=z.reshape(len(z),) + pem=nncalculepe(lastep,Lastep,mu,Mu,T) + t=np.linspace(0.0,1.0,len(z)) + valpem=np.array([pem(s) for s in t]) + vg=np.array([guess(s) for s in t]) + + plt.subplot(2,1,2) + #plt.title("Square wave birth rate") + plt.ylabel(r"${p_e(t_0 T,T)}$",fontsize=22) + plt.plot(t,z,label="simulation de Bacaer " +str(nbperiodes)+" periodes",color=couleurcontrole) + plt.plot(t,valpem,label="calcul integral exact",dashes=[3,3],color='black') + plt.xlabel(r"\textbf{Introduction time} ${t_0}$",fontsize=22) + #pas besoin de mettre le guess je pense + #plt.plot(t,vg,label=r"guess : $(1-\lambda/\mu)^+$",linestyle=':') + #plt.legend() + plt.savefig("approximationbacaer1d.pdf",dpi=300) + print("difference max",np.abs(valpem-z).max()) + #on n'a pas zero car les tableaux n'ont pas meme dimension + return((z,valpem)) + + +def nvfbaca(Lambda,Mu,T,tv=False,nbtau=10,toute=True): + r""" la on essaie de plotter toute la solution de l'ode de bacaer""" + tau=nbtau*T + nbpts=10*tau + d=len(Mu) + zzero=np.array([1.0 for i in range(d)]) + + def msisi(x,t): + #s=tau-(t/T) + s=(tau-t)/T + return(np.array([-x[i]*Mu[i](s) + (1-x[i])*((np.array([Lambda[i][j](s)*x[j] for j in range(d)])).sum()) for i in range(d)])) + + timeint=np.arange(0,tau+1/nbpts,tau/nbpts) + z=np.array(odeint(msisi,zzero,timeint)) + if toute: + plt.ylabel(r"${Y^{(\tau)}}$",fontsize=22) + plt.plot(timeint,z,label=r"$Y^{(\tau)}$",color=couleurcontrole) + plt.xlabel(r"\textbf{Approximation time : } ${s=\tau-t_0}$",fontsize=22) + i=int(T*nbpts/tau) + #print("T=",T,"tau=",tau,"z.shape",z.shape,"i=",i) + res=z[-i:] + if tv: + plt.plot(timeint[-2*i:],z[-2*i:],label="z") + plt.legend() + return(res[::-1]) +