In [2]:
%matplotlib inline
In [3]:
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import polyhedra as poly
import manifold_reflected_brownian_motion as mrbm
import bg_run_script as bgrs
import bga_4_0 as bga
import statistics as sts
import bg_join_script as bgj
bga = reload(bga)
bgj = reload(bgj)
In [4]:
from matplotlib import rc
rc('font', **{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
fig_width = 10
fig_height = 6
my_figsize = (fig_width, fig_height)
fig_dpi = 200
my_dpi = fig_dpi
In [5]:
#poly_name = "tetrahedron"
poly_name = "octahedron"
#poly_name = "icosahedron"
output_str = "final"
run_str = "test3"
In [6]:
verts, face_inds, cents = getattr(poly, poly_name)()
V, E, F, S, species, f_types, adj_list, dual = bga.get_poly(poly_name)
ints, ids, paths, shell_int, shell_paths, edges, shell_edge = bga.get_bg_ss(poly_name)
Rs = bga.generate_rotations(adj_list)
num_ints = len(ints)
In [7]:
Ss, Ts = bga.get_degeneracies(ints, edges, adj_list, Rs=Rs)
rs = bga.get_rs(ints, adj_list, Rs)
Es = bga.get_open_edges(ints, adj_list)
#Es = -bga.get_closed_edges(ints, adj_list)
E_j, E_k = bgj.transision_energies(Es, edges)
In [8]:
### STATIONARY DISTRIBUTION
x_min = 0
x_max = 2.0
#my_figsize = (6.6,5)
#my_dpi = 100
betas = np.linspace(x_min, x_max, 100)
pi = np.zeros((num_ints,len(betas)))
for k, beta in enumerate(betas):
pi[:,k] = np.exp(-beta*Es)/rs
pi[0,k] = 0.0
pi[:,k] /= sum(pi[:,k])
In [9]:
plt.figure(figsize=my_figsize, dpi=my_dpi)
for j in range(1,num_ints):
plt.plot(betas, pi[j,:].T, label=j)
plt.ylim([0.0,0.2])
plt.title("Octahedron Stationary Distribution vs. Inverse Temperature")
plt.xlabel(r"Inverse Temperature $\beta$") #$\beta$")
plt.ylabel(r"Stationary Probability $\pi_j$")
plt.legend()
plt.savefig("octahedron_pi.eps", figsize=my_figsize, dpi=my_dpi)
In [10]:
plt.figure(figsize=my_figsize, dpi=my_dpi)
for j in range(1,num_ints):
plt.plot(betas, np.log(pi[j,:].T), label=j)
plt.ylim([-8,0.0])
plt.title("Octahedron Log Stationary Distribution vs. Inverse Temperature")
plt.xlabel(r"Inverse Temperature $\beta$") #$\beta$")
plt.ylabel(r"Stationary Log Probability $\log\pi_j$")
plt.legend()
plt.savefig("octahedron_log_pi.eps", figsize=my_figsize, dpi=my_dpi)
In [11]:
int_dict = bgj.load_ints(poly_name, run_str, output_str)
edge_rate_inds = bgj.find_sim_edge_relations(int_dict,
edges,
ints,
face_inds,
V,
dual,
verts,
adj_list,
Rs)
In [12]:
epsilons = 2.0**np.linspace(np.log2(0.2), np.log2(1.0), 6)
#epsilons = np.linspace(0.2,1.0,6)
tau_beta = np.zeros((num_ints, len(betas), len(epsilons)))
tau_ratio = np.zeros((num_ints, len(betas), len(epsilons)))
for k, beta in enumerate(betas):
for j, epsilon in enumerate(epsilons):
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
taus = bga.get_taus(Q, num_ints)
tau_beta[:,k,j] = taus
tau_ratio[:,k,j] = -taus*Q[-1,-1]
In [13]:
plt.figure(figsize=my_figsize, dpi=my_dpi)
for j, epsilon in enumerate(epsilons):
plt.plot(betas, np.log(tau_beta[0,:,j]), label=round(epsilon,2))
plt.ylim([2.0,8.0])
plt.title("Octahedron Log Expected Formation Time vs. Inverse Temperature")
plt.xlabel(r"Inverse Temperature $\beta$") #$\beta$")
plt.ylabel(r"Log Formation Time $\log\tau$")
plt.legend(title=r"$\epsilon$", loc=4)
plt.savefig("octahedron_tau.eps", figsize=my_figsize, dpi=my_dpi)
#plt.plot(betas, np.log(tau_beta[0,:,:]))
In [14]:
beta = 0.5
epsilon = 0.5
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
P = bga.get_dist(Q, 2.0, 100)
In [15]:
fig = plt.figure(figsize=my_figsize, dpi=my_dpi)
for k in range(len(ints)):
plt.plot(np.linspace(0.0,2.0,100), P[:,k], label=str(k+1))
plt.ylim(0.0, 0.4)
plt.title(r"Octahedron Finite Time Distributions: $\beta = $"+str(beta)+r" $ \epsilon = $" + str(epsilon) )
plt.xlabel(r"Time") #$\beta$")
plt.ylabel(r"Occupation Probability")
#plt.legend(title=r"Intermediate", loc=4)
plt.legend(loc=1)
plt.savefig("octahedron_finite_dist.eps", figsize=my_figsize, dpi=my_dpi)
In [16]:
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
In [17]:
rjs = np.array([float(rs[e[0]]) for e in edges])
bi_rates = np.array(Ss)/rjs*np.exp(-beta*np.array(E_jk))
print bi_rates
print Ss
print rjs
print E_jk
In [18]:
print Ss, rjs
print Ss/rjs
In [19]:
pos_arr = np.array([[ 1., 0.5],
[ 2., 0.5],
[ 3., 0.5],
[ 4., 0.8],
[ 4., 0.6],
[ 4., 0.4],
[ 4., 0.2],
[ 5., 2.0/3.0],
[ 5., 1.0/3.0],
[ 6., 0.75],
[ 6., 0.5],
[ 6., 0.25],
[ 7., 0.5],
[ 8., 0.5]])
#node_scale = 0.1
#pos_arr[:,1] -= 0.5
#pos_arr[:,1] *= node_scale
#pos_arr[:,1] += 0.5
pos = {k+1: pos_arr[k,:] for k in range(len(ints)-1)}
In [20]:
nx_figsize = (12,5)
nx_dpi = 200
In [21]:
# extract nodes from graph
nodes = set([n1 for n1, n2 in edges[1:]] + [n2 for n1, n2 in edges[1:]])
# create networkx graph
#G=nx.DiGraph()
G=nx.Graph()
# add nodes
for node in nodes:
G.add_node(node)
# add edges
for edge in edges[1:]:
G.add_edge(edge[0], edge[1])
# draw graph
#pos = nx.shell_layout(G)
#nx.draw(G, pos)
In [22]:
node_size_scale = 10000.0
edge_thickness_scale = 100.0
node_color = 'w'
#node_size = node_size_scale*pi[1:,10]
#edge_thickness = 1.0
#edge_thickness = edge_thickness_scale*bi_rates[1:]
edge_alpha = 0.5
edge_color = 'blue'
node_text_size = 12
In [23]:
epsilon = 0.5
beta = 0.8
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
rjs = np.array([float(rs[e[0]]) for e in edges])
bi_rates = np.array(Ss)/rjs*np.exp(-beta*np.array(E_jk))
In [24]:
e_labels = {}
for k, edge in enumerate(edges):
if 0 in edge:
continue
e_labels[(edge[0],edge[1])] = str(bi_rates[k])[:6]
print k, edge
print e_labels
In [25]:
fig = plt.figure(figsize=nx_figsize, dpi=nx_dpi)
a = nx.draw_networkx_nodes(G,
pos,
node_color=node_color)
b = nx.draw_networkx_edges(G,
pos,
alpha=edge_alpha,
edge_color=edge_color)
c = nx.draw_networkx_labels(G,
pos,
font_size=node_text_size)
d = nx.draw_networkx_edge_labels(G,
pos,
edge_labels=e_labels,
label_pos=0.65)#,
#label_pos=0.5,
#font_size=10,
#font_color='k',
#font_family='sans-serif',
#font_weight='normal',
#alpha=1.0,
#bbox=None,
#ax=None,
#rotate=True)
plt.xticks(range(1,len(ints[0])+1), range(1,len(ints[0])+1))
plt.yticks([],[])
plt.tick_params(
axis='x',
which='both',
bottom='off',
top='off',
labelbottom='on')
plt.xlabel("Number of Faces")
plt.title("Octahedron Combinatorial Configuration Space")
plt.savefig("octahedron_ccs.eps", figsize=my_figsize, dpi=my_dpi)
#plt.show()
In [26]:
pi[1:,10]
Out[26]:
In [27]:
##
#### PI AND Q
##
betas = [0.7, 0.8, 0.9]
epsilons = [0.2, 0.5, 1.0]
fig, axarr = plt.subplots(3, 3, figsize=(30,18), sharex=True, sharey=True)
plt.suptitle("Octahedron Stationary Distributions and Transitions Rates", fontsize=24)
for k, beta in enumerate(betas):
#pi = np.zeros((num_ints)
pi = np.exp(-beta*Es)/rs
pi[0] = 0.0
pi /= sum(pi)
for j, epsilon in enumerate(epsilons):
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
rjs = np.array([float(rs[e[0]]) for e in edges])
bi_rates = np.array(Ss)/rjs*np.exp(-beta*np.array(E_jk))
node_size_scale = 10000.0
edge_thickness_scale = 400.0
node_color = 'w'
node_size = node_size_scale*pi[1:]
#edge_thickness = 1.0
edge_thickness = edge_thickness_scale*bi_rates[1:]
edge_alpha = 0.5
edge_color = 'blue'
node_text_size = 12
#fig = plt.figure(figsize=nx_figsize, dpi=nx_dpi)
axarr[j,k].plot([],[])
a = nx.draw_networkx_nodes(G,
pos,
node_size=node_size,
node_color=node_color,
ax=axarr[j,k])
b = nx.draw_networkx_edges(G,
pos,
width=edge_thickness,
alpha=edge_alpha,
edge_color=edge_color,
ax=axarr[j,k])
#c = nx.draw_networkx_labels(G,
# pos,
# font_size=node_text_size)
plt.xticks(range(1,len(ints[0])+1), range(1,len(ints[0])+1))
plt.yticks([],[])
plt.xticks([],[])
plt.tick_params(
axis='x',
which='both',
bottom='off',
top='off',
labelbottom='off')
#axarr[1,1]
for k, beta in enumerate(betas):
axarr[2,k].set_xlabel(r"$\beta =$ " + str(beta), fontsize=18)
for j, epsilon in enumerate(epsilons):
axarr[j,0].set_ylabel(r"$\epsilon =$ " + str(epsilon), fontsize=18)
fig.subplots_adjust(hspace=0)
fig.subplots_adjust(wspace=0)
plt.savefig("octahedron_pi_Q_grid.eps", figsize=my_figsize, dpi=my_dpi)
In [28]:
##
#### PATH OCCUPATION PROBABILITY
##
betas = [0.7, 0.8, 0.9]
epsilons = [0.2, 0.5, 1.0]
fig, axarr = plt.subplots(3, 3, figsize=(30,18), sharex=True, sharey=True)
plt.suptitle("Octahedron Stationary Distributions and Transitions Rates", fontsize=24)
for k, beta in enumerate(betas):
#pi = np.zeros((num_ints)
pi = np.exp(-beta*Es)/rs
pi[0] = 0.0
pi /= sum(pi)
for j, epsilon in enumerate(epsilons):
forward_rates, backward_rates, E_jk = bgj.get_rates(epsilon,
beta,
int_dict,
edges,
edge_rate_inds,
E_j,
E_k,
Ss,
Ts)
Q = bgj.rates_to_Q(forward_rates, backward_rates, edges, num_ints)
rjs = np.array([float(rs[e[0]]) for e in edges])
bi_rates = np.array(Ss)/rjs*np.exp(-beta*np.array(E_jk))
node_size_scale = 10000.0
edge_thickness_scale = 400.0
node_color = 'w'
node_size = node_size_scale*pi[1:]
#edge_thickness = 1.0
edge_thickness = edge_thickness_scale*bi_rates[1:]
edge_alpha = 0.5
edge_color = 'blue'
node_text_size = 12
#fig = plt.figure(figsize=nx_figsize, dpi=nx_dpi)
axarr[j,k].plot([],[])
a = nx.draw_networkx_nodes(G,
pos,
node_size=node_size,
node_color=node_color,
ax=axarr[j,k])
b = nx.draw_networkx_edges(G,
pos,
width=edge_thickness,
alpha=edge_alpha,
edge_color=edge_color,
ax=axarr[j,k])
#c = nx.draw_networkx_labels(G,
# pos,
# font_size=node_text_size)
plt.xticks(range(1,len(ints[0])+1), range(1,len(ints[0])+1))
plt.yticks([],[])
plt.xticks([],[])
plt.tick_params(
axis='x',
which='both',
bottom='off',
top='off',
labelbottom='off')
#axarr[1,1]
for k, beta in enumerate(betas):
axarr[2,k].set_xlabel(r"$\beta =$ " + str(beta), fontsize=18)
for j, epsilon in enumerate(epsilons):
axarr[j,0].set_ylabel(r"$\epsilon =$ " + str(epsilon), fontsize=18)
fig.subplots_adjust(hspace=0)
fig.subplots_adjust(wspace=0)
#plt.savefig("octahedron_pi_Q_grid.eps", figsize=my_figsize, dpi=my_dpi)
In [28]: