Fig_18.ipynb

Kernel: Python 3 (system-wide)

import numpy as np
import matplotlib.pyplot as plt

from scipy.linalg import lstsq

from scipy import signal, stats
from scipy.stats import t as t_dist

plt.style.use('../jfm.mplstyle')

labels   = np.loadtxt("../Experimental_Dataset/consolidated_amplifications.tsv",usecols=(0,),dtype=str)
raw_data = np.loadtxt("../Experimental_Dataset/consolidated_amplifications.tsv",usecols=(1,2,3,4,5),dtype=float)

length = raw_data[:,0]
amp_ups_wm = raw_data[:,1] #-1
amp_ups_os = raw_data[:,2] #-1
amp_max_wm = raw_data[:,3] #-1
amp_max_os = raw_data[:,4] #-1

mask_still = np.array([ 'g00' in label for label in labels ])
mask_flow  = ~mask_still
mask_h50   = np.array([ 'h50' in label for label in labels ])
mask_h80   = np.array([ 'h80' in label for label in labels ])

mask_h50_still = np.logical_and(mask_h50,mask_still)
mask_h50_flow  = np.logical_and(mask_h50,mask_flow )
mask_h80_still = np.logical_and(mask_h80,mask_still)
mask_h80_flow  = np.logical_and(mask_h80,mask_flow )

tinv = lambda alpha, df: abs(t_dist.ppf(alpha/2, df))

x = length[mask_still]
y = amp_max_os[mask_still]

result = stats.linregress(x,y)
m, b = result.slope, result.intercept
m_conf = tinv(0.05,len(x)-2)*result.stderr
b_conf = tinv(0.05,len(x)-2)*result.intercept_stderr

fit_amp_max = m*x + b
fit_amp_max_u = (m+m_conf)*x + b + b_conf
fit_amp_max_l = (m-m_conf)*x + b - b_conf

x = length[mask_still]
y = amp_ups_os[mask_still]

result = stats.linregress(x,y)
m, b = result.slope, result.intercept
m_conf = tinv(0.05,len(x)-2)*result.stderr
b_conf = tinv(0.05,len(x)-2)*result.intercept_stderr

fit_amp_ups = m*x + b
fit_amp_max_u = (m+m_conf)*x + b + b_conf
fit_amp_max_l = (m-m_conf)*x + b - b_conf

fig, ax = plt.subplots(figsize=(5.0,1.85))

plt.plot(length[mask_still],fit_amp_ups,'k--')

kwargs = {
    'markersize':3,
}

plt.plot(length[mask_h50_still],amp_ups_os[mask_h50_still],'k^',label="no flow",**kwargs)
plt.plot(length[mask_h80_still],amp_ups_os[mask_h80_still],'ko',label="no flow",**kwargs)

plt.plot(length[mask_h50_flow],amp_ups_os[mask_h50_flow],'r^',label="with flow",**kwargs)
plt.plot(length[mask_h80_flow],amp_ups_os[mask_h80_flow],'ro',label="with flow",**kwargs)

plt.ylim(1,1.60)
plt.xlabel(r"$\lambda/b_0$")
plt.ylabel(r"$a_c/a_i$")

plt.tight_layout()
plt.savefig('Fig_18.pdf')
plt.show()