L-system で描いた再帰曲線の例

obelisk.hatenablog.com
前回のエントリで Ruby で実装してみた L-system を使って、実際に再帰曲線を描いてみます。Gem 'kaki-lsystem' が必要です。

C曲線。

コード。

require 'kaki/lsystem'

l = Lsystem.new(600, 500, "C curve")
l.move(-130, 100)
l.prologue do
  clear(color(0x2b5, 0x60ff, 0x5d89))
  color(0xfc9e, 0xe0aa, 0x7ac7)
end
l.set("-") {left(45)}
l.set("+") {right(45)}
l.set("F") {forward(8)}
l.init("F")
l.rule("F", "+F--F+")
l.draw(10)

以下、require 'kaki/lsystem' は書きません。

コッホ・アイランド。

コード。

l = Lsystem.new(600, 600, "Koch island")
l.move(-150, -150)
l.prologue do
  clear(color(0xffff, 0xf600, 0xf540))
  color(0x28a5, 0xa3c8, 0xca7)
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("F") {forward(8)}
l.init("F-F-F-F")
l.rule("F", "F+FF-FF-F-F+F+FF-F-F+F+FF+FF-F")
l.draw(2)

ペンローズ・タイル。

コード。

l = Lsystem.new(600, 600, "Penrose Tiling")
l.prologue do
  clear(color(0xe013, 0xfee3, 0xfee3))
  color(0xffff, 0x5feb, 0xe51f)
end
l.set("-") {left(36)}
l.set("+") {right(36)}
walk = proc {forward(10)}
l.set("M", &walk)
l.set("N", &walk)
l.set("O", &walk)
l.set("P", &walk)
l.set("A", &walk)
l.init("[N]++[N]++[N]++[N]++[N]")
l.rule("M", "OA++PA----NA[-OA----MA]++")
l.rule("N", "+OA--PA[---MA--NA]+")
l.rule("O", "-MA++NA[+++OA++PA]-")
l.rule("P", "--OA++++MA[+PA++++NA]--NA")
l.rule("A", "")
l.draw(5)

Sierpinski arrowhead。続けると「シェルピンスキーのギャスケット」になります。

コード。

n = 6
length = n.zero? ? 500 : 500.0 / 2 ** n

l = Lsystem.new(520, 500, "Sierpinski arrowhead")
l.move(-250, -240)
l.prologue do
  clear(color(0xffff, 0xef5b, 0xf677))
  color(0x12de, 0x5320, 0xb4a)
end
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("L") {}
l.set("R") {}
l.set("F") {forward(length)}
l.init("LF")
l.rule("L", "-RF+LF+RF-")
l.rule("R", "+LF-RF-LF+")
l.rule("F", "")
l.draw(n)

　
　

コード。

l = Lsystem.new(500, 500)
l.move(50, -230)
l.prologue do
  clear(color(0xd467, 0xfbab, 0xd0c7))
  color(0xe713, 0x5742, 0xff23)
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("F") {forward(10)}
l.init("F-F-F-F")
l.rule("F", "FF-F-F-F-F-F+F")
l.draw(3)

　
　
草っぽいの。

コード。

l = Lsystem.new(500, 500)
l.move(-150, -240)
l.prologue do
  clear(color(0xf7d3, 0xffff, 0xc99b))
  color(0x30e3, 0x8636, 0x10d2)
end
l.dir = Vector[0, 1]
l.set("-") {left(10)}
l.set("+") {right(15)}
l.set("F") {forward(13)}
l.set("X") {}
l.set("Z") {}
l.init("X")
l.rule("F", "FX[FX[+XF]]")
l.rule("X", "FF[+XZ++X-F[+ZX]][-X++F-X]")
l.rule("Z", "[+F-X-F][++ZX]")
l.draw(4)

ペアノ・ゴスペル曲線。

コード。

l = Lsystem.new(600, 550, "Peano-Gosper curve")
l.move(80, 250)
l.prologue do
  clear(color(0x5968, 0x5968, 0x5968))
  color(0x8024, 0xffff, 0x653)
