first commit

simulatorCOV.py

+
+import numpy as np
+import scipy.stats as st
+from collections import Counter
+import random
+import matplotlib.pyplot as plt
+import time
+import pandas as pd
+
+sane,exp,synt,asynt,rec,vax = 0,1,2,3,4,5
+
+n_edgs_adn_adult,n_edgs_adn_child = 3,7
+gamma,xmin = -2.09,0.1
+
+n_family = 50000
+
+mean_n_class = 20
+child_per_class = 24
+I0 = 0.01
+m0 = 18
+
+
+TE = 6.4
+TI = 5
+
+
+mu = 1-np.exp(-1/TI)
+
+
+def power_law_activity(rd,N,alpha,xmin=0.01):
+    xmax = 1.
+
+    r = rd.uniform(0,1,size=N)
+    ax = np.power((np.power(xmax,alpha+1)-np.power(xmin,alpha+1))*r+np.power(xmin,alpha+1),(1./(alpha+1)))
+    return ax
+
+
+def IsoFamily(status_adult,status_child,people_family,to_fam):
+    n_adult,n_child = status_adult.shape[0],status_child.shape[0]
+
+    avail = np.zeros(status_adult.shape[0]+status_child.shape[0]).astype(bool)
+
+    iso_fam = np.unique(people_family[np.concatenate((status_adult==synt,status_child==synt))])
+
+    z = to_fam[iso_fam]
+    if len(z)>0:
+        avail[np.concatenate(z)] = True
+
+    avail_adult,avail_child = ~avail[:n_adult],~avail[n_adult:]
+    return avail_adult,avail_child
+
+def IsoClassFamily(status_adult,status_child,people_family,child_class,to_fam,to_class):
+    avail_adult,avail_child0 = IsoFamily(status_adult,status_child,people_family,to_fam)
+
+    n_child = status_child.shape[0]
+
+    avail_child = np.zeros(status_child.shape[0]).astype(bool)
+
+    iso_class = np.unique(child_class[status_child==synt])
+
+    z = to_class[iso_class]
+    if len(z)>0:
+        avail_child[np.concatenate(z)] = True
+
+
+    avail_child = (~avail_child)*avail_child0
+    return avail_adult,avail_child
+    
+def Select_Avail(status_adult,status_child,people_family,child_class,typ,to_fam,to_class):
+    
+    if typ=='family':
+        return IsoFamily(status_adult,status_child,people_family,to_fam)
+    if typ=='class':
+        return IsoClassFamily(status_adult,status_child,people_family,child_class,to_fam,to_class)
+    
+def Family(n_fam):
+    s = [1]*int(0.358*n_fam) #Singoli
+
+    cf = 3+st.poisson.rvs(0.79,size=int(n_fam*0.215)) # Coppie con minori
+    sf = 2+st.poisson.rvs(0.79,size=int(n_fam*0.06)) # Singoli con minori
+    ath = 2+st.poisson.rvs(0.41,size= int( n_fam*(1-(0.358+0.215+0.06))) ) #Altri adulti>2
+
+    n_child = (cf-2).sum()+(sf-1).sum()
+    n_adult = len(s)+2*len(cf)+len(sf)+(ath).sum()
+
+    single_adult = np.arange(len(s))
+
+    cp_adult = single_adult.max()+1+np.repeat( np.arange(cf.shape[0]),2)
+
+    cp_child = single_adult.max()+1+np.array([i for i,n in enumerate(cf-2) for _ in range(n)])
+
+    sf_adult = cp_adult.max()+1+np.arange(sf.shape[0])
+
+    sf_child = cp_adult.max()+1+np.array([i for i,n in enumerate(sf-1) for _ in range(n)])
+
+    ath_adult = sf_adult.max()+1+np.array([i for i,n in enumerate(ath) for _ in range(n)])
+
+    people_familiy = np.concatenate((single_adult,cp_adult,sf_adult,ath_adult,cp_child,sf_child))
+    
+    size_family = np.bincount( people_familiy )
+    
+    to_fam = np.column_stack((people_familiy,np.arange(n_adult+n_child)))
+    to_fam = to_fam[np.argsort(to_fam[:,0])]
+    to_fam = np.array(np.split(to_fam[:,1], np.unique(to_fam[:,0], return_index = True)[1])[1:])
+
+    
+    return n_adult,n_child,people_familiy,size_family,to_fam
+
+def CreateSchool(n_child,mean_n_class,child_per_class):
+    n_class = int(n_child/child_per_class)
+
+    n_school = int(round(n_child/(child_per_class*mean_n_class),0))
+
+    indx_chilid = np.arange(n_child).astype(int)
+
+    school = np.array([np.array_split(h,mean_n_class) for h in np.array_split(indx_chilid,n_school)])
+
+    possible_adn_link_child = [list(np.concatenate(np.delete(sc,i))) for sc in school for i in range(len(sc))]
+
+    child_class = np.concatenate([np.ones(len(s))*i for i,s in 
+                                  enumerate(np.concatenate(school))]).astype(int)
+    
+    n_class = len(possible_adn_link_child)
+    
+    
+    size_class = np.bincount( child_class )
+    
+    
+    to_class = np.column_stack((child_class,np.arange(n_child)))
+    to_class = to_class[np.argsort(to_class[:,0])]
+    to_class = np.array(np.split(to_class[:,1], np.unique(to_class[:,0], return_index = True)[1])[1:])
+
+    
+    return child_class,possible_adn_link_child,n_class,size_class, to_class
+
+def CongionClass(rd,status_child,child_sane_av,child_class,child_inf_av,lmd_child,n_class):
+    inf_class = np.bincount( child_class[child_inf_av] )
+
+    #probabiliy of infection per classe
+    p_class = 1-(1-lmd_child)**inf_class
+    p_class = np.pad(p_class,(0,n_class-p_class.shape[0]))
+    
+    infx = ((status_child[child_sane_av]==exp)+rd.binomial( 1,p_class[child_class[child_sane_av]] ).astype(bool) )*exp
+    
+    return infx
+
+def ADN_Child(rd,child_sane_av,child_inf_av,A_child,possible_adn_link_child,child_class,n_edgs_adn_child):
+    
+    active_child = np.where((child_sane_av+child_inf_av)*rd.binomial(1,A_child))[0]
+    if len(active_child)>0:
+        #Crea links
+        adnlink = np.array([random.sample(possible_adn_link_child[i],n_edgs_adn_child) 
+                            for i in child_class[active_child]])
+
+        origin = (np.tile(active_child,(n_edgs_adn_child,1)).T).flatten()
+        dest = adnlink.flatten()
+
+        #Contagia quelli a rischio
+        risk_cont = np.concatenate( (origin[np.where(child_sane_av[origin]*child_inf_av[dest])[0]],
+                       dest[np.where(child_sane_av[dest]*child_inf_av[origin])[0]] ))
+
+        idx,n = np.unique( risk_cont , return_counts=True)
+
+        return idx,n
+        
+    else:
+        return np.array([]),[]
+
+
+def CheckStatus(x,avail):
+    x_inf = ((x==synt)+(x==asynt))
+    x_exp = (x==exp)
+    x_sane = (x==sane)
+
+    x_inf_av = (x_inf*avail).astype(bool)
+    x_exp_av = (x_exp*avail).astype(bool)
+    x_sane_av = (x_sane*avail).astype(bool)
+    
+    return x_inf,x_exp,x_sane,x_inf_av,x_exp_av,x_sane_av
+
+
+def CongionFamility(rd,people_family, adult_inf, child_inf, n_family,adult_sane,child_sane,
+                    lmd_child,lmd_adult,status_adult,status_child):
+    n_child,n_adult = child_inf.shape[0],adult_inf.shape[0]
+
+    inf_adult_familiy = np.bincount(people_family[:n_adult][adult_inf])
+    inf_adult_familiy = np.pad(inf_adult_familiy,(0,n_family-inf_adult_familiy.shape[0]))
+
+    inf_child_familiy = np.bincount(people_family[n_adult:][child_inf])
+    inf_child_familiy = np.pad(inf_child_familiy,(0,n_family-inf_child_familiy.shape[0]))
+
+
+    p_family = 1-(1-lmd_adult)**inf_adult_familiy*(1-lmd_child)**inf_child_familiy
+    #p_child_family = 1-(1-lmd_child)**(inf_adult_familiy+inf_child_familiy)
+    
+    
+    ad_to_inf = ((status_adult[adult_sane]==exp)+rd.binomial(1,p_family[
+        people_family[:n_adult][adult_sane]]).astype(bool))*exp
+    
+    ch_to_inf = ((status_child[child_sane]==exp)+rd.binomial(1,p_family[
+        people_family[n_adult:][child_sane]]).astype(bool))*exp
+    
+    return ad_to_inf,ch_to_inf
+
+
+def ADN_Adult(rd,adult_exp_av,adult_sane_av,adult_inf_av,A_adult,possible_adult_link,people_family,n_edgs_adn_adult):
+    n_adult = adult_exp_av.shape[0]
+
+    # Contagio ADN_adulti
+    active_adult = np.where((adult_sane_av+adult_inf_av)*rd.binomial(1,A_adult))[0]
+    if len(active_adult)>0:
+        #Crea links
+        adnlink = np.array([random.sample(possible_adult_link,n_edgs_adn_adult)
+                            for i in people_family[:n_adult][active_adult]])
+
+
+        origin = (np.tile(active_adult,(n_edgs_adn_adult,1)).T).flatten()
+        dest = adnlink.flatten()
+
+
+        msk = (people_family[origin]!=people_family[dest])
+
+        origin,dest = origin[msk],dest[msk]
+
+        #Contagia quelli a rischio
+        risk_cont = np.concatenate( (origin[np.where(adult_sane_av[origin]*adult_inf_av[dest])[0]],
+                       dest[np.where(adult_sane_av[dest]*adult_inf_av[origin])[0]] ))
+
+        idx,n = np.unique( risk_cont , return_counts=True)
+        return idx,n
+        
+    else:
+        return np.array([]),[]
+
+def Update(rd,status_x,x_exp,x_inf,vd,vn,mu):
+    to_inf = rd.binomial(1, vd+vn , size=x_exp.sum()).astype(bool)
+
+    to_det = rd.binomial(1,vd/(vd+vn),size=to_inf.sum()).astype(bool)
+
+    whoexp = np.where(x_exp)[0]
+
+
+    status_x[whoexp[to_inf][to_det]] = synt
+    status_x[whoexp[to_inf][~to_det]] = asynt
+
+
+
+    whoinf = np.where(x_inf)[0]
+
+
+    to_rec = rd.binomial(1, mu , size=x_inf.sum()).astype(bool)
+
+    status_x[whoinf[to_rec]] = rec
+    
+    return 
+
+def GET_lmdPart(R0,size_family,people_family,n_adult,child_class,n_child,d):
+    from scipy.integrate import quad
+    w = quad(lambda x: x**gamma,xmin,1)[0]
+    Am = quad(lambda x: x*(x**gamma),xmin,1)[0]/w
+
+
+    k_adult_fam = (size_family[people_family[:n_adult]]-1).mean()
+    k_child_fam = (size_family[people_family[n_adult:]]-1).mean()
+
+
+    n_edgs_adn_adult = int(round((m0-k_adult_fam)/(1+Am),0))
+
+    k_adult = 2*Am*n_edgs_adn_adult + k_adult_fam
+
+
+    k_child_class = (np.unique(child_class,return_counts=True)[1]-1).mean()*(5/7)#weektime
+    k_child = 2*Am*n_edgs_adn_child + k_child_class
+
+
+    wc,wa = (n_child/(n_child+n_adult)),(n_adult/(n_child+n_adult))
+
+    kmean = (k_adult*wa + d*k_child*wc)
+    adn_full = int( round(kmean/(2*Am),0) )
+
+    lmd_adult = R0/( TI* kmean)
+    lmd_child = d*lmd_adult
+    
+    return lmd_adult,lmd_child,n_edgs_adn_adult,adn_full
+
+def GET_SERIES(Sa,Sc):
+
+    Qc = pd.DataFrame.from_dict(Sc)
+    Qa = pd.DataFrame.from_dict(Sa)
+
+    if not vax in Qc.columns:
+        Qc[vax]=0
+
+    if not vax in Qa.columns:
+        Qa[vax]=0
+
+    Qa[Qa.isna()] = 0
+    Qc[Qc.isna()] = 0
+
+    Qa = Qa[[sane,exp,synt,asynt,rec,vax]].astype(int)
+    Qa.columns=['sane','exp','synt','asynt','rec','vax']
+
+    Qc = Qc[[sane,exp,synt,asynt,rec,vax]].astype(int)
+    Qc.columns=['sane','exp','synt','asynt','rec','vax']
+    
+    return Qa,Qc
+
+
+
+def RUN(x):
+
+	seed,qA,qC,R0,pvax,d,vaxstrat,policy = x
+
+	rd = np.random.RandomState(seed)
+
+	vdA = qA*(1-np.exp(-1/TE))
+	vnA = (1-qA)*(1-np.exp(-1/TE))
+
+	vdC = qC*(1-np.exp(-1/TE))
+	vnC = (1-qC)*(1-np.exp(-1/TE))
+
+	for _ in range(100):
+		try:
+			n_adult,n_child,people_family,size_family,to_fam = Family(n_family)
+
+			A_child = power_law_activity(rd,n_child,gamma,xmin)
+			A_adult = power_law_activity(rd,n_adult,gamma,xmin)
+
+			child_class,possible_adn_link_child,n_class,size_class,to_class =  CreateSchool(n_child,mean_n_class,
+																							child_per_class)
+
+			possible_adult_link = list(np.arange(n_adult))
+
+
+			lmd_adult,lmd_child,n_edgs_adn_adult,adn_full =  GET_lmdPart(R0,size_family,
+																		 people_family,n_adult,
+																		 child_class,n_child,d)
+		except:
+			continue
+		break
+
+
+
+	nvax = int(pvax*(n_adult+n_child))
+
+
+
+	inf_class = np.zeros(n_class,dtype='int')
+
+	status_adult = np.zeros(n_adult).astype(int)
+	status_child = np.zeros(n_child).astype(int)
+
+	if nvax>0:
+
+		if vaxstrat=='random_adult':
+			status_adult[np.array(random.sample(range(n_adult),nvax))] = vax
+		if vaxstrat=='priority':
+			status_adult[np.argsort(size_family[people_family[:n_adult]])[-nvax:]] = vax
+		if vaxstrat=='All':
+			vax_adult = int(rd.binomial(nvax,n_adult/(n_adult+n_child)))
+			vax_child = int(nvax -vax_adult)
+			
+			status_adult[np.array(random.sample(list(range(n_adult)),vax_adult))] = vax
+			status_child[np.array(random.sample(list(range(n_child)),vax_child))] = vax
+
+
+	status_adult[(status_adult!=vax)] = rd.binomial(1,I0,size=n_adult-(status_adult==vax).sum())
+	status_child[(status_child!=vax)] = rd.binomial(1,I0,size=n_child-(status_child==vax).sum())
+
+
+	Sc,Sa = [],[]
+
+	t = 0
+	while True:
+
+		Sc.append( Counter(status_child) )
+		Sa.append( Counter(status_adult) )
+		if Sc[-1][exp]+Sc[-1][asynt]+Sc[-1][synt]+Sa[-1][exp]+Sa[-1][asynt]+Sa[-1][synt]==0:
+			break
+
+		avail_adult,avail_child = Select_Avail(status_adult,status_child,people_family,
+											   child_class,policy,to_fam,to_class)
+
+		child_inf,child_exp,child_sane,child_inf_av,child_exp_av,child_sane_av =  CheckStatus(status_child,
+																							  avail_child)
+
+		adult_inf,adult_exp,adult_sane,adult_inf_av,adult_exp_av,adult_sane_av =  CheckStatus(status_adult,
+																							  avail_adult)
+
+
+		if t%7<5:
+			status_child[child_sane_av] = CongionClass(rd,status_child,child_sane_av,
+													   child_class,child_inf_av,lmd_child,n_class)
+
+			idx,n = ADN_Child(rd,child_sane_av,child_inf_av,A_child,
+							  possible_adn_link_child,child_class,n_edgs_adn_child)
+			if len(n)>0:
+				status_child[idx[rd.binomial(1,1-(1-lmd_child)**n).astype(bool)]] = exp
+
+		status_adult[adult_sane],status_child[child_sane] =  CongionFamility(rd,people_family,
+																			 adult_inf, child_inf, 
+																			 n_family,adult_sane,child_sane,
+																			lmd_child,lmd_adult,
+																			 status_adult,status_child)
+
+
+
+		idx,n = ADN_Adult(rd,adult_exp_av,adult_sane_av,adult_inf_av,A_adult,possible_adult_link,people_family,n_edgs_adn_adult)
+		if len(n)>0:
+			status_adult[idx[rd.binomial(1,1-(1-lmd_adult)**n).astype(bool)]] = exp
+
+
+		Update(rd,status_child,child_exp,child_inf,vdC,vnC,mu)
+
+		Update(rd,status_adult,adult_exp,adult_inf,vdA,vnA,mu)
+		t+=1
+
+
+
+	Qa,Qc = GET_SERIES(Sa,Sc)
+	return Qa,Qc