kupkasmale/WebWorld
1
1import panel as pn2import numpy as np3import pandas as pd4import hvplot.pandas5from numba import jit6import holoviews as hv7from functools import partial8 9import hvplot.networkx as hvnx10import networkx as nx11 12 13@jit(nopython = True)14def gaussianC(C):15 16 if C != 1.0:17 Cj = np.random.uniform()18 if C<Cj:19 y = np.random.normal()20 else:21 y = 0.022 return y23 24 else:25 26 return np.random.normal()27 return y 28 29#@pn.cache()30def fixed_parameters():31 32 N = 60 #0 33 R = 1e6 #134 c = 0.2 #235 ld = 0.1 #336 Ng = 1e4 #437 L = 10 #538 K = 500 #639 b = 5e-3 #740 41 return np.array([N,R,c,ld,Ng,L,K,b])42 43def start_sim(parvalues,slides,button,tabs,revent):44 button.disabled = True45 46 #print(SN)47 Nsp = int(parvalues[0])48 R = slides[0].value #parvalues[1]49 c = slides[1].value #parvalues[2]50 ld = parvalues[3]51 Ngen = int(parvalues[4])52 L = int(parvalues[5])53 K = int(parvalues[6])54 b = parvalues[7]55 56 #print(c,R)57 58 TR = -1*np.ones((Nsp,L),dtype=np.int32)59 M = Mab(K);60 #print(K,L)61 G = np.zeros((Nsp, Nsp)); #Rates62 S = np.zeros((Nsp, Nsp)); #Scores63 A = np.zeros((Nsp, Nsp)); #Competition64 F = np.zeros((Nsp, Nsp)); #Foraging Efforts65 N = np.zeros(Nsp); #Populations66 67 #Initial Population: No68 #Minimum population allowed: Nd69 No = 1.0; Nd = 1.0;70 71 #Initial Setup:72 prm = np.array([K,L,R,ld,b]);73 so, s1, sio, fio = initial_setting(prm,M);74 75 N[0] = R/ld; N[1] = No;76 S[1,0] = sio; F[1,0] = fio; A[1,1] = 1.0;77 TR[0] = so; TR[1] =s1;78 79 dF = 1.0; DN = 1.0; DT = 0.2;80 81 SN = pd.DataFrame(columns=['SN'])82 Snumb = []83 #Stp = []84 #SN['time'] = Stp85 SN['SN'] = Snumb86 87 figa = (88 SN['SN']89 ).hvplot(90 title="Number of Species",91 ylabel="Number of Species",92 xlabel="Evolutionary event",93 autorange = "y",94 color="indigo",95 )96 plotns = pn.pane.HoloViews(figa, sizing_mode="stretch_both", name="PlotNS")97 tabs[0] = plotns98 99 100 101 for k in range(Ngen):102 103 while DN>1e-2:104 while dF>0.1:105 F,G,dF = foraging_strategy(N,F,A,S,b)106 N, chk, sj, DN = pop_update(N,G,ld,DT,Nd)107 if chk == True:108 N,S,F,A,TR = pop_check(N,S,F,A,TR,sj)109 dF = 1.0110 DN=1.0111 112 Fq = F113 Nq = N114 115 F,S,A,TR,N = speciation(M,TR,A,S,F,K,L,c,No,N);116 117 NA = nspecies(N)118 parvalues[8].value = NA119 parvalues[9].value = k120 Snumb.append(NA)121 #Stp.append(k)122 123 if k%100 == 0:124 #SN['SN'] = []125 #SN['SN'] = Snumb126 figa = (127 pd.DataFrame(Snumb,columns=['SN'])128 ).hvplot(129 title="Number of Species",130 ylabel="Number of Species",131 xlabel="Evolutionary event",132 autorange = "y",133 color="indigo",134 )135 plotns = pn.pane.HoloViews(figa, sizing_mode="stretch_both", name="PlotNS")136 tabs[0] = plotns137 138 if NA == Nsp-1:139 break140 141 figa = (142 pd.DataFrame(Snumb,columns=['SN'])143 ).hvplot(144 title="Number of Species",145 ylabel="Number of Species",146 xlabel="Evolutionary event",147 #xlim=(0, Ng),148 autorange = "y",149 color="indigo",150 )151 plotns = pn.pane.HoloViews(figa, sizing_mode="stretch_both", name="PlotNS")152 tabs[0] = plotns153 154 FwG, pos,sz,cn = get_foodweb(Fq,Nq)155 #(Fq,Nq)156 nodes = hvnx.draw_networkx_nodes(FwG, pos,with_labels=True,node_size=sz,157 node_color=cn)158 egx = FwG.edges()159 weights = [FwG[u][v]['edge_width'] for u,v in egx ]160 print(weights)161 162 edges = hvnx.draw_networkx_edges(FwG, pos, node_size=sz,edge_color='red',163 arrowstyle="->",arrowsize=1,alpha=0.3, edge_width=weights)164 165 figb = nodes * edges166 plotfw = pn.pane.HoloViews(figb, sizing_mode="stretch_both", name="PlotFW")167 tabs[1] = plotfw168 169 button.disabled = False170 171def Mab(K):172 173 A=np.ones((K,K))174 175 Iup = np.triu_indices(K);176 Idn = np.tril_indices(K);177 178 N = len(Iup[0])179 U = np.array([gaussianC(1.0) for j in range(N)])180 181 A[Iup] = U182 for i,j in zip(Idn[0],Idn[1]):183 A[i,j] = - A[j,i]184 185 np.fill_diagonal(A,0)186 return A187 188 189@jit(nopython = True)190def initial_setting(prm,M):191 192 K = prm[0]; L = prm[1]; R = prm[2];193 ld = prm[3]; b = prm[4];194 S1o = 0.0;195 bo = b/ld;196 while (S1o <= bo):197 so = gen_string(L,K);198 s1 = gen_string(L,K);199 S1o = Sij(so,s1,M)/L;200 f1o = 1.0; 201 202 return so, s1, S1o, f1o203 204@jit(nopython = True)205def gen_string(L,K):206 207 u = K*np.ones(int(L))208 l = 0209 210 while l<L:211 nj = np.random.randint(0,int(K))212 if nj not in u:213 u[l]=nj214 l+=1215 return u.astype(np.int32)216 217 218@jit(nopython = True)219def Sij(a,b,M):220 221 Sab = 0.0222 for ia in a:223 S = M[ia,b]224 Sab+=np.sum(S)225 226 if Sab > 0.0:227 return Sab228 229 else:230 return 0.0231 232@jit(nopython = True)233def intersect(a,b):234 235 u = np.intersect1d(a,b);236 d = len(u)/len(a)237 return d238 239@jit(nopython = True)240def comp_score(qij,c):241 return c+(1-c)*qij242 243 244@jit(nopython = True)245def foraging_strategy(N,F,A,S,b):246 247 Fn = np.zeros(F.shape)248 G = np.zeros(F.shape)249 nk = np.where(N>0)[0]250 jn = np.max(nk)251 252 for i in nk:253 if(i>0):254 Sq = S[i,:]255 Fq = F[i,:]256 pk = np.where(Sq>0)[0]257 258 for j in pk:259 Sk = np.where(S[:,j]>0)[0]260 261 sn = 0.0;262 for k in Sk:263 sn+=A[k,i]*S[k,j]*F[k,j]*N[k]264 G[i,j]=Sq[j]*Fq[j]*N[j]/(b*N[j]+sn)265 266 for i in nk:267 if (i>0):268 SGij = np.sum(G[i,:])269 Fn[i,:] = G[i,:]/SGij270 jn = np.where(Fn[i,:]>0.0)271 for j in jn[0]:272 if Fn[i,j]<1e-6:273 Fn[i,j] = 1e-6274 275 276 DF = np.max(np.abs((F-Fn).ravel()))277 278 return Fn, G, DF279 280 281@jit(nopython = True)282def pop_update(N,G,ld,DT,Nd):283 284 rem = False;285 286 Np = np.zeros(N.shape)287 la = (1 - DT)*N;288 lb = ld*DT*np.sum(G,axis=1)*N289 V = np.zeros(N.shape)290 291 nk = np.where(N>0)[0]292 293 for j in nk:294 if j>0:295 V[j]=np.sum(G[:,j]*N)296 297 Np = la + lb - DT*V298 Np[0] = N[0]299 300 U = Np[(Np!=0) & (Np<Nd)]301 if len(U)>0:302 sj = np.where(((Np!=0) & (Np<Nd)))[0]303 rem = True304 Np[sj] = 0.0305 306 DN = np.max(np.abs(N-Np))307 308 return Np,rem, sj, DN309 310@jit(nopython = True)311def pop_check(N,S,F,A,TR,jn):312 Ln = len(TR[0,:])313 for jd in jn:314 N[jd] = 0.0315 S[jd,:] = 0.0316 S[:,jd] = 0.0317 F[jd,:] = 0.0318 F[:,jd] = 0.0319 A[jd,:] = 0.0320 A[:,jd] = 0.0321 TR[jd,:] = -1*np.ones(Ln,dtype=np.int32)322 323 return N,S,F,A,TR 324 325 326@jit(nopython = True)327def speciation(M,TR,A,S,F,K,L,c,No,N):328 Tn = TR[:,0]329 330 nsj = np.min(np.where(Tn==-1)[0]) #Allocated new species 331 Na = np.where(Tn>-1)[0]332 Nb = Na[Na!=0] #excludes external species333 ns = np.random.choice(Nb) #parent species334 maxj = np.max(Nb) #index of largest non zero species335 336 snew = False; NFT = False;337 new_traits = TR[ns,:].copy()338 339 while(snew==False):340 while(NFT == False):341 nf = np.random.randint(0,K)342 if nf not in new_traits:343 jn = np.random.randint(0,L)344 new_traits[jn] = nf345 NFT = True346 n=0347 for q in Na:348 a = TR[q]349 dij = np.intersect1d(a,new_traits)350 lj = len(dij)351 if (lj == L):352 n+=1353 354 if n==0:355 snew = True356 else:357 NFT = False358 359 ##Initialising scores and pop360 TR[nsj] = new_traits361 N[nsj] = No362 N[ns] -= 1 363 A[nsj][nsj] = 1.0364 365 ia = TR[nsj]366 367 for q in Na:368 ja = TR[q]369 370 A[nsj][q] = comp_score(intersect(ia,ja),c)371 A[q][nsj] = A[nsj][q]372 373 if q == 0:374 375 S[nsj][0] = Sij(ia,ja,M)376 if S[nsj][0] > 0.0:377 378 if F[ns][0] > 0.0:379 F[nsj][0] = F[ns][0]380 else:381 F[nsj][0] = np.random.uniform(1e-6,1)382 383 if q>0:384 S[nsj][q] = Sij(ia,ja,M)385 S[q][nsj] = Sij(ja,ia,M)386 387 if S[nsj][q] > 0.0:388 if F[ns][q] > 0.0:389 F[nsj][q] = F[ns][q]390 else:391 F[nsj][q] = np.random.uniform(1e-6,1)392 393 if S[q][nsj]>0.0:394 if F[q][ns] > 0.0:395 F[q][nsj] = F[q][ns]396 else:397 F[q][nsj] = np.random.uniform(1e-6,1)398 399 return F,S,A,TR,N400 401 402@jit(nopython = True)403def nspecies(N):404 405 Na = np.where(N>0)[0]406 #ns =np.max(Na)407 NA= len(Na)408 409 return NA-1410 411def get_foodweb(F,N):412 413 F[F<=0.01] = 0.0414 #Positive Population indices415 jall = np.where(N>0)[0]416 #All positive pop indices without the external species 417 jy = jall[jall>0] 418 #Basal Species indices419 jo = np.where(F[:,0]>0.0)[0] 420 tlvs = []421 tlvs.append(jo)422 spn = False423 M=0424 #Levels425 while spn == False:426 jx = np.setdiff1d(jy,jo) #Speciesn in jy but not in jo427 ln=[]428 for q in jx:429 Un = np.where(F[q,:]>0.0)[0]430 for p in jo:431 if p in Un:432 ln.append(q)433 break434 ln=list(set(ln))435 tlvs.append(ln)436 un = list(jo)+ln437 jo = np.asarray(un)438 M+=1439 440 if M ==5:441 spn = True442 443 DN = []444 DR = []445 DRL = []446 for q in tlvs:447 if len(q)>0:448 Plog = np.log(1+N[q])449 Nlv = len(Plog)450 DN.append(Plog)451 DR.append(np.max(Plog))452 DRL.append(Nlv)453 454 Dy = np.max(DR)455 456 Xn = []; Yn = [];457 for q in range(len(DN)):458 X = [(6*j+1)*DR[q]-DR[q] for j in range(DRL[q])]459 Y = [(q+1)*8*Dy for j in range(DRL[q])]460 Xn.append(X);Yn.append(Y)461 462 Lmx = [] 463 for q in Xn: 464 Lmx.append(q[len(q)-1])465 466 lj = np.max(Lmx)467 468 for q in range(len(Xn)):469 if len(Xn[q]) > 1:470 Xn[q] = [lj/Lmx[q]*k for k in Xn[q]]471 472 FW = pd.DataFrame(columns=['n','X','Y','R'])473 sn = [];xn = []; yn = []; rn = []474 for q in range(len(tlvs)):475 for p in range(len(tlvs[q])):476 sn.append(tlvs[q][p])477 xn.append(Xn[q][p])478 yn.append(Yn[q][p])479 rn.append(DN[q][p])480 481 FW['n'] = sn; FW['X'] = xn; FW['Y'] = yn; FW['R'] = rn482 483 tlevs = []484 485 for q in tlvs:486 if len(q)>0:487 tlevs.append(q)488 489 490 G = nx.DiGraph()491 pos = {}492 n=0493 sizes = []494 colormap = []495 edgesw = []496 497 for q in tlevs:498 if n == 0:499 basal = FW['n'].isin(q)500 FBAS = FW[basal]501 RY = np.max(list(FBAS['R']))502 503 for i in q:504 505 nprop = FBAS.loc[FBAS['n']==i]506 rn = nprop['R'].values[0]507 xn = nprop['X'].values[0]508 yn = nprop['Y'].values[0]509 510 a = 'er'+str(i)511 b = str(i)512 G.add_node(a)513 sizes.append(100.0)514 pos[a] = (xn,0.0)515 colormap.append('yellow')516 G.add_node(b)517 sizes.append(30*rn)518 colormap.append('blue') 519 pos[b] = (xn,yn)520 G.add_edge(a,b,edge_width = 5*F[i,0])521 n+=1522 else:523 524 TRL = FW['n'].isin(q)525 FRL = FW[TRL]526 527 for i in q:528 a = str(i)529 nprop = FRL.loc[FRL['n']==i]530 rn = nprop['R'].values[0]531 xn = nprop['X'].values[0]532 yn = nprop['Y'].values[0]533 G.add_node(a)534 sizes.append(30*rn)535 colormap.append('blue') 536 pos[a] = (xn, yn)537 538 for q in tlevs:539 for i in q:540 jn = np.where(F[i,:]>0)[0]541 if len(jn)>0:542 for j in jn :543 if j!=0:544 a=str(j);b=str(i)545 G.add_edge(a,b,edge_width=5*F[i,j])546 547 return G,pos,sizes,colormap