算法_Simulated-Annealing_Python
本文用到的包
import random
import numpy as np
from numpy import linalg as LA
import pylab as plt
import networkx as nx
import sys
Simulated-Annealing (模拟退火) 是一种常见的优化方法。这里我们用它将一个任意维度的网络嵌入二维空间。为了检验算法效果，我们先制造一个二维网络A，从这个网络上我们可以得到任意两点之间的欧式距离。我们重新生成一个同等规模的随机网络C，利用C网络两点之间距离与A网络对应距离的差异，来促使C不断朝A演化，经过一千步得到B。可以看出效果非常好，几乎完全恢复了已有网络的结构。
图1. 使用SA方法恢复网络结构
SA的思路来源于统计力学中的波尔兹曼分布。该分布假设，一个多粒子系统的态出现的概率与其总能量反比，同时这种反比关系受到温度约束。
![Eq. 1][1]
[1]: http://latex.codecogs.com/svg.latex?P({state})\propto{e^{-\frac{E}{kT}}}
上式中E代表系统能量，k是一个非常小的常数，T是系统的温度。我们在其他系统中，如果将要最小化的某个目标（往往是sum of squares）看做能量，则玻尔兹曼分布给我们建议了一种梯度下降的具体方法，因为从公式1我们可以推得
![Eq. 2][2]
[2]: http://latex.codecogs.com/svg.latex?\frac{P({state2})}{P({state1})}=e^{\frac{E_1-E_2}{kT}}
于是，在实际优化过程中，我们先可以对系统做一小调整，在本例中就是比照目标网络对初始网络中的节点位置进行微调，计算这个调整带来的“能量”差异，然后通过Eq.2来判断这个调整是否应该被接受。
与死板地接受所有带来能量下降的调整的优化方法不同，模拟退火可以保证系统有一定的自由度：以小概率接受短期看不利，但长期可能有利的变化，使得系统不会被锁死在一个变化轨道上。所谓模拟退火，模拟指的就是早期保持较高的稳定，给系统较大的自由度接受短期不利的变化；后期降低温度，使得系统稳定地朝一个方向演化。
将上述思想转化为Python函数如下：
# -----------simulated annealing functions----------------------
# refreshing results
def flushPrint(variable):
sys.stdout.write('\r')
sys.stdout.write('%s' % variable)
sys.stdout.flush()
# calculate variance contributed by a node or all nodes
def varaince(newPos,Dis,node,nodePosition):
if len(nodePosition) == 2:
return sum([(LA.norm(newPos[j]-nodePosition)-Dis[(min(node,j),max(node,j))])**2 for j in newPos if j!=node])
else:
return 2*sum([(LA.norm(newPos[j]-newPos[i])-Dis[(i,j)])**2 for i,j in Dis])
def moveOneNode(newPos,Dis):
k = 1.0/len(newPos) #factor to normalize force
[i] = random.sample(newPos,1)
newpos_i = np.array(list(newPos[i]))# prevent Python from connecting the value of newPos[i] with newpos_i
varianceBefore = varaince(newPos,Dis,i,newpos_i)
for j in newPos:
newDis = LA.norm(newPos[j]-newpos_i)
if j!=i and newDis!=0:
oldDis = Dis[(min(i,j),max(i,j))]
deltaDis = newDis - oldDis # magnitude and direction of the force j imposing on i
newpos_i += (newPos[j]-newPos[i])*k*deltaDis/newDis # j pushes/pulls i
varianceAfter = varaince(newPos,Dis,i,newpos_i)
deltaE = varianceAfter-varianceBefore
return i,newpos_i,deltaE
def SA(newPos,Dis,Nsimulate,BoltzmanConstant,Temperature,coolingSpeed):
v=varaince(newPos,Dis,"null","null")
V=[v]
for k in range(Nsimulate):
flushPrint(k)
i,newpos_i,deltaE = moveOneNode(newPos,Dis)
if np.random.rand()<np.exp(-deltaE/(Temperature*BoltzmanConstant)):
newPos[i] = newpos_i
Temperature *= coolingSpeed
v += deltaE
V.append(v)
return newPos,V
现在制造两个网络，一个用于作为优化目标，另一个作为初始状态
# -----------initialization----------------
# generating a random network in 2D space #
N = 20
Pos = dict((i,np.random.rand(2)) for i in range(N))
Dis = {}
E = nx.Graph()
for i in Pos:
for j in Pos:
if i<j:
dis = LA.norm(Pos[j]-Pos[i])
Dis[(i,j)] = dis
if dis < 0.4:
E.add_edge(i,j,weight=dis)
# generate new 2d nodes
L = max(Dis.values()) #character length of the universe
newPos = dict((i,np.random.rand(2)*L) for i in range(N))
newPosCopy = newPos.copy() # copy it for plotting initial state
本例中因为系统规模较小，可以设定参数进行如下演化
newPos,V=SA(newPos,Dis,Nsimulate=1000, BoltzmanConstant=.01, Temperature=1, coolingSpeed = .9)
先定义绘制网络的函数
def plotG(ax,positions,edges):
#nx,ny = zip(*nodes.values())
for i in positions:
xi,yi = positions[i]
ax.scatter(xi,yi,s=300,facecolor='RoyalBlue',edgecolor='gray')
ax.text(xi,yi,i,ha='center',va='center',size=10,color='white')
for i,j in edges:
x1,y1 = positions[i]
x2,y2 = positions[j]
ax.plot([x1,x2],[y1,y2],marker='',linestyle='-',color='gray',alpha=0.3)
最后通过以下代码得到图1
fig = plt.figure(figsize=(10, 10))
ax1 = fig.add_subplot(221)
plotG(ax1,Pos,E.edges())
plt.text(.9,.9,'A',size=14)
ax2 = fig.add_subplot(222)
plt.text(.9,1.1,'B',size=14)
plotG(ax2,newPos,E.edges())
ax3 = fig.add_subplot(223)
plt.text(.9,.9,'C',size=14)
plotG(ax3,newPosCopy,E.edges())
ax3 = fig.add_subplot(224)
plt.text(900,37,'D',size=14)
plt.plot(V,marker='.',linestyle='-',color='RoyalBlue')
plt.xlabel('Simulation steps')
plt.ylabel('Variance')
#
plt.tight_layout()
plt.show()
#plt.savefig('/Users/csid/Desktop/SA.png',transparent=True)