end
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("F") {forward(9)}
l.set("X") {}
l.set("Y") {}
l.init("X")
l.rule("X", "X+YF++YF-FX--FXFX-YF+")
l.rule("Y", "-FX+YFYF++YF+FX--FX-Y")
l.draw(4)

再帰曲線を描く言語「L-system」を Ruby で実装した

L-system って何でしょうか。Wikipedia にはこうあります。

L-system（エルシステム、Lindenmayer system）は形式文法の一種で、植物の成長プロセスを初めとした様々な自然物の構造を記述・表現できるアルゴリズムである。自然物の他にも、反復関数系（Iterated Function System; IFS）のようないわゆる自己相似図形やフラクタル図形を生成する場合にも用いられる。

https://ja.wikipedia.org/wiki/L-system

つまりは、再帰曲線を描く言語ですね。例としてまず「コッホ曲線」を描いてみましょう。
まず、前に一定の長さ進むコマンドを「F」とします。で、開始手続きも「F」としましょう。これだけで

と直線が一本引かれます。
つぎに、再帰的な「置き換え」の規則を定めます。それを「F→F-F++F-F」とします。これは開始手続きの「F」を「F-F++F-F」で置き換えるということになります。「-」は左60度、「+」は右60度の回転とします。よって、新たな手続きも「F-F++F-F」となります。これを図示すると

となります。
置き換えを続けましょう。今度は前の手続き「F++F-F」の F をすべて F++F-F に置き換えます。すると「F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F」となります（ややこしいですね）。これは

となります。
もう一度置き換えてみます。今度は「F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F-F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F-F-F++F-F-F-F++F-F++F-F++F-F-F-F++F-F」ですよ。これは

となりますね。まさしく再帰的な「コッホ曲線」が描けています。

さて、これはこんなコードで描けます。Gem化しているので、Gem 'kaki-lsystem' をインストールして下さい。

require 'kaki/lsystem'

n = 3
length = 500.0 / 3 ** n

l = Lsystem.new(550, 200)
l.move(-250, -70)
l.set("-") {left(60)}
l.set("+") {right(60)}
l.set("F") {forward(length)}
l.init("F")
l.rule("F", "F-F++F-F")
l.draw(n)

init("F") は開始手続きの指定、rule("F", "F-F++F-F") は置き換えの規則の指定です。draw(3) で変換を 3回施したものを描画します。

set("-") {left(60)} はコマンド「-」が何をするのか、ブロック内に記述します。つまり left(60) は左60度回転ということですね。forward(20) ならまっすぐ 20 進むということです（ちなみに forward(20, false) なら、描画せずに 20 進みます）。このブロック内は、じつは Gem 'oekaki' のタートルグラフィックスの機能を使っているので、そのメソッドを書きます。よく使うのは left(deg), right(deg), forward(length, draw=true) あたりでしょうか。これについては

• タートルグラフィックスで再帰曲線をいろいろ描いてみる（Ruby） - Camera Obscura
• File: README — Documentation for kaki-lsystem (0.0.3)

などを参考にして下さい。なお、set() は使用したすべてのコマンドについて設定して下さい。何もしない場合でも空のブロックが必要です。また、rule() で変換規則が指定されていないコマンドは、そのまま持ち越されます。

置き換えるコマンドを二種類以上にすることもできます。例えば

require 'kaki/lsystem'

l = Lsystem.new(400, 350, "Dragon curve")
l.move(-50, -50)
l.prologue do
  clear(color(0xf7d3, 0xffff, 0xc99b))    #背景色
  color(0x143f, 0x832d, 0xe783)           #ペンの色
end
l.set("-") {left(90)}
l.set("+") {right(90)}
l.set("A") {forward(20)}
l.set("B") {forward(20)}
l.init("A")
l.rule("A", "A+B+")
l.rule("B", "-A-B")
l.draw(7)

