Drkcore

17062009 q-chem

Molekel

Molekelいい感じ。

molekel

17062009 q-chem GAMESS

GAMESSのPROJCT=.F.

ありがちなエラー。

 *** ERROR! THERE ARE NOT 5 OR 6 TRANS/ROT MODES NUM

そういう場合には

$statpt projct=.f. $end

と入れとく

16062009 zsh

zshの本

zshの本というのが気になるのだけど、実物をめくってみないと買えないかな。

zshの本 (エッセンシャルソフトウェアガイドブック)
広瀬雄二
技術評論社 / 3360円 ( 2009-06-17 )

16062009 Python

グラフ分割問題をsimulated annealingで

グラフつくるとことは禁断探索法と一緒。

メタヒューリスティクスの数理
久保幹雄,J. P. ペドロソ
共立出版 / ￥ 3,675 ()
在庫あり。

def evaluate(nodes, adj, sol, alpha):
    s = [0 for i in nodes]
    d = [0 for i in nodes]
    bal = 0
    for i in nodes:
        if sol[i] == 0:
            bal += 1
        else:
            bal -= 1
        for j in adj[i]:
            if sol[j] == sol[i]:
                s[i] += 1
            else:
                s[i] -= 1
    cost = 0
    for i in nodes: cost += d[i]
    cost /= 2
    cost += alpha*abs(bal)
    return cost, bal, s, d

def find_move_rnd(n, sol, alpha, s, d, bal):
    istar = random.randint(0,n-1)
    part = sol[istar]
    if part == 0 and bal > 0 or part == 1 and bal < 0:
        penalty = -2*alpha
    else:
        penalty = 2*alpha
    delta = s[istar] - d[istar] + penalty
    return istar, delta

def update_move(adj, sol, s, d, bal, istar):
    part = sol[istar]
    sol[istar] = 1-part
    s[istar], d[istar] = d[istar], s[istar]
    for j in adj[istar]:
        if sol[j] == part:
            s[j] -= 1
            d[j] += 1
        else:
            s[j] += 1
            d[j] -= 1
    if part == 0:
        bal -= 1
    else:
        bal += 1
    return bal

def estimate_temperature(n, sol, s, d, bal, initprob):
    ntrials = 10 * len(sol)
    nsucc = 0
    deltaZ = 0.0
    for i in range(0,ntrials):
        part = random.randint(0,1)
        istar, delta = find_move_rnd(n, sol, alpha, s, d, bal)
        if delta > 0:
            nsucc += 1
            deltaZ += delta
    if nsucc != 0:
        deltaZ /= nsucc
    T = -deltaZ/math.log(initprob)
    return T

def metropolis(T, delta):
    if delta <= 0 or random.random() < math.exp(-(delta)/T):
        return True
    else:
        return False

def annealing(nodes, adj, sol, initprob, tempfactor, sizefactor, cutoff, freezelim, minpercent, alpha):
    n = len(nodes)
    z,bal,s,d = evaluate(nodes, adj, sol, alpha)
    if bal == 0:
        solstar, zstar = list(sol), z
    T = estimate_temperature(n, sol, s, d, bal, initprob)
    nfrozen = 0
    part = 0
    while nfrozen < freezelim:
        changes, trials = 0,0
        while trials < n*sizefactor and changes < n*cutoff:
            trials += 1
            istar, delta = find_move_rnd(n, sol, alpha, s, d, bal)
            if metropolis(T, delta):
                changes += 1
                bal = update_move(adj, sol, s, d, bal, istar)
                z += delta
                if bal == 0:
                    if z < zstar and part == 0:
                        solstar, zstar = list(sol), z
                        nfrozen = 0
        T *= tempfactor
        if float(changes)/trials < minpercent:
            nfrozen += 1
    return solstar, zstar


if __name__ == "__main__":
    import networkx as nx
    import matplotlib.pyplot as plt

    num_nodes = 50
    nodes, adj = rnd_graph_fast(num_nodes,0.3)
    G=nx.Graph()
    for i in range(num_nodes):
        G.add_node(i)
    for i in range(num_nodes-1):
        for j in adj[i]:
            if i < j:
                G.add_edge(i,j)
    pos=nx.spring_layout(G)

    initprob = .4    # initial acceptance probability
    sizefactor = 16  # for det. # tentatives at each temp.
    tempfactor = .95 # cooling ratio
    freezelim = 5    # max number of iterations with less that minpercent acceptances
    minpercent = .02 # fraction of accepted moves for being not frozen
    alpha = .05      # imballance factor
    N = len(nodes)   # neighborhood size
    L = N*sizefactor # number of tentatives at current temperature 

    sol = construct(nodes)
    z,bal,s,d = evaluate(nodes, adj, sol, alpha)

    sol,z = annealing(nodes, adj, sol, initprob, tempfactor, N, L, freezelim, minpercent, alpha)
    zp,bal,sp,dp = evaluate(nodes, adj, sol, alpha)

    rnodelist = []
    bnodelist = []
    for i in range(len(sol)):
        if sol[i] == 1:
            rnodelist.append(i)
        else:
            bnodelist.append(i)
    nx.draw(G,pos,nodelist=rnodelist,node_color='r')
    nx.draw(G,pos,nodelist=bnodelist,node_color='b')
    plt.savefig("path.png")

graph