In Part 2 we introduced the 'language' in which a system must be described in order for it to be interpretable by cadCAD and some of the basic concepts of the library:
In the previous example, we observed how two robotic arms acting in parallel could result in counterintuitive system level behavior despite the simplicity of the individual robotic arm policies. In this notebook we'll introduce the concept of networks. This done by extending from two boxes of marbles to n boxes which are the nodes in our network. Furthermore, there are are going to be arms between some of the boxes but not others forming a network where the arms are the edges.
The robot and the marbles
balls
) with an integer number of marbles in each; a network of robotic arms is capable of taking a marble from their one of their boxes and dropping it into the other one.from cadCAD.engine import ExecutionMode, ExecutionContext, Executor
from cadCAD.configuration import Configuration
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
T = 50 #iterations in our simulation
n=7 #number of boxes in our network
m= 2 #for barabasi graph type number of edges is (n-2)*m
g = nx.barabasi_albert_graph(n, m)
G = g.to_directed()
for ed in g.edges:
G.remove_edge(ed[1],ed[0])
k = len(G.edges)
def local_energy_function(node, G, balls):
src_balls = balls[node]
#barrier force field
v = -np.log(1+src_balls)
#laplacian force field
nbrs = [e[1] for e in G.edges if e[0]==node]
for nbr in nbrs:
dst_balls = balls[nbr]
v= v+(src_balls - dst_balls)**2
return v
balls = np.zeros(n,)
balance = np.zeros(n,)
energy = np.zeros(n,)
for node in G.nodes:
rv = np.random.randint(1,25)
G.nodes[node]['initial_balls'] = rv
balls[node] = rv
for node in G.nodes:
energy[node] = local_energy_function(node, G, balls)
scale=100
nx.draw_kamada_kawai(G, node_size=balls*scale,labels=nx.get_node_attributes(G,'initial_balls'))
initial_conditions = {'balls':balls, 'balance': balance, 'energy':energy}
initial_conditions
{'balance': array([0., 0., 0., 0., 0., 0., 0.]), 'balls': array([ 6., 3., 12., 6., 9., 17., 4.]), 'energy': array([159.05408985, 79.61370564, 42.43505064, 128.05408985, -2.30258509, 166.10962824, -1.60943791])}
#input the deltas along the edges and update the boxes
#mechanism: edge by node dimensional operator
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# We make the state update functions less "intelligent",
# ie. they simply add the number of marbles specified in _input
# (which, per the policy function definition, may be negative)
def conserved_flow(delta_balls, s):
new_balls = s['balls']
for e in G.edges:
move_ball = delta_balls[e]
src = e[0]
dst = e[1]
if (new_balls[src] >= move_ball) and (new_balls[dst] >= -move_ball):
new_balls[src] = new_balls[src]-move_ball
new_balls[dst] = new_balls[dst]+move_ball
return new_balls
def update_balls(params, step, sL, s, _input):
delta_balls = _input['delta']
new_balls = conserved_flow(delta_balls, s)
key = 'balls'
value = new_balls
return (key, value)
def update_balance(params, step, sL, s, _input):
delta_balls = _input['delta']
new_balls = conserved_flow(delta_balls, s)
for node in G.nodes:
new_local_energy = local_energy_function(node, G, new_balls)
balance[node] = balance[node]-(new_local_energy-s['energy'][node])
key = 'balance'
value = balance
return (key, value)
def update_energy(params, step, sL, s, _input):
delta_balls = _input['delta']
new_balls = conserved_flow(delta_balls, s)
for node in G.nodes:
energy[node] = local_energy_function(node, G, new_balls)
key = 'energy'
value = energy
return (key, value)
# this time lets make three kinds of robots
def greedy_robot(src_balls, dst_balls):
#robot wishes to accumlate balls at its source
#takes half of its neighbors balls
if src_balls < dst_balls:
delta = -np.floor(dst_balls/2)
else:
delta = 0
return delta
def fair_robot(src_balls, dst_balls):
#robot follows the simple balancing rule
delta = np.sign(src_balls-dst_balls)
return delta
def giving_robot(src_balls, dst_balls):
#robot wishes to gice away balls one at a time
if src_balls > 1:
delta = 1
else:
delta = 0
return delta
#in the previous version the robots were assigned to the edges
#moving towards an agent based model formulation we assign the stratgies
#instead to the nodes
robot_strategies = [greedy_robot,fair_robot, giving_robot]
for node in G.nodes:
nstrats = len(robot_strategies)
rv = np.random.randint(0,nstrats)
G.nodes[node]['strat'] = robot_strategies[rv]
for e in G.edges:
owner_node = e[0]
G.edges[e]['strat'] = G.nodes[owner_node]['strat']
#Policy: node by edge dimensional operator
#input the states of the boxes output the deltas along the edges
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# We specify the robotic networks logic in a Policy Function
# unlike previous examples our policy controls a vector valued action, defined over the edges of our network
def robotic_network(params, step, sL, s):
delta_balls = {}
for e in G.edges:
src = e[0]
src_balls = s['balls'][src]
dst = e[1]
dst_balls = s['balls'][dst]
#transfer balls according to specific robot strat
srat = G.edges[e]['strat']
delta_balls[e] = srat(src_balls,dst_balls)
return({'delta': delta_balls})
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# In the Partial State Update Blocks, the user specifies if state update functions will be run in series or in parallel
partial_state_update_blocks = [
{
'policies': { # The following policy functions will be evaluated and their returns will be passed to the state update functions
'robotic_network': robotic_network
},
'variables': { # The following state variables will be updated simultaneously
'balls': update_balls,
'energy': update_energy,
'balances': update_balance,
}
}
]
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Settings of general simulation parameters, unrelated to the system itself
# `T` is a range with the number of discrete units of time the simulation will run for;
# `N` is the number of times the simulation will be run (Monte Carlo runs)
# In this example, we'll run the simulation once (N=1) and its duration will be of 10 timesteps
# We'll cover the `M` key in a future article. For now, let's leave it empty
simulation_parameters = {
'T': range(T),
'N': 1,
'M': {}
}
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# The configurations above are then packaged into a `Configuration` object
config = Configuration(initial_state=initial_conditions, #dict containing variable names and initial values
partial_state_update_blocks=partial_state_update_blocks, #dict containing state update functions
sim_config=simulation_parameters #dict containing simulation parameters
)
exec_mode = ExecutionMode()
exec_context = ExecutionContext(exec_mode.single_proc)
executor = Executor(exec_context, [config]) # Pass the configuration object inside an array
raw_result, tensor = executor.main() # The `main()` method returns a tuple; its first elements contains the raw results
df = pd.DataFrame(raw_result)
single_proc: [<cadCAD.configuration.Configuration object at 0x150fa7e438>]
balls_list = [b for b in df.balls]
energy_list = [b for b in df.energy]
balance_list = [b for b in df.balance]
total_energy = [sum(b) for b in df.energy]
total_balance = [sum(b) for b in df.balance]
plt.plot(balance_list, energy_list, 'o-', alpha=.4)
[<matplotlib.lines.Line2D at 0x151898def0>, <matplotlib.lines.Line2D at 0x151899a160>, <matplotlib.lines.Line2D at 0x151899a2b0>, <matplotlib.lines.Line2D at 0x151899a400>, <matplotlib.lines.Line2D at 0x151899a550>, <matplotlib.lines.Line2D at 0x151899a6a0>, <matplotlib.lines.Line2D at 0x151899a7f0>]
plt.semilogy(balance_list, balls_list, 'o-', alpha=.4)
plt.legend(['Box #'+str(node)+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]" for node in range(n)], ncol = 5,loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x1518a22358>
plt.loglog(balls_list, energy_list, 'o-', alpha=.4)
[<matplotlib.lines.Line2D at 0x1518bcbb70>, <matplotlib.lines.Line2D at 0x1518c371d0>, <matplotlib.lines.Line2D at 0x1518c37320>, <matplotlib.lines.Line2D at 0x1518c37470>, <matplotlib.lines.Line2D at 0x1518c375c0>, <matplotlib.lines.Line2D at 0x1518c37710>, <matplotlib.lines.Line2D at 0x1518c37860>]
plt.plot(df.timestep.values, balls_list)
plt.xlabel('iteration')
plt.ylabel('Number of Marbles')
plt.title('Marbles in each box')
plt.legend(['Box #'+str(node)+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]" for node in range(n)], ncol = 5,loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x15187bb828>
plt.semilogy(df.timestep.values, energy_list)
plt.xlabel('iteration')
plt.ylabel('Local Engergy')
plt.title('Energy Levels')
plt.legend(['Box #'+str(node)+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]" for node in range(n)], ncol = 5,loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x1518ed8860>
plt.plot(df.timestep.values, balance_list)
plt.xlabel('iteration')
plt.ylabel('Net funds for each Robot')
plt.title('Net funds for each Robot')
plt.legend(['Box #'+str(node)+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]" for node in range(n)], ncol = 5,loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x1519259470>
end_state_balls = np.array([b for b in balls_list[-1]])
#avg_balls = np.array([np.mean(b) for b in balls_list])
for node in G.nodes:
G.nodes[node]['final_balls'] = end_state_balls[node]
#G.nodes[node]['avg_balls'] = avg_balls[node]
cmap = plt.cm.jet
Nc = cmap.N
Ns = len(robot_strategies)
d = int(Nc/Ns)
k = len(G.edges)
strat_color = []
for e in G.edges:
for i in range(Ns):
if G.edges[e]['strat']==robot_strategies[i]:
color = cmap(i*d)
G.edges[e]['color'] = color
strat_color = strat_color+[color]
nx.draw_kamada_kawai(G, node_size=end_state_balls*scale, labels=nx.get_node_attributes(G,'final_balls'), edge_color=strat_color)
for node in G.nodes:
print("box #"+str(node)+" has "+str(end_state_balls[node])+" marbles: "+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]")
box #0 has 28.0 marbles: [greedy_robot] box #1 has 9.0 marbles: [greedy_robot] box #2 has 0.0 marbles: [fair_robot] box #3 has 2.0 marbles: [giving_robot] box #4 has 3.0 marbles: [giving_robot] box #5 has 1.0 marbles: [giving_robot] box #6 has 14.0 marbles: [greedy_robot]
rolling_avg_balls = np.zeros((T+1, n))
for t in range(T+1):
for node in G.nodes:
for tau in range(t):
rolling_avg_balls[t,node] = (tau)/(tau+1)*rolling_avg_balls[t, node]+ 1/(tau+1)*balls_list[tau][node]
plt.plot(range(len(rolling_avg_balls)),rolling_avg_balls)
plt.xlabel('iteration')
plt.ylabel('number of balls')
plt.title('time average balls in each box')
plt.legend(['Box #'+str(node)+"["+str(G.nodes[node]['strat']).split(" ")[1]+"]" for node in range(n)], ncol = 5,loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x15195b50b8>
avg_balls = np.zeros(n)
for node in G.nodes:
#store in graph
G.nodes[node]['avg_balls'] = int(10*(rolling_avg_balls[-1][node]))/10
#store as vector
avg_balls[node] = G.nodes[node]['avg_balls']
#need both for plotting
nx.draw_kamada_kawai(G, node_size=avg_balls*scale, labels=nx.get_node_attributes(G,'avg_balls'))
print(rolling_avg_balls[-1])
[27.56 8.88 2.3 1.12 3.24 1.3 12.6 ]
cmap = plt.cm.jet
Nc = cmap.N
Nt = len(simulation_parameters['T'])
dN = int(Nc/Nt)
cmaplist = [cmap(i*dN) for i in range(Nt)]
for t in simulation_parameters['T']:
state = np.array([b for b in balls_list[t]])
nx.draw_kamada_kawai(G, node_size=state*scale, alpha = .4/(t+1), node_color = cmaplist[t])
nx.draw_kamada_kawai(G, node_size=avg_balls*scale, labels=nx.get_node_attributes(G,'avg_balls'), edge_color=strat_color)
print(end_state_balls)
[28. 9. 0. 2. 3. 1. 14.]
avg_balls
array([27.5, 8.8, 2.3, 1.1, 3.2, 1.2, 12.6])