では、開始手続きが A で、置き換えの規則は「A→A+B+」「B→-A-B」と二種類になっています。これを実行すると

となります。いわゆる「ドラゴン曲線」ですね。他のメソッドも簡単に説明しておきましょう。move(x, y) はペンの開始位置を指定します。デフォルトは (0, 0) で、これはキャンバスの中心です。なお、数学のグラフと同じで y方向は上がプラスの向きです。prologue {...} はブロック内で初期処理をします。背景色や色の指定をすることが多いです。なお、ここでは残念ながら left(), right() が使えない（というか2度呼ばれてしまうので正確に指定できません）ので、ペンの最初の方向はアクセサーメソッド dir= で指定することになります（あとで説明します）。

コマンド「[」「]」でスタックのプッシュ、ポップを表します。つまり [ で状態を保存し、 ] で元の状態に復帰します。例えば

require 'kaki/lsystem'

l = Lsystem.new(400, 640)
l.move(50, -310)
l.prologue {color(0x28a5, 0xa3c8, 0xca7)}
l.dir = Vector[0, 1]
l.set("-") {left(15)}
l.set("+") {right(15)}
l.set("F") {forward(10)}
l.init("F")
l.rule("F", "F[+F-F-F]F[--F+F+F]")
l.draw(4)

こんな風です。dir= で開始時の方向をベクトルで指定しています。ただし、ベクトルは単位ベクトル（長さが 1 ）になるようにして下さい。この結果は

となります。左からそれぞれ 1, 2, 4 回の繰り返しになります。

ブロックはクロージャなので、変数を使ったサンプル。「ヒルベルト曲線」の色を描画しながら替えていきます。

require 'kaki/lsystem'

n, f = 0, 0
col = [[0xffff, 0, 0], [0, 0xffff, 0]]

l = Lsystem.new(500, 500, "Hilbert curve")
l.move(-230, 230)
l.set("+") {left(90)}
l.set("-") {right(90)}
l.set("A") {}
l.set("B") {}
l.set("F") do
  if n.zero?
    f = (f + 1) % 2
    color(*col[f])
  end
  forward(15)
  n = (n + 1) % 16
end
l.init("A")
l.rule("A", "-BF+AFA+FB-")
l.rule("B", "+AF-BFB-FA+")
l.draw(5)

結果。

　

Gem 'kaki-lsystem' のソースコードは Gist にあります。この Gem は内部で Gem 'oekaki' を使っています。
https://gist.github.com/obelisk68/1c3d278cb5fe2778a11ba719d44af743

※参考にしたページ（ありがとうございます！）
L-system
L-system入門 - 人工知能に関する断創録
HiiragiCompany -L-System Tips-
L-systemを利用した植物の描画 | Makoto Hirahara
　
　
※L-system の描画例
obelisk.hatenablog.com

ヒープソート（C言語、Go言語）

（追記）ヒープを作るのに、何だかふつうとはちがう（オリジナルの？）アルゴリズムを採用してしまったようですね。これでもちゃんと動作しますが、ふつうの方法よりも多少効率が悪いようです。最後にふつうの実装も載せます。

---

　
ヒープソートですが、まず「（二分木）ヒープ」を作ります。（ここでいうヒープは「ヒープ領域」とは関係ありません。）

ヒープは二分木構造をもっていて、どこを取っても2つの（あるいは最後だけ1つの）子の値よりも、親の値の方が大きいようになっているものです。ヒープはいろいろな使い道があるそうです。まずはヒープを作ってみます。結構むずかしいです。

1. 親のインデックスを parent とすると、左側の子のインデックスは parent * 2 + 1 になります。
2. 子が存在しない場合は、親のインデックスをひとつ減らして最初から同じことをします。
3. 子が存在するばあいは、左の子と右の子で値の大きい方を選び、それを larger とします。
   • 親の方が larger 以上ならば、親のインデックスをひとつ減らして最初から同じことをします。
   • そうでなければ 親と larger の位置を交換します。そして、その降りてきた親を新しい親として、最初から同じことをします。

