Training an allosteric network

Begin with a few package imports:

from IPython.display import HTML

import numpy as np
import matplotlib.pyplot as plt

from packing_utils import *
from allosteric_utils import *
from plot_imports import *

Network generation

Generate a contact network from a disk packing simulation to use as the elastic network:

n = 64 # number of disks/nodes
seed = 0

net = Packing(n, rfac=0.8, seed=seed)
net.params['contact'] = 0.02 # reduce the default contact repulsion
net.generate()

progress: 100%|#####################################################################################################################| 1000.00/1000.00 [00:03<00:00]

Use the network graph to initialize an Allosteric Class object.

%matplotlib notebook

allo = Allosteric(net.graph, dim=2)
allo.plot()

Automatically select one source and one target:

allo.add_targets(seed=12)
allo.add_sources(seed=12)

Save the untrained network to a file:

allo.save('allo_{:d}.txt'.format(n))

Training

Specify a source and target strain to apply, training time, and training method.

es = 0.2                # source strain of 20%
et = 0.2                # target strain of 20%
ka = 100.               # stiffness of spring for applied strain
duration = 4e7          # training time
frames = 200            # number of output frames
train = 2               # 1 trains bond rest lengths, 2 trains bond stiffnesses
method = 'learning'     # learning rule to use, 'aging' (directed aging) or 'learning' (coupled learning)

# compute the mean square bond length as a normalization factor
l2 = np.mean([edge[2]['length']**2 for edge in allo.graph.edges(data=True)])

eta = 1e-1              # nudge factor
alpha = 1e-3/l2         # learning rate
vmin = 1e-3             # minimum allowed stiffness
vsmooth = 1e-2          # stiffness value at which to begin smooth ramp down to vmin

Initially, all bond stiffnesses are equal to 1:

%matplotlib notebook

k = allo.distribution_plot(kind='stiffness', vmin=0, vmax=2)

Since the network is untrained, applying a strain at the source (blue curve) produces no strain at the target (red curve):

%matplotlib notebook

allo.reset_init() # set to equilibrium
sol = allo.solve(duration=duration, frames=frames, T=0, applied_args=(es, 0, ka))
allo.strain_plot()

progress: 100%|############################################################################################################9| 40000000.00/40000000.00 [00:03<00:00]

Train the system:

sol = allo.solve(duration=duration, frames=frames, T=0, applied_args=(es, et, ka),
           train=train, method=method, eta=eta, alpha=alpha, vmin=vmin, vsmooth=vsmooth)

progress: 100%|#############################################################################################################| 40000000.00/40000000.00 [00:05<00:00]

We can view the increase in the target strain over time:

%matplotlib notebook

allo.strain_plot()

Visualization

View the final distribution of stiffnesses:

%matplotlib notebook

k = allo.distribution_plot(kind='stiffness', vmin=0, vmax=2, nbins=25)

Also visualize which bonds weakened (reds) and which bonds strengthened (blues):

%matplotlib notebook

allo.reset_init() # reset network to its equilibrium state for visualization

cmap = continuous_cmap([pal['red'], pal['yellow'], np.array([0.8,0.8,0.8]),
                        pal['blue'], pal['purple']], [0,0.1,0.5,0.9,1])
allo.color_plot(cmap, vmin=0, vmax=2)

Animation

Apply a sinusoidal strain at the source and monitor the target:

%matplotlib notebook

allo.reset_init()
period = 2e7

allo.solve(duration=duration, frames=frames, T=period, applied_args=(es, 0, ka))
allo.strain_plot()

progress: 100%|#############################################################################################################| 40000000.00/40000000.00 [00:00<00:00]

The following cell produces a movie of the network as it undergoes the sinusoidal strain:

%matplotlib inline

ani = allo.animate()
HTML(ani.to_html5_video())

Save the trained network:

allo.save('allo_{:d}_trained.txt'.format(n))