fermi-function update rule for heterogenouse

coevolve.py

@@ -1,8 +1,9 @@
 # -*- coding: utf-8 -*-
+# -*- Author: shaodan -*-
+# -*-  2015.07.11 -*-
 
 import numpy as np
 import networkx as nx
-import matplotlib.pyplot as plt
 
 
 class CoEvolveRule(object):
@@ -13,7 +14,7 @@ def rewire(self, graph, s_e, anchor):
         pass
 
 
-class Rewire(CoEvolveRule):
+class Preference(CoEvolveRule):
 
     def __init__(self, order=4):
         super(self.__class__, self).__init__(order)
@@ -44,7 +45,7 @@ def rewire(self, graph, s_e, anchor):
             graph.add_edge(anchor, new)
         return old, new
 
-    def rewire_new(self, G, s_e, anchor):
+    def rewire_one(self, G, s_e, anchor):
         # rewire only one link
         if s_e == 0:
             p = np.ones(self.N)

draw.py

@@ -2,30 +2,59 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import networkx as nx
 from matplotlib import cm
 from numpy.random import randn
 
-# Make plot with vertical (default) colorbar
-fig, ax = plt.subplots()
+def degree_histogram():
+    G = nx.gnp_random_graph(100,0.02)
 
-data = np.clip(randn(250, 250), -1, 1)
+    degree_sequence=sorted(nx.degree(G).values(),reverse=True) # degree sequence
+    #print "Degree sequence", degree_sequence
+    dmax=max(degree_sequence)
 
-cax = ax.imshow(data, interpolation='nearest', cmap=cm.coolwarm)
-ax.set_title('Gaussian noise with vertical colorbar')
+    plt.loglog(degree_sequence,'b-',marker='o')
+    plt.title("Degree rank plot")
+    plt.ylabel("degree")
+    plt.xlabel("rank")
 
-# Add colorbar, make sure to specify tick locations to match desired ticklabels
-cbar = fig.colorbar(cax, ticks=[-1, 0, 1])
-cbar.ax.set_yticklabels(['< -1', '0', '> 1'])# vertically oriented colorbar
+    # draw graph in inset
+    plt.axes([0.45,0.45,0.45,0.45])
+    Gcc=sorted(nx.connected_component_subgraphs(G), key = len, reverse=True)[0]
+    pos=nx.spring_layout(Gcc)
+    plt.axis('off')
+    nx.draw_networkx_nodes(Gcc,pos,node_size=20)
+    nx.draw_networkx_edges(Gcc,pos,alpha=0.4)
 
-# Make plot with horizontal colorbar
-fig, ax = plt.subplots()
+    # plt.savefig("degree_histogram.png")
+    plt.show()
 
-data = np.clip(randn(250, 250), -1, 1)
 
-cax = ax.imshow(data, interpolation='nearest', cmap=cm.afmhot)
-ax.set_title('Gaussian noise with horizontal colorbar')
+def unknown():
+    # Make plot with vertical (default) colorbar
+    fig, ax = plt.subplots()
 
-cbar = fig.colorbar(cax, ticks=[-1, 0, 1], orientation='horizontal')
-cbar.ax.set_xticklabels(['Low', 'Medium', 'High'])# horizontal colorbar
+    data = np.clip(randn(250, 250), -1, 1)
 
-plt.show()
+    cax = ax.imshow(data, interpolation='nearest', cmap=cm.coolwarm)
+    ax.set_title('Gaussian noise with vertical colorbar')
+
+    # Add colorbar, make sure to specify tick locations to match desired ticklabels
+    cbar = fig.colorbar(cax, ticks=[-1, 0, 1])
+    cbar.ax.set_yticklabels(['< -1', '0', '> 1'])# vertically oriented colorbar
+
+    # Make plot with horizontal colorbar
+    fig, ax = plt.subplots()
+
+    data = np.clip(randn(250, 250), -1, 1)
+
+    cax = ax.imshow(data, interpolation='nearest', cmap=cm.afmhot)
+    ax.set_title('Gaussian noise with horizontal colorbar')
+
+    cbar = fig.colorbar(cax, ticks=[-1, 0, 1], orientation='horizontal')
+    cbar.ax.set_xticklabels(['Low', 'Medium', 'High'])# horizontal colorbar
+
+    plt.show()
+
+
+degree_histogram()

egt.py

@@ -17,12 +17,12 @@
 # G = nx.random_regular_graph(4, 1000)
 # G = nx.convert_node_labels_to_integers(nx.grid_2d_graph(10,10))
 # G = nx.star_graph(10)
