Project Euler 185

プロジェクトオイラー
http://projecteuler.net/index.php

Q185.
省略

あれこれ考えて分からなかったので、焼きなまし法を使ってみた。ある整数のエネルギーを各整数のcorrectの真値との差の平方の総和として、エネルギー最小問題とした。しかし、これは大雑把にしかわからないので、エネルギーが12以下の数が求まったら、その近くをしらみつぶしに調べてみた。

from itertools import combinations, imap

from random import randint, seed, random

from math import exp
def encode(n):

    a = [ ]

    while n:

        a.append(int(n % 10))

        n /= 10

    return a
def decode(a):

    n = 0

    for e in reversed(a):

        n = n * 10 + e

    return n
def E(ary, b):

    return sum(imap(lambda a: (a[1]

        - sum(imap(lambda x: x[0] == x[1], zip(a[0], ary)))) ** 2, b))
def calc_status_at_random(n):

    return [ randint(0, 9) for i in range(n) ]
def transit_status_at_random(s):

    i = randint(0, N - 1)

    d = randint(0, 8)

    if d < s[i]:

        s[i] = d

    else:

        s[i] = d + 1

    return s[:]     # copy
def gen_neighbor(s0):

    s = s0[:]

    for k in range(N):

        d0 = s[k]

        for d in range(10):

            if d != d0:

                s[k] = d

                yield s

        s[k] = d0
def find_detail(b, stat):

    s_set = set()

    s_set.add(tuple(stat))

    

    E_limit = 12

    while E_limit >= 0:

        new_set = set()

        for s0 in s_set:

            E0 = E(s0, b)

            for s1 in gen_neighbor(list(s0)):

                E1 = E(s1, b)

                if E1 == 0:

                    return decode(s1)

                if E1 <= E_limit:

                    t1 = tuple(s1)

                    if not t1 in s_used:

                        new_set.add(t1)

                        s_used.add(t1)

                        print s1, E1, E_limit

        s_set = new_set

        

        if E_limit <= 10:

            E_limit -= 1

        else:

            E_limit -= 2

    return -1
def find(b):

    s0 = calc_status_at_random(N)

    E0 = E(s0, b)

    T = 7.0

    while T > 0.1:

        s1 = transit_status_at_random(s0)

        E1 = E(s1, b)

        if E1 == 0:

            return decode(s1)

        if random() < exp(-(E1 - E0) / T):

            if E1 <= 13:

                print E1, s1, T

                n = find_detail(b, s1)

                if n != -1:

                    return n

            s0 = s1

            E0 = E1

        

        T -= 1e-5
a = {

    5616185650518293: 2,

    3847439647293047: 1,

    5855462940810587: 3,

    9742855507068353: 3,

    4296849643607543: 3,

    3174248439465858: 1,

    4513559094146117: 2,

    7890971548908067: 3,

    8157356344118483: 1,

    2615250744386899: 2,

    8690095851526254: 3,

    6375711915077050: 1,

    6913859173121360: 1,

    6442889055042768: 2,

    2321386104303845: 0,

    2326509471271448: 2,

    5251583379644322: 2,

    1748270476758276: 3,

    4895722652190306: 1,

    3041631117224635: 3,

    1841236454324589: 3,

    2659862637316867: 2,

}
seed(0)

b = [ ]

for n, correct in a.iteritems():

    b.append([ encode(n), correct ])
N = len(b[0][0])

M = len(b)

s_used = set()

print find(b)