Bird Simulation Evaluation Script

Imports & Preparations



In [1]:

    
import numpy as np
import scipy as sp
import birds
import argparse
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.animation import FuncAnimation
from matplotlib.collections import PathCollection
from IPython.display import HTML
from scipy.optimize import curve_fit
#%matplotlib ipympl
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from matplotlib import rcParams
rcParams.update({'figure.autolayout': True})

Figure output Settings



In [2]:

    
figpath = '../img/'
figwidth = 4 #figure width in inches
figsize = (figwidth,figwidth*2.5/4)

Run with default Settings



In [ ]:

    
frames = 1000
birds.param_record = False
birds.trace = None
birds.flock = birds.Flock()
fig = plt.figure(figsize=(5, 5*birds.flock.args['height']/birds.flock.args['width']), facecolor="white")
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], aspect=1, frameon=False)
birds.collection = birds.MarkerCollection(birds.flock.args['n'])
ax.add_collection(birds.collection._collection)
ax.set_xlim(0, birds.flock.args['width'])
ax.set_ylim(0, birds.flock.args['height'])
ax.set_xticks([])
ax.set_yticks([])

animation = FuncAnimation(fig, birds.update, interval=10, frames=frames)
HTML(animation.to_html5_video())

Find moving phase

Run with varying Eta



In [5]:

    
def avg_with_error(f, avg_time, va_t = False, eta_index = 0, prerun_time = 500):
    var = np.zeros(avg_time)
    va_avg = 0
    for t in range(avg_time):
        va_tmp = f.get_va()
        va_avg += va_tmp
        var[t] = va_tmp
        f.run()
    va_avg /= avg_time
    if va_t is not False:
        va_t[eta_index, prerun_time:] = var
    var = np.sum((var-va_avg)**2 / avg_time)
    return va_avg, var



In [6]:

    
res = 30
prerun_time = 500
averaging_time = 1000
repeat = 3
rho=4
Eta = np.linspace(0.,7.,res)
N = [50,100,400]

va_t = np.zeros((6, prerun_time+averaging_time)) 

va = np.zeros((len(N),res))
vas = np.zeros(repeat)
errorbars = np.zeros_like(va)
variance = np.zeros(repeat)
for c,n in enumerate(N):
    for i,eta in enumerate(Eta):
        for j in range(repeat):
            f = birds.Flock(n=n,eta=eta,rho=rho)
            record_time = (n==100 and i%5==0 and j==0)
            for t in range(prerun_time):
                f.run()
                if record_time:
                    va_t[int(i/5),t]=f.get_va()
                    
            if record_time:
                va_avg, vari = avg_with_error(f, averaging_time, va_t, int(i/5), prerun_time)
            else:
                va_avg, vari = avg_with_error(f, averaging_time)
            vas[j] = va_avg
            variance[j] = vari
        va[c][i] = vas.sum()/repeat
        errorbars[c][i] = np.sqrt(variance.sum()/repeat)
        
print('Has been run.')









    



Has been run.



In [26]:

    
plt.figure(figsize=(10,8))
x = np.linspace(0,1500,1500)
for e,vt in enumerate(va_t):
    plt.plot(x, vt, label='Eta = '+ '%1.2f' % Eta[e*5])
plt.legend()
plt.xlim([0, 40])
plt.xlabel('$t$')
plt.ylabel('$v_a$')
plt.show()
len(x)









    



/usr/lib/python3.6/site-packages/matplotlib/figure.py:1999: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "






    












    Out[26]:





1500

Fit $\eta_c$



In [ ]:

    
prec = 0.05

# Initial square/lin fit to determinine p0-parameters
def finvsquare(x, sqwidth):
    return 1-(x*sqwidth)**2
def flin(x, m, b):
    return m*x+b

# find array index of eta_c as first approx to split fit in linear- and phase-relation-parts
vamin = np.argmin(va[0]) 
    
# initial linear fit for better guess of eta_c
linparams = sp.optimize.curve_fit(flin,Eta[vamin:len(Eta)],va[0][vamin:len(va[0])],p0=[0.2,-0.5],sigma=errorbars[0][vamin:len(errorbars[0])])
valin = flin(Eta,*linparams[0])

# set fit parameters for square/lin fit
eta_c_ind = np.where(np.absolute(va[0] - valin) > prec)[0][-1]
x = Eta[0:eta_c_ind]
y = va[0][0:eta_c_ind]
err = errorbars[0][0:eta_c_ind]

# fit it with square + line
sqwidth = sp.optimize.curve_fit(finvsquare,x,y,p0=1/6,sigma=err)
linparams = sp.optimize.curve_fit(flin,Eta[eta_c_ind:len(Eta)],va[0][eta_c_ind:len(va[0])],p0=[0.2,-0.5],sigma=errorbars[0][eta_c_ind:len(errorbars[0])])

# calculate second approx for eta_c from fit
p,q = linparams[0][0]/sqwidth[0][0]**2,(linparams[0][1]-1)/sqwidth[0][0]**2
eta_c_0 = -p/2 + np.sqrt((p/2)**2 - q)

# recalculate fit parameters
eta_c_ind = np.where(Eta < eta_c_0)[0][-1]
x = Eta[0:eta_c_ind]
y = va[0][0:eta_c_ind]
err = errorbars[0][0:eta_c_ind]

# fit beta
xlog = np.log(eta_c_0-x)
ylog = np.log(y)
errlog = np.log(err)