終了条件は、親のインデックスが負になった場合です。開始時の親は、配列の最後のその親になります。これのインデックスは、配列の長さを len とすると、(len - 2) / 2 で表されます。

これを C で表現してみましょう。

void build_heap(int ar[], int last) {
    int child;
    int parent = (last - 1) / 2;
    
    while (parent >= 0) {
        child = parent * 2 + 1;
        if (child > last) {
            parent--;
            continue;
        }
        if (child != last) && ar[child + 1] > ar[child]) child++;
        if (ar[parent] >= ar[child]) {
            parent--;
            continue;
        }
        swap(ar, parent, child);
        parent = child;
    }
}

swap() は値の交換、last は配列の右端のインデックスです。これでヒープを作ることができました。
例えば [4, 1, 6, 2, 9, 7, 3, 8] からヒープを作ると、[9, 8, 7, 4, 1, 6, 3, 2] となります。こういう二分木ですね。

9 -- 8 -- 4 -- 2
       -- 1
  -- 7 -- 6
       -- 3

9 の子が「8, 7」、8 の子が「4, 1」、7 の子が「6, 3」、4 の子が「2」ということです。どこを見ても親の方が大きくなっています。

なお、これは上の方の値が大きくなるので、正確には「最大ヒープ」というものです。上の方ほど小さいヒープは、「最小ヒープ」と呼ばれます。
　

では、これを使ってソートしてみましょう。

1. まずヒープを作ります。
2. ヒープの頭（配列の先頭）が最大値なので、これと配列の最後を入れ替えます。
3. 配列の最後（最大値）を除いたものを再帰的にソートして、除かれたものを最後に付け加えます。

これがヒープソートです。C で表現するとこんな感じです。

void heap_sort(int ar[], int len) {
    if (len <= 1) return;
    build_heap(ar, len - 1);
    swap(ar, 0, len - 1);
    heap_sort(ar, len - 1);
}

　
全体のコードとしてはこうなります。
heap_sort.c

#include <stdio.h>

void swap(int ar[], int a, int b) {
    int tmp = ar[a];
    ar[a] = ar[b];
    ar[b] = tmp;
}

void build_heap(int ar[], int last) {
    int child;
    int parent = (last - 1) / 2;
    
    while (parent >= 0) {
        child = parent * 2 + 1;
        if (child > last) {
            parent--;
            continue;
        }
        if (child != last && ar[child + 1] > ar[child]) child++;
        if (ar[parent] >= ar[child]) {
            parent--;
            continue;
        }
        swap(ar, parent, child);
        parent = child;
    }
}

void heap_sort(int ar[], int len) {
    if (len <= 1) return;
    build_heap(ar, len - 1);
    swap(ar, 0, len - 1);
    heap_sort(ar, len - 1);
}

int main() {
    int ar[] = {4, 1, 6, 2, 9, 7, 3, 8};
    int len = sizeof ar / sizeof(int);
    
    heap_sort(ar, len);
    
    printf("[");
    for (int i = 0; i < len; i++) {
        printf("%d, ", ar[i]);
    }
    printf("\b\b]\n");
    
    return 0;
}

Go言語版。C の場合とほぼ同じです。
heap_sort.go

package main
import "fmt"

func build_heap(ar []int) []int {
    last := len(ar) - 1
    parent := (last - 1) / 2
    for parent >= 0 {
        child := parent * 2 + 1
        if child > last {
            parent--
            continue
        }
        if child != last && ar[child + 1] > ar[child] {child++}
        if ar[parent] >= ar[child] {
            parent--
            continue
        }
        ar[parent], ar[child] = ar[child], ar[parent]
        parent = child
    }
    return ar
}

func heap_sort(ar []int) []int {
    len := len(ar)
    if len <= 1 {return ar}
    ar = build_heap(ar)
    ar[0], ar[len - 1] = ar[len - 1], ar[0]
    return append(heap_sort(ar[:len - 1]), ar[len - 1])
}

