Skip to content

Instantly share code, notes, and snippets.

@obelisk68
Created February 14, 2019 14:18
Show Gist options
  • Save obelisk68/f5628705b33107c972ca3ba068d061a8 to your computer and use it in GitHub Desktop.
Save obelisk68/f5628705b33107c972ca3ba068d061a8 to your computer and use it in GitHub Desktop.
「最短ヌクレオチド連鎖問題」焼きなまし法
def dist(a, b)
case
when a[1, 2] == b[0, 2] then 1
when a[-1] == b[0] then 2
else 3
end
end
L = 64
start_t = 171
end_t = 66
table = %w(UUU UUC UCU CUU UUA UAU AUU UUG UGU GUU UCC CCU CUC UCA CAU AUC
UCG CGU GUC UAC ACU CUA UAA AAU AUA UAG AGU GUA UGC GCU CUG UGA
GAU AUG UGG GGU GUG CCC CCA CAC ACC CCG CGC GCC CAA AAC ACA CAG
AGC GCA CGA GAC ACG CGG GGC GCG AAA AAG AGA GAA AGG GGA GAG GGG)
nucleotide = "UCAG"
table = nucleotide.chars.repeated_permutation(3).map(&:join)
edge = table.each_cons(2).map {|a, b| dist(a, b)} + [0]
min = Float::INFINITY
1.step do |i|
a, b = rand(1...L), rand(1...L)
next if a >= b
prev = edge.sum + 3
prev_d = edge[a - 1] + edge[a] + edge[b - 1] + edge[b]
table[a], table[b] = table[b], table[a]
a_d1, a_d2 = dist(table[a - 1], table[a]), dist(table[a], table[a + 1])
if b == L - 1
b_d1, b_d2 = dist(table[b - 1], table[b]), 0
else
b_d1, b_d2 = dist(table[b - 1], table[b]), dist(table[b], table[b + 1])
end
temp = start_t + (end_t - start_t) * prev / end_t.to_f
edge[a - 1], edge[a] = a_d1, a_d2
edge[b - 1], edge[b] = b_d1, b_d2
nxt = edge.sum + 3
prob = Math.exp((prev - nxt) / temp)
if prev_d > a_d1 + a_d2 + b_d1 + b_d2 or prob > rand
if nxt < min
puts "prob = #{prob}"
min = nxt
puts "i = #{i}, length = #{nxt}"
puts table.join(" ")
end
else
table[a], table[b] = table[b], table[a]
edge[a - 1], edge[a] = a_d1, a_d2
edge[b - 1], edge[b] = b_d1, b_d2
end
end
#http://obelisk.hatenablog.com/entry/2018/03/14/003028
#参考
#http://shindannin.hatenadiary.com/entry/20121224/1356364040
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment