First Example, Interaction in 2D Lattice

home > Research

Tips. Skim through this tutorial if this is your first time reading this series and come back after you gain some familiarity. This post serves as a slightly less trivial motivation to study interacting system. Tutorial 2 will start with trivial 1D non-interacting chain. It will ease you into the use of this package of codes.

The main simulation script galw2D.py  is run on Python with package kero installed. Our system is as the following: initially every individual in the grid performs action P, except for one, who is allowed to perform any of the actions P, a1 or a2 with probability 0.8, 0.1 and 0.1 respectively. At one time step, each individual chooses one of its available actions, and then we update the effect of each action taken and proceed to the next time step. If an individual X takes and action a1, then its neighbors NN(X) will gain a certain probability of performing a1 in the next iteration. Likewise for a2, as shown below in figure 1.

densen.JPG

Figure 1. Lattice evolving from time step n0 to n0+1.

We can create a lot more variations. For this system, we also add that if X performs a1, then the neighbors NN(X) will lose some tendency to perform a2. Likewise if X performs a2, NN(X) will less likely perform a1 in the next step.

Analysis

We suggest how analysis can be performed. Note that all individuals have probability (P, a1, a2) = (1,0,0) except for the particular member “5” with (P, a1, a2) = (0.8, 0.1, 0.1). When X performs a1, NN(X) will be modified such that (P, a1, a2) -> (P, a1 + 0.01, a2 – 0.01). If a1 > 1.0,  since probability cannot exceed 1.0, it will be set to 1.0, while if a2 < 0 it will be set to 0. Then the probability vector will be normalized, so that P + a1 + a2 = 1. If X performs a2, then NN(X) a2 will be increased and a1 reduced similarly.

As long as “5” performs P at the beginning, nothing happens, and the system remains the same. Let n0 be the least time step when “5” perform anything other than P, and assume that it performs a1. At n0+1, neighbors of “5”, denoted NN(“5”) has 0.01/(0.01+1) approximately 0.01 chance of performing a1. Note that if any NN(“5”) performs a1 at n0+1, then the probability of “5” performing “a2” will be zero at n0+2 and the system will propagate a1 to the entire grid in the long run.

Now we compute the probability any NN(“5”) performs a1 at n0. Let us call the even N1, N2, N3, N4, each corresponding to a neighbor, then we need

P(N1\cup N2\cup N3\cup N4) 
= P(N1)+P(N2)+... -P(N1\cap N2)- ... -P(N1\cap N2\cap N3\cap N4\cap)
\approx 0.34

This applies to the case when “5” performs a2 at n0 and any of NN(“5”) performs a2 at n0+1. For now, we can see there is a significant probability the system collapses to either all a1 or all a2 in the long run. Even for such simple cases, we cannot fully account for non-obvious transition in between. I think that the way to go is to identify critical points like this, although the trouble may come when the system is chaotic and it may be interesting to note that a system might be circular in their states.


1. Set up system

In this system, we make a 4 by 3 grid of interacting individuals. This is done by setting m=4, n=3. In this tutorial we show explicitly how the probability is assigned to the network via a variable with dictionary data type, called set_of_fields[“prob”], which is the property of network we call net2D. We can add any other field such as set_of_fields[“something else”] if we want to add more interactions next time.

The variable action_index_set is required: for each action P, a1 and a2 we give it numerical label 0, 1 and 2 for many purposes. Use 0, 1, 2,… convention for consistency.

2. Run simulation

Set the duration of simulation using variable Tset. Here we run it over 1000 time steps. In this example we only have 1 field, which is the probability field set_of_fields[“prob”] deciding what is the probability of an action each member will take in one time step. At the section # Add interaction in this region, not only does each member select a random choice of action, it will influence other members, as described before. All the effects of interaction should be put in this part of the code.

3. Results

Figure 2 shows that the system has collapsed into all a1 state, i.e. at any point in time in the future each member is likely to perform a1.

The python package kero includes standard visualization methods we can use to plot the evolution of the network in time. Visit the links for a more technical explanation of each function and the plot.

galw2D.JPG

Figure 2. (a) 2D lattice of the network. The color intensity of the lattice reflects the probability of action a1 (action index=1) is chosen at the current state. The number labels are the name of network members. (b) The actions taken by each member (color-coded) in the last 20 time steps. (c) The mean number of actions taken by 3 selected network members “1”, “5” and “6”. As shown, the system has settled to the state where all members are likely performing a1 in the long run.


Further works

We would like to further develop more user-friendly interface and study more complex system that can be applied to real situation.


Codes

The codes import peripheral.py from here.

galw2D.py

import kero.multib.nDnet as nd
import matplotlib.pyplot as plt
from peripheral import *

fig=plt.figure()