func main() {
    ar := heap_sort([]int{4, 1, 6, 2, 9, 7, 3, 8})
    fmt.Println(ar)
}

ふつうの実装

C言語。
heap_sort1.c

#include <stdio.h>

void heapify(int ar[], int parent, int last) {
    int child, v;
    
    v = ar[parent];
    while (1) {
        child = parent * 2 + 1;
        if (child > last) break;
        if (child != last && ar[child + 1] > ar[child]) child++;
        if (v >= ar[child]) break;
        ar[parent] = ar[child];
        parent = child;
    }
    ar[parent] = v;
}

void build_heap(int ar[], int last) {
    for (int i = (last - 1) / 2; i >= 0; i--)
        heapify(ar, i, last);
}

void heap_sort(int ar[], int len) {
    if (len <= 1) return;
    build_heap(ar, len - 1);
    int tmp = ar[0]; ar[0] = ar[len - 1]; ar[len - 1] = tmp;
    heap_sort(ar, len - 1);
}

int main() {
    int ar[] = {4, 1, 6, 2, 9, 7, 3, 8};
    int len = sizeof ar / sizeof(int);
    
    heap_sort(ar, len);
    
    printf("[");
    for (int i = 0; i < len; i++) {
        printf("%d, ", ar[i]);
    }
    printf("\b\b]\n");
    
    return 0;
}

　
Go言語版。
heap_sort1.go

package main
import "fmt"

func build_heap(ar []int) []int {
    heapify := func (parent int) {
        last := len(ar) - 1
        v := ar[parent]
        for {
            child := parent * 2 + 1
            if child > last {break}
            if child != last && ar[child + 1] > ar[child] {child++}
            if v >= ar[child] {break}
            ar[parent] = ar[child]
            parent = child
        }
        ar[parent] = v
    }
    
    for i:= (len(ar) - 2) / 2; i >= 0; i-- {heapify(i)}
    return ar
}

func heap_sort(ar []int) []int {
    len := len(ar)
    if len <= 1 {return ar}
    ar = build_heap(ar)
    ar[0], ar[len - 1] = ar[len - 1], ar[0]
    return append(heap_sort(ar[:len - 1]), ar[len - 1])
}

func main() {
    ar := heap_sort([]int{4, 1, 6, 2, 9, 7, 3, 8})
    fmt.Println(ar)
}

※参考
http://www.th.cs.meiji.ac.jp/assets/researches/2005/omoto/heapsort.html

バブルソート（C言語、Go言語）

バブルソートはわかりやすいソートです。隣どうしを比較して、前の方が値が大きければ値を交換します。それを繰り返し行うことでソートします。ただ、わかりやすく実装も簡単ですが、効率はよくありません。ソートのアルゴリズムとしてはまず使わないのではないでしょうか。アルゴリズム的にもつまらないですね。

バブルソートの名前は、値が泡（バブル）のように上へ登っていくからでしょうね。値が下がっていくようにも実装できますが、せっかくなので名前どおりにしてみました（笑）。

C言語版。
bubble_sort.c

#include <stdio.h>

void bubble_sort(int ar[], int len) {
    int i, j, tmp;
    
    for (i = len - 2; i >= 0; i--)
        for (j = 0; j <= i; j++)
            if (ar[j] > ar[j + 1]) {
                tmp = ar[j]; ar[j] = ar[j + 1]; ar[j + 1] = tmp;
            }
}

int main () {
    int ar[] = {3, 1, 5, 2, 9, 7, 0};
    int len = sizeof ar / sizeof(int);
    
    bubble_sort(ar, len);
    
    printf("[");
    for (int i = 0; i < len; i++) {
        printf("%d, ", ar[i]);
    }
    printf("\b\b]\n");
    
    return 0;
}

　
Go言語版。C の場合とほとんど同じで、工夫の余地はありません。
bubble_sort2.go

