「最短ヌクレオチド連鎖問題」焼きなまし法

simulated_annealing.rb

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