そんで、これがコード。ほとんど 1-pass で書いたので(略)。肝心のグラフ出力部を作っていないけど、座標とリンクは取れるからどうとでも。アルゴリズムは簡単で、リンクを二次元行列で表現する。点 i,j に対して距離は d_{i,j} = d{j,i}, d_{i,i} = 0 だからリンク行列は対角要素が 0 の対称行列になる。この行列に i 番目の要素以外が 0 の縦ベクトル (0, ... , 0,  1, 0, ..., 0) を右からかけると点 i に関するリンクベクトルが得られる。

#!/usr/bin/python
# coding: UTF-8

from numpy import *
import random
import math
import sys

DIM = 10
K = 1.0
SCALE = DIM * K
DISP_SCALE = K
BETA = 10.0

class GeoUtil:
    @classmethod
    def getDist(cls, x1, y1, x2, y2):
        return math.sqrt((x1-x2)**2 + (y1-y2)**2)

    @classmethod
    def getTrialDisplacement(cls):
        return random.uniform(-DISP_SCALE, DISP_SCALE), random.uniform(-DISP_SCALE, DISP_SCALE)

class Points:
    def __init__(self):
        self.linkMtrx = self.initLink()
        self.pts = [Point() for x in range(DIM)]

    def initLink(self):
        rndMin, rndMax = 0, 10
        mtrx = zeros((DIM, DIM))
        for l in range(DIM):
            for c in range(l):
                mtrx[l,c] = mtrx[c,l] = random.randint(rndMin, rndMax)
        return mtrx

    def getDist(self, i, j):
        (xi, yi) = self.pts[i].getCoord()
        (xj, yj) = self.pts[j].getCoord()
        return GeoUtil.getDist(xi, yi, xj, yj)

    def getDistMtrx(self):
        mtrx = zeros((DIM, DIM))
        for l in range(DIM):
            for c in range(l):
                mtrx[l,c] = mtrx[c,l] = self.getDist(l,c)
        return mtrx

    def getDisplacementMtrx(self):
        return self.linkMtrx - self.getDistMtrx()

    def getTotalEnergy(self):
        energy = 0.0
        for i in range(DIM):
            energy += self.getLocalEnergy(i)
        return energy

    def getLocalEnergy(self, i):
        Dmtrx = self.getDisplacementMtrx()
        localEnergy = 0.0
        Dlist = dot(Dmtrx, probeVec(i))
        for d in Dlist:
            localEnergy += d**2
        return localEnergy

    def tryMoving(self, i):
        oldLocalE = self.getLocalEnergy(i)
        self.pts[i].updateCoord()
        newLocalE = self.getLocalEnergy(i)

        dE = oldLocalE - newLocalE
        if dE < 0:
            if random.uniform(0,1) > math.exp(-BETA * math.fabs(dE)):
                self.pts[i].revertCoord()

    def updateSystem(self):
        for i in range(DIM):
            self.tryMoving(i)
        return self.getTotalEnergy()

class Point:
    s_num = -1
    num = 0
    def __init__(self):
        Point.s_num += 1
        self.num = Point.s_num
        self.x, self.y = random.random()*SCALE, random.random()*SCALE

    def getCoord(self):
        return self.x, self.y

    def setCoord(self, x, y):
        self.x, self.y = x, y

    def updateCoord(self):
        self.oldX, self.oldY = self.getCoord()
        dx, dy = GeoUtil.getTrialDisplacement()
        self.setCoord(self.oldX + dx, self.oldY + dy)

    def revertCoord(self):
        self.setCoord(self.oldX, self.oldY)

    def getLinkedPoints(self, mtrx):
        return dot(mtrx, probeVec(self.num))

def probeVec(index):
    vec = zeros((DIM))
    vec[index] = 1
    return vec

def main():
    pts = Points()
    for i in range(501):
        e = pts.updateSystem()/(DIM**2)
        if i % 5 == 0:
            print i, e

if __name__ == '__main__':
    argvs = sys.argv
    argc = len(argvs)
    if argc != 2:
        print 'Usage: # python %s num' % argvs[0]
        quit()
    DIM = int(argvs[1])

    main()
~