package main
import "fmt"

func bubble_sort(ar []int) []int {
    for i := len(ar) - 2 ; i >= 0; i-- {
        for j := 0; j <= i; j++ {
            if ar[j] > ar[j + 1] {
               ar[j], ar[j + 1] = ar[j + 1], ar[j]
            }
        }
    }
    return ar
}

func main() {
    fmt.Println(bubble_sort([]int{5, 2, 1, 8, 6, 3, 0}))
}

おまけ。Ruby 版。

class Array
  def bubble_sort
    (size - 2).downto(0) do |i|
      (i + 1).times do |j|
        self[j], self[j + 1] = self[j + 1], self[j] if self[j] > self[j + 1]
      end
    end
    self
  end
end

なお、これらの実装では与えられた配列が破壊的に変更されます。

ジュリア集合を描いてみる（Ruby）

ジュリア集合（Julia set）を Ruby で描いてみました。

　
描画には自作の Gem 'oekaki' を使っています。
oekaki | RubyGems.org | your community gem host
GTK+でお絵かきしてみた（Ruby） - Camera Obscura
　
コード。
julia_set.rb

require 'oekaki'

Width, Height = 600, 400

Oekaki.app width: Width, height: Height, title: "Julia set" do
  draw do
    clear
    red   = color(0xffff, 0, 0)
    green = color(0, 0xffff, 0)
    blue  = color(0, 0, 0xffff)
    
    a, jmid = 4.0 / Width, Height / 2
    Width.times do |i|
      jmid.times do |j|
        x, y = a * i - 2.0, a * j
        catch(:jump) do
          16.times do
            x2, y2 = x * x, y * y
            d = x2 + y2
            break if d < 1e-10
            d = d * d
            x, y = (1.0 / 3) * (2 * x + (x2 - y2) / d), (2.0 / 3) * y * (1 - x / d)
            ya = sqrt(3) / 2
            
            drawing = proc do |x1, y1, c1, c2|
              if x1 * x1 + y1 * y1 < 0.0025
                point(i, jmid + j, c1)
                point(i, jmid - j, c2)
                throw(:jump)
              end
            end
            
            drawing[x - 1, y, red, red]
            drawing[x + 0.5, y + ya, green, blue]
            drawing[x + 0.5, y - ya, blue, green]
          end
        end
      end
    end
  end
end

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

これおもしろい本です。おすすめ！

『Cプログラミング診断室』を読んだ

改訂新版 Cプログラミング診断室

たまたま BOOK OFF にあったので買って読んでみました。かなり古い本で、一度品切れになってから著者がネット上で全文を公開して、じつはこのブログの過去記事でそのサイトが紹介してあります。
obelisk.hatenablog.com
ただ、自分はネットでは長文が読めない体質（？）なので、紙媒体で読んでみてよかったです。

内容は過去記事にも書きましたが、C言語で書かれたいわゆる「クソコード」（笑）を眺めて、怒りにうち震えつつリファクタリングするという、身も蓋もない本です。この改訂新版が出てからでも既に 15年が経っていますから、C言語はさすがに長く使われていますね。中身は非常におもしろかったし、ためにもなりました。本書ではいろいろなことが指摘されていますが、特によく出てくる怒り（笑）としては、「変数や関数の名前はとても大切」「長い関数を書くな」「コピペをするな、重複したコードを書くな」「ハードコーディングをするな」「何でも自分で書かないで、よいライブラリを使え」というところでしょうか。これらは C に限らないですね。

また、著者は「閉じこもっていないで、よき仲間・助言者を持て」としきりに強調します。本だけではダメだ、独学では初心者プログラマから脱却することはむずかしいと。これなどは自分としては、ため息が出ますね。僕はプログラミングは完全に独学なので、おそらくはひどいコードを書いているのだろうなと。でも、田舎に住んでいるとなかなかむずかしいです。思うに、田舎に住んでいて誰も相談するひとはいないけれど、素人プログラミングしてみたいという自分みたいな人は一定数いるのではないでしょうか。そういう人たちが頼りにできる場所はいまでも少ないですね。このブログは密かに、レヴェルは低いけれどプログラミングしてみたいという（自分のような）人に少しでも参考にならないかと思ってやっているところもあります。

で、そういう人のためにも本書は役立ちますよ。具体的に「クソコード」（本当にひどいものです）を掲載して、なるたけ添削して下すっているのですから。

ただ、つらつら思ったのですが、C言語はこれからも使われ続けるでしょうが、C言語つらいなあとも感じました。本書でもすべて C で書く必要はないとして、（古い本なので）これからは Perl がいいのではないかとあります。いまなら Ruby とか Python とかになるでしょうか。とにかく C は、きちんと書けるようになるまですごく時間がかかる。それなので、Go言語などは C の代替として考えられたのでしょうね。僕のような初心者プログラマは、C よりも Go の方がずっとラクです。低レヴェル（ハードウェアに近いレイヤという意味です）プログラミングは Go に置き換わればいいと思うのだが、そうはならないだろうな。いや、脱線しましたが、とにかく本書はおもしろかったです。ある程度 C のコードが読めるようになったら、是非おすすめです。

ちなみに、どうでもいいですが、このブログの最初の方のコードは自分でも碌なものではないと思うものが多いです。そういうのに検索で人がきたりするので、ときどき書き直したりしています。でも、独学素人プログラマは、プログラミング用のブログをもつのは絶対いいことだと思っています。そこで恥をさらせばいいのです。

プログラミング言語C 第2版 ANSI規格準拠

「最短ヌクレオチド連鎖問題」が解けない

高校生のとき、友人にこんな問題を出されました。適当に再構成してみます。

タンパク質を構成するアミノ酸は 20種類あるよね。DNA または RNA はそのアミノ酸の設計図で、いわゆる「コドン」を指定している。コドンというのはヌクレオチドの塩基３個から成るもので、ヌクレオチドの塩基は 4種類ある。RNA の場合なら U（ウラシル）、C（シトシン）、A（アデニン）、G（グアニン） だ。これら 4種類の ３個づつの重複順列でコドンは出来ているから、全部で 4×4×4 = 64 種類あることになる。これらコドンによってアミノ酸が指定されるのだよね。

さて、ここからは遊びだ。64種類のコドンをすべて一度づつ使って、それらを適当な順番ですべて連結したとする。このとき、実際にはあり得ないのだが、重複するものは圧縮できるものとする。例えば UUU と UUA というコドンをこの順につなげたとき、UUUUUA の UU の部分が重複しているから、圧縮できて UUUA になる。この次に例えば UAG が来れば、また 2つ圧縮できて UUUAG となる。重複は一度つなげるごとに最大でこのような 2つで、1つの場合もまたまったく重複しない場合もある。

このとき、適当な順番ですべての種類のコドンを連結して圧縮したなら、そのヌクレオチド連鎖の長さの最小値はいくらになるだろう？

　
たぶん、高校の生物の先生が冗談で出した問題だったと思う。そのとき自分はすでに PC（富士通 FM-7！）をもっていたので、機械語でそれを解くプログラムを考えてみたのだが、そのときはうまくいかなかった。学校ではかしこい生徒がかなり短い例を作ったらしいが、その答えは我々にはわかっていない。そしてそれが正解かどうかもわからない。

あれからもう何十年も過ぎたけれど、この問題はいまでも時々思い出す。そしていま素人プログラミングで遊んでいるので、つらつら考えてみるのだが、どうもわからない。これを仮に「最短ヌクレオチド連鎖問題」ということにして、この問題が自分には解けないのだ。

で、じつはそれを解く Ruby プログラムを考えてみたのだが、それはたぶんうまくない。誰かが正解を求めるコードを書いてくれないだろうか、そんな気もあってここに記した。コドンの種類は 64通りということはわかっているから、答えは 66 より小さくはならないことはわかっている。さて、どうであろうか。