# --------------------------------------------------------
# 1. Set up system
# --------------------------------------------------------
m=4
n=3
net2D = nd.Network()
net2D.build_2D_grid(m,n)
# ----------- Create Probability Vector Field -----------
# action_set = ["a1","a2"]
# action_index_set = [0,1] # specify number representation of each action
# action_set, action_index_set=net2D.attach_a_simple_probability_vector_field()
def attach_a_probability_vector_field(one_member_to_modify):
	net2D.set_of_fields["prob"] = {}
	action_set = ["P","a1","a2"]
	action_index_set = [0,1,2]
	action_probability_set = [1,0,0] # members will perform only action P, 0% a1 and a2
	for v in net2D.member_set:
		net2D.set_of_fields["prob"][v] = [action_set,action_probability_set]
	# For selected member, identified by one_member_to_modify, it has
	#  80% chance of doing P, and 10% a1,a2 each
	net2D.set_of_fields["prob"][one_member_to_modify][1] = [0.8,0.1,0.1]
	return action_set, action_index_set

action_set, action_index_set=attach_a_probability_vector_field("5")
# print(net2D.set_of_fields["prob"])
# --------------------------------------------------------
# 2. Run Simulation
# --------------------------------------------------------
start1 = time.time()
Tset=range(1000) # Simulation length
action_history = {}
for v in net2D.member_set:
	action_history[v]=[]
for t in Tset:
	for v in net2D.member_set:
		this_action_set = net2D.set_of_fields["prob"][v][0] # eg [a1,a2]
		this_probability_set = net2D.set_of_fields["prob"][v][1] # eg [0.5,0.5]
		action,action_index = nd.choose_action_continuous_probability(this_action_set, this_probability_set, N=1)
		action_history[v].append([action_index,action]) 
	# ---------------------------------------
	# Add interaction in this region
	# ---------------------------------------
	debuff_factor=0.01 # check this out, switch between 1e-6 and 1e-2
	for v in net2D.member_set:
		a_set = action_history[v][-1]
		neighbors = net2D.member_set[v]
		# print(v, " : ", a_set," -> ", neighbors)
		if a_set[1]=="a1":
			for k in neighbors:
				temp = net2D.set_of_fields["prob"][k][1]
				temp[1]=temp[1]+0.01
				temp[2]=temp[2]-debuff_factor
				if temp[2]<0: 					
					temp[2]=0 				
				if temp[1]>1.0:
					temp[1]=1
				sigma_temp=sum(temp)
				net2D.set_of_fields["prob"][k][1] = [x/sigma_temp for x in temp]


		elif a_set[1]=="a2":	
			for k in neighbors:
				temp = net2D.set_of_fields["prob"][k][1]
				temp[1]=temp[1]-debuff_factor
				temp[2]=temp[2]+0.01
				if temp[1]<0: 					
					temp[1]=0 				
				if temp[2]>1.0:
					temp[2]=1
				sigma_temp=sum(temp)
				net2D.set_of_fields["prob"][k][1] = [x/sigma_temp for x in temp]
time_keeper("Simulation",start1)
# --------------------------------------------------------
# 3. Results
# --------------------------------------------------------

# --------------------------------------------------------
ax2 = fig.add_subplot(222,projection='3d')
# --------------------------------------------------------
start2=time.time()
processor=nd.Network_processor()
processor.add_member_and_history(net2D.member_set,Tset,action_history)
processor.map_action_index(action_set,action_index_set)

processor.plot_action_over_time(
		ax2,
		plot_interval=(len(Tset)-24,None),
		title="Action-time",
		)
time_keeper("Action-time",start2)

# --------------------------------------------------------
ax3 = fig.add_subplot(223)
# --------------------------------------------------------
start3=time.time()
mean_over_N=200
interval_mean_history=processor.compute_interval_mean(mean_over_N=mean_over_N)
# print(interval_mean_history)
member_set_to_display= ["1","5","6"] # net2D.member_set # 
action_index_set_to_display=action_index_set # [0,1] # 

processor.plot_interval_mean_over_time(
		ax3,
		interval_mean_history,
		member_set_to_display,
		action_index_set_to_display,
		plotting_interval= (100,None), # (1,None)
		plotting_interval_step_length=40,
		title="Mean over last " + str(mean_over_N) +" steps"
		)
time_keeper("Mean-time",start3)

# --------------------------------------------------------
ax = fig.add_subplot(221,projection=None)
# --------------------------------------------------------
start4=time.time()
action_index_to_display=1
time_step=None
interval_mean_at_t=processor.extract_interval_mean_at_a_time_step(interval_mean_history,time_step=time_step)
# print("interval_mean_at_t:\n",interval_mean_at_t)
CS=[interval_mean_at_t[action_index_to_display][v][0] for v in interval_mean_at_t[action_index_to_display]]
# print(CS)
net2D.visualize_3D_graph(ax,color_set=CS,
						this_figure=fig,
						no_of_points_per_edge=4,
						title="2D grid: action index="+str(action_index_to_display),
						cmap=plt.cm.hot
						)	
time_keeper("Graph",start4)
plt.show()

Simulation is done using kero version: 0.5.1.