-G = nx.random_graphs.watts_strogatz_graph(1000, 4, 0.3)
-# G = nx.random_graphs.barabasi_albert_graph(1000, 5, 10)
+# G = nx.random_graphs.watts_strogatz_graph(1000, 4, 0.3)
+G = nx.random_graphs.barabasi_albert_graph(100, 5, 10)
 # G = nx.random_graphs.powerlaw_cluster_graph(1000, 10, 0.2)
 # G = nx.convert_node_labels_to_integers(nx.davis_southern_women_graph())
 # douban = nx.read_edgelist(os.path.dirname(os.path.realpath(__file__))+'/dataset/ASU/Douban-dataset/data/edges.csv', delimiter=',', nodetype=int, data=False)
-# G = nx.relabel_nodes(douban, {len(douban): 0}, copy=False)
+# G = nx.relabel_nodes(douban, {len(douban): 0}, copy=False)  #数据从1开始标号，需要转换为0开始记号
 
 # 网络结构绘图
 # pos=nx.spring_layout(G)
@@ -39,17 +39,19 @@
 # u = update.BirthDeath()
 # u = update.DeathBirth()
 u = update.Fermi()
+# u = update.HeteroFermi()
 
 # 演化
 e = evolution.Evolution(G, g, u)
-e.evolve(20000)
+e.evolve(2000)
 
 
 # 重复实验，得到关系图
 def repeat(times):
     a = [0] * times
     for i in xrange(100):
-        e.evolve(1000, 100)
+        e.evolve(20000, 20000)
+        print(i)
         a[i] = e.proportion[-1]
 
 # 共演

evolution.py

@@ -68,6 +68,46 @@ def evolve(self, turns, profile=None):
             if self.generation % profile == 0:
                 print('turn:'+str(self.generation))
 
+    # 时间同步演化
+    def evolve_syn(self, turns, profile=None):
+        # 初始化
+        self.strategy = np.random.randint(2, size=self.size)
+        self.fitness = np.empty(self.size, dtype=np.double)
+        # 演化记录
+        self.proportion = [0] * turns
+        # 输出间隔
+        if profile is None:
+            profile = turns/10
+            if profile < 1:
+                profile = 10
+        # 循环
+        death = None
+        for i in xrange(turns):
+            self.game.play(self.population, self.strategy, self.fitness)
+            (birth, death) = self.rule.update(self.population, self.fitness)
+
+            if self.has_mut and np.random.random() <= 0.01:
+                new_strategy = np.random.randint(2)
+            else:
+                new_strategy = self.strategy[birth]
+
+            # 统计
+            if i == 0:
+                self.proportion[0] = (self.strategy == 0).sum()
+            else:
+                self.proportion[i] = self.proportion[i - 1] + self.strategy[death] - new_strategy
+
+            # 更新策略
+            if self.strategy[death] == new_strategy:
+                death = []
+            else:
+                self.strategy[death] = new_strategy
+
+            # 记录总演化轮数
+            self.generation += 1
+            if self.generation % profile == 0:
+                print('turn:'+str(self.generation))
+
     def show(self):
         plt.figure(1)
         plt.plot(self.proportion)
@@ -157,7 +197,7 @@ def evolve(self, turns, profile=None):
             self.strategy[death] = new_s
             self.evolve_strategies[death] = new_s_e
 
-            self.coevolve.rewire_new(self.population, self.evolve_strategies[death], death)
+            self.coevolve.rewire_one(self.population, self.evolve_strategies[death], death)
 
             if (i+1)%profile == 0:
                 print('turn:'+str(i+1))

game.py

@@ -1,4 +1,6 @@
 # -*- coding: utf-8 -*-
+# -*- Author: shaodan -*-
+# -*-  2015.07.11 -*-
 
 import numpy as np
 import networkx as nx

update.py

@@ -1,4 +1,6 @@
 # -*- coding: utf-8 -*-
+# -*- Author: shaodan -*-
+# -*-  2015.07.11 -*-
 
 import numpy as np
 import networkx as nx
@@ -62,7 +64,29 @@ def update(self, graph, fitness):
         birth, death = np.random.randint(size, size=2)
         while birth == death or (not graph.has_edge(birth, death)):
             birth, death = np.random.randint(size, size=2)
-        if 1 / (1+ np.exp((fitness[birth]-fitness[death])/self.K)) < np.random.random():
+        if 1 / (1+np.exp((fitness[death]-fitness[birth])/self.K)) > np.random.random():
+            death = birth
+        return birth, death
+
+
+class HeteroFermi(UpdateRule):
+
+    def __init__(self, delta):
+        # delta = max(T, R) - min(S,P) > 0
+        # for pd delta = T-S
+        # for sd delta = T-P
+        # for sh delta = R-S
+        self.delta = delta
+
+    def update(self, graph, fitness):
+        size = len(graph)
+        # choice random pair in graph
+        birth, death = np.random.randint(size, size=2)
+        while birth == death or (not graph.has_edge(birth, death)):
+            birth, death = np.random.randint(size, size=2)
+        degree_b = graph.degree(birth)
+        degree_d = graph.degree(death)
+        if (fitness[death]-fitness[birth]) / (self.delta * max(degree_b, degree_d)) > np.random.random():
             death = birth
         return birth, death