なお、時間がどれだけかかってもよいというなら簡単である。用意されたコドンは完全に対称的なので、始まりは何でもよい。なので 64 - 1 = 63 の階乗とおりの順列を考え、それを圧縮して比較し、最短を出すことができる。しかし、63 の階乗というのは

1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000

なのである。これを仮に Ruby で解くとすると、自分の PC では 100万回のループでほぼ 1秒かかるというのが大まかな目安だから、おおよそ

62868097266756724509010233015027211363038232186303546045988064422017052054

年もかかってしまうのだ（細かい数値に意味はないけれど）。さて、それをどうするかということである。

おそらく、いちばん有望な解き方はこうである。まず、与えられた幾らかのコドンを連結圧縮して、それが最短になっているとする。このとき、まだ使われていないコドンをそこにひとつ追加して、そこで最短になるような連結の仕方（アルゴリズム）を与えることができるか。もしそれが可能ならば、再帰的な処理で答えを確実に求めることができる。問題は、そのようなアルゴリズムが存在するか、だ。

追記

いまのところ自分の見つけた最短連鎖の長さは 80 である。その順は例えば

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

である。これはたぶん最短でない。考え方としては、できるだけ 2個の圧縮が続くように、でなければ 1個の圧縮、となるように適当に並べただけである。最短の保証はない。

再追記（8/25）

Haskell で brute force するコードを書いてみた。これで求められる筈だけれど、もちろん永遠に終わらない。
shortest_chain_of_nucleotide.hs

import Data.List

main :: IO ()
main = do
    print $ head $ f $ map joinCodons $ permutations codons
        where nucleotide = "UCAG"
              codons = [[a, b, c] | a <- nucleotide, b <- nucleotide, c <- nucleotide]
              f a = [(length x, x) | x <- a, (minimum $ map length a) == length x]

joinCodons :: (Eq a) => [[a]] -> [a]
joinCodons (x:y:ys) = if null ys
                      then head st
                      else joinCodons st
                          where st = [compress x y] ++ ys

compress :: (Eq a) => [a] -> [a] -> [a] 
compress x y = if double x y
               then x ++ drop 2 y
               else if single x y
                    then x ++ drop 1 y
                    else x ++ y

double :: (Eq a) => [a] -> [a] -> Bool
double x y = drop (length x - 2) x == take 2 y

single :: (Eq a) => [a] -> [a] -> Bool
single x y = last x == head y

追記

これは本質的に「最短ハミルトン路」問題。

上の長さ 80 の例から単純な乱択アルゴリズムで長さ 72 まで短縮した例。

i = 2677, length = 72
UUU UUC UCU CUU UUA UAU AUU UUG UGU GUC UCC CCU CUC UCA CAU AUC
UCG CGU GUA UAC ACU CUA UAA AAU AUA UAG AGU GUG UGC GCU CUG UGA
GAU AUG UGG GGU GUU CCC CCA CAC ACC CCG CGC GCA CAA AAC ACA CAG
AGC GCG CGA GAC ACG CGG GGC GCC AAA AAG AGA GAA AGG GGA GAG GGG

Ruby コード。

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

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)
edge = table.each_cons(2).map {|a, b| dist(a, b)} + [0]

1.step do |i|
  a, b = rand(1...L), rand(1...L)
  next if a >= b
  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
  if prev_d > a_d1 + a_d2 + b_d1 + b_d2
    edge[a - 1], edge[a] = a_d1, a_d2
    edge[b - 1], edge[b] = b_d1, b_d2
    puts "i = #{i}, length = #{3 + edge.sum}"
    p table
  else
    table[a], table[b] = table[b], table[a]
  end
end

　
「焼きなまし法」でも挑戦してみたけれど、全然ダメ。一回での距離の動きの幅が小さすぎる。
「最短ヌクレオチド連鎖問題」焼きなまし法 · GitHub
（2019/2/14）