In [135]:
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
def plot_3d_wireframe(filename):
data = np.load(filename)
beta = data['beta']
num_of_strings = data['num_of_strings']
L = data['L']
frames = data['frames']
distance_list = data['distance_list']
path_length = data['path_length']
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
hist, xedges, yedges = np.histogram2d(distance_list, path_length, bins=25)
xpos, ypos = np.meshgrid(xedges[:-1] + (xedges[1] - xedges[0]) / 2.,
yedges[:-1] + (yedges[1] - yedges[0]) / 2.)
zpos = hist.T
ax.plot_wireframe(xpos, ypos, zpos, rstride=1)
ax.plot(xpos[0], xpos[0], lw=2)
ax.set_xlim(xedges[0], xedges[-1])
ax.set_ylim(yedges[0], yedges[-1])
ax.set_xlabel('Distance')
ax.set_ylabel('Path length')
ax.set_title('Path length and distances between two points in the cluster'
+ r'($\beta = %2.2f$)' % beta)
plt.show()
In [12]:
result_data_path = {
'0': "./results/data/distances/beta=0.00_161012_171430.npz",
'5': "./results/data/distances/beta=5.00_161012_171649.npz",
'10': "./results/data/distances/beta=10.00_161012_172119.npz",
'15': "./results/data/distances/beta=15.00_161012_172209.npz",
'20': "./results/data/distances/beta=20.00_161012_172338.npz",
'0-2': "./results/data/distances/beta=0.00_161015_153311.npz",
'5-2': "./results/data/distances/beta=5.00_161015_153838.npz",
'10-2': "./results/data/distances/beta=10.00_161015_154048.npz",
'15-2': "./results/data/distances/beta=15.00_161015_154136.npz",
'20-2': "./results/data/distances/beta=20.00_161015_154419.npz"
}
In [2]:
plot_3d_wireframe(result_data_path['0'])
In [3]:
plot_3d_wireframe(result_data_path['5'])
In [4]:
plot_3d_wireframe(result_data_path['10'])
In [5]:
plot_3d_wireframe(result_data_path['15'])
In [6]:
plot_3d_wireframe(result_data_path['20'])
In [4]:
plot_3d_wireframe(result_data_path['0-2'])
In [7]:
plot_3d_wireframe(result_data_path['5-2'])
In [8]:
plot_3d_wireframe(result_data_path['10-2'])
In [9]:
plot_3d_wireframe(result_data_path['15-2'])
In [10]:
plot_3d_wireframe(result_data_path['20-2'])
In [67]:
def plot_tortuosity(filename):
data = np.load(filename)
beta = data['beta']
num_of_strings = data['num_of_strings']
L = data['L']
frames = data['frames']
distance_list = data['distance_list']
path_length = data['path_length']
num_of_pairs = 300
x = np.array(path_length)
x = x.reshape((num_of_strings, int(x.shape[0] / num_of_strings / num_of_pairs), num_of_pairs))
y = np.array(distance_list)
y = y.reshape((num_of_strings, int(y.shape[0] / num_of_strings/ num_of_pairs), num_of_pairs))
y = y / x
y_ave = np.average(y, axis=2)
L = x[0, :, 0]
print L
lambda_ave = np.average(y_ave, axis=0)
print lambda_ave
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(L, lambda_ave)
ax.set_xlabel(r'Path length $L$')
ax.set_ylabel(r'Tortuosity $T$')
ax.set_title('Tortuosity $T = \lambda_{\text{avg}} / L$ '
+ r'($\beta = %2.2f$)' % beta)
plt.show()
In [68]:
plot_tortuosity(result_data_path['0-2'])
In [69]:
def calc_tortuosity_for_each_beta(filename):
data = np.load(filename)
beta = data['beta']
num_of_strings = data['num_of_strings']
L = data['L']
frames = data['frames']
distance_list = data['distance_list']
path_length = data['path_length']
num_of_pairs = 300
x = np.array(path_length)
x = x.reshape((num_of_strings, int(x.shape[0] / num_of_strings / num_of_pairs), num_of_pairs))
y = np.array(distance_list)
y = y.reshape((num_of_strings, int(y.shape[0] / num_of_strings/ num_of_pairs), num_of_pairs))
y = y / x
y_ave = np.average(y, axis=2)
L = x[0, :, 0]
lambda_ave = np.average(y_ave, axis=0)
return L, lambda_ave
In [85]:
fig, ax = plt.subplots()
for i in [0, 5, 10, 15, 20]:
fn = result_data_path['%d-2' % i]
ax.plot(*calc_tortuosity_for_each_beta(fn), ls='', marker='.', label=r'$\beta = %2.2f$' % float(i))
ax.set_xlabel(r'Path length $L$')
ax.set_ylabel(r'Tortuosity $T$')
ax.set_title('Tortuosity $T = \lambda_{\mathrm{avg}} / L$')
ax.legend(loc='best')
plt.show()
In [86]:
fig, ax = plt.subplots()
for i in [0, 5, 10, 15, 20]:
fn = result_data_path['%d-2' % i]
ax.loglog(*calc_tortuosity_for_each_beta(fn), ls='', marker='.', label=r'$\beta = %2.2f$' % float(i))
ax.set_xlabel(r'Path length $L$')
ax.set_ylabel(r'Tortuosity $T$')
ax.set_title('Tortuosity $T = \lambda_{\mathrm{avg}} / L$')
ax.legend(loc='best')
plt.show()
collapseできないか?
In [88]:
from optimize import Optimize_powerlaw
In [89]:
def fit_to_powerlaw_for_beta(i):
fn = result_data_path['%d-2' % i]
L, T = calc_tortuosity_for_each_beta(fn)
optimizer = Optimize_powerlaw(args=(L, T), parameters=[0., -0.5])
result = optimizer.fitting()
fig, ax = plt.subplots()
ax.loglog(L, T, ls='', marker='.', label=r'$\beta = %2.2f$' % float(i))
ax.loglog(L, optimizer.fitted(L), lw=2, label=r'$\alpha = %f$' % result['D'])
ax.set_xlabel(r'Path length $L$')
ax.set_ylabel(r'Tortuosity $T$')
ax.set_title('Tortuosity $T = \lambda_{\mathrm{avg}} / L$')
ax.legend(loc='best')
plt.show()
In [90]:
fit_to_powerlaw_for_beta(0)
In [91]:
fit_to_powerlaw_for_beta(5)
In [92]:
fit_to_powerlaw_for_beta(10)
In [93]:
fit_to_powerlaw_for_beta(15)
In [94]:
fit_to_powerlaw_for_beta(20)
In [117]:
from ipywidgets import interact, interactive, fixed
from IPython.html import widgets
In [136]:
Ls, Ts = {}, {}
for i in [0, 5, 10, 15, 20]:
fn = result_data_path['%d-2' % i]
L, T = calc_tortuosity_for_each_beta(fn)
Ls[i] = L
Ts[i] = T
In [144]:
@interact(c=(-15., 5., 0.01), d=(-15., 5., 0.01))
def collapsed(c=-10., d=-8.):
fig, ax = plt.subplots()
for i in [5, 10, 15, 20]:
L, T = Ls[i], Ts[i]
ax.plot(float(i) / (L**c), T * (L**(-d)) , ls='', marker='.', label=r'$\beta = %2.2f$' % float(i))
ax.set_xlabel(r'$\beta / L^{c}$')
ax.set_ylabel(r'$T L^{-d}$')
ax.set_title(r'scaling ($c = %2.2f, d = %2.2f$)' % (c, d))
ax.legend(loc='best')