# logarithmic fit
def fphase_temp_log(x, beta, offset):
    return beta * x + offset
tempparams = sp.optimize.curve_fit(fphase_temp_log,xlog,ylog,p0=[0.5,0],sigma=errlog)

# final fit for eta_c and beta using former fits as start values
def fphase(x, eta_c, beta, offset):
    return (eta_c - x)**beta * np.e**offset
def fphaselog(x, eta_c, beta, offset):
    return np.log(eta_c - x)*beta + offset
phaseparams = sp.optimize.curve_fit(fphaselog,x,ylog,p0=[eta_c_0,*tempparams[0]],sigma=errlog)

Plot $v_a$ over $\eta$



In [ ]:

    
plt.figure(figsize=figsize)
for c,n in enumerate(N):
    plt.errorbar(Eta, va[c], yerr=errorbars[c], fmt='.', label="N="+str(n))
x = np.linspace(Eta[0],Eta[-1])
xphase = x[np.where(phaseparams[0][0] > x)]
plt.plot(x,finvsquare(x,sqwidth[0]), label='Square Fit', linewidth=0.5)
plt.plot(xphase,fphase(xphase,*phaseparams[0]), label='Phase relation Fit')
plt.plot(x,flin(x,*linparams[0]), label='linear Fit', linewidth=0.5)

plt.xlabel("$\\eta$")
plt.ylabel("$v_a$")
plt.xlim([0,5.5])
plt.ylim([0,1])
plt.legend()

plt.savefig(figpath+'va_over_eta.eps')

Plot $v_a$ over $(\eta_c - \eta)/\eta_c)$



In [ ]:

    
plt.figure(figsize=figsize)
plt.ylim(ymin=0.2)
plt.xlim(xmin=0.01)
eta_c = 4.5  # This is a guess of the critical eta value. There must be a better way of determining it
for c,n in enumerate(N):
    plt.plot( (eta_c-Eta)/eta_c, va[c],'.',label="N="+str(n))
ca = plt.gca()
ca.set_xscale('log')
ca.set_yscale('log')
plt.xlabel("$\\frac{\\eta_c - \\eta}{\\eta_c}$")
plt.ylabel("$v_a$")



x = (eta_c-Eta)/eta_c
select = x > 0
x = np.log(x[select])
y = np.log(va[-1][select])
coef = np.polyfit(x,y,1)
plt.plot(np.logspace(-1,0,10), np.logspace(-1,0,10)**coef[0], label="$\\beta=$"+str(coef[0]))
plt.legend();

# if you come up with a better naming scheme PLEASE change this
plt.savefig(figpath+'va_over_etac_minus_eta_over_etac.eps')
print('eta_c='+str(eta_c))

Run with varying density



In [15]:

    
res = 15
time = 1000
averaging_time = 500
repeat = 5
eta = .3
Rho = np.logspace(-3,-0, res)
N = [100]

va = np.zeros((len(N), res))
vas = np.zeros(repeat)
errorbars = np.zeros_like(va)
variance = np.zeros(repeat)
for c,n in enumerate(N):
    for i,rho in enumerate(Rho):
        for j in range(repeat):
            f = birds.Flock(n=n, eta=eta, rho=rho)
            for t in range(time):
                f.run()
            va_avg, vari = avg_with_error(f, averaging_time)
            vas[j] = va_avg
            variance[j] = vari
        va[c][i] = vas.sum()/repeat
        errorbars[c][i] = np.sqrt(variance.sum()/repeat)



In [16]:

    
plt.figure(figsize=figsize)
for c,n in enumerate(N):
    plt.errorbar(Rho, va[c], yerr=errorbars[c], fmt='.', label="N="+str(n))

plt.xlabel("$\\rho$")
plt.ylabel("$v_a$")
plt.legend()
plt.title("Alignment dependance on density");

plt.savefig(figpath+'va_over_rho.eps')









    



/usr/lib/python3.6/site-packages/matplotlib/figure.py:1999: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "

Run with varying angle

This has to be done for low rho and eta. This should be evident from the graphs above as a high value of eta reduces the alignment going to almost zero.

A high density causes more alignment, and thus if we attempt running with higher densities, they all align anyway.



In [9]:

    
res = 20
time = 1000
averaging_time = 1000
repeat = 3
eta = 0
Angle = np.linspace(1,180,res,dtype=int)
Rho= [0.01,0.1,1] # np.logspace(-3, 0, 5)
n = 100

va = np.zeros((len(Rho), res))
vas = np.zeros(repeat)
errorbars = np.zeros_like(va)
variance = np.zeros(repeat)
for c,rho in enumerate(Rho):
    for i,angle in enumerate(Angle):
        for j in range(repeat):
            f = birds.Flock(n=n, eta=eta, rho=rho, angle=angle)
            for t in range(time):
                f.run()
            va_avg, vari = avg_with_error(f, averaging_time)
            vas[j] = va_avg
            variance[j] = vari
        va[c][i] = vas.sum()/repeat
        errorbars[c][i] = np.sqrt(variance.sum()/repeat)

plt.figure()
for c,rho in enumerate(Rho):
    plt.errorbar(Angle/360*2*np.pi, va[c], yerr=errorbars[c], fmt='.', label="$\\rho$="+str(np.round(rho,decimals=4)))
plt.xlabel("$\\theta_{cone}$")
plt.ylabel("$v_a$")
plt.legend()
plt.title("Alignment dependance on angle");

plt.savefig(figpath+'va_over_angle.eps')



In [ ]: