2013年11月10日日曜日

((Rで) 書く ((もっとRっぽい) Lisp) インタプリタ)

((Rで) 書く (Lisp) インタプリタ)の続きです。

以前の記事を見た@kohskeさんから、 環境はRの環境オブジェクトを使えば良いのでは、とアドバイスを頂きました。

やってみたら自前で環境を実装しなくて良い分簡潔になったし、より面白い使い方ができるようになりました。

コード

前回同様Gistに上げておきました。 変なところがあったら教えてください。

主な変更点は

  • Rの環境オブジェクトを流用
  • Lispの関数の呼び出し方法を変更: proc(exps) から do.call(proc, exps) に

遊び方

前回と同様です。

Rインタプリタを起動してコードを読み込ませます。
repl() を実行するとLispの対話式インタプリタが起動します。

$ ls
lispr.R
$ R -q
> source("lispr.R")
> repl()
lispr> (+ 1 2)
3
lispr> (define add2 (lambda (x) (+ x 2)))
#<closure>
lispr> (add2 40)
42
lispr>

関数 repl の引数 parent に渡した環境がLispの環境の親になります。(デフォルトだとRのグローバル環境)
よって今回はLispの環境からRの環境を参照できるので、LispでRの関数が使えます。

lispr> (: 1 10)
1 2 3 4 5 6 7 8 9 10
lispr> (rnorm 3)
0.536481224524994 -0.547993231580984 -2.14041393248752
lispr> (seq 0 10 2)
0 2 4 6 8 10
lispr> (begin (plot cars) (lines (lowess cars)) 1)
1
lispr> 

S式でRを使える!楽しい!

感想

やはり言語処理系を作るのは楽しいです。

2013年8月9日金曜日

((Rubyで) 書く (Lisp) インタプリタ)

最近Rubyを始めました。

手始めにLispインタプリタを書いてみます。

元ねたはPeter Norvigの (How to Write a (Lisp) Interpreter (in Python)) (日本語訳: ((Pythonで) 書く (Lisp) インタプリタ))です。

(ちなみに以前、Rでも同じことをやりました: ((Rで) 書く (Lisp) インタプリタ))

コード

Gistに上げておきました。 変なところがあったら教えてください。

遊び方

$ ruby lisp.rb 
lisp.rb> (+ 1 2 3)
6
lisp.rb> (define add2 (lambda (x) (+ x 2)))
#<Proc:0x00000001647958@lisp.rb:62 (lambda)>
lisp.rb> (add2 40)
42
lisp.rb> (define map (lambda (f l) (if (null? l) (quote ()) (cons (f (car l)) (map f (cdr l))))))
#<Proc:0x0000000164f1a8@lisp.rb:62 (lambda)>
lisp.rb> (map add2 (list 10 20 30))
(12 22 32)
lisp.rb> 

感想

Rubyを使ってみて思ったことをメモしておきます。

  • Perlっぽい

    • メソッドのかっこを省略できる
    • ifとかwhileを後置できる
  • Lisp(関数型)っぽい

    • Symbol
    • いろんなものが値を返す(ifとか)
    • returnを省略できる
  • オブジェクト指向

    • 何でもオブジェクト!
    • mapとかの高階関数がメソッドなのはなんだか違和感があるけど、 メソッドチェインができるのは良いかも。

割と良い感じです。もうちょっと使ってみます。 ワンライナーとかちょっとしたスクリプトとかPerlのかわりに使ってみようかな。

参考文献

同じことをやっている方は他にもいます。
元ねたのページ のコメントで紹介されていた 以下の実装が簡潔で素敵です。 わからないところはこちらをチラ見しながら作りました。

(追記)

その後少し手を加えて、REPLで式の途中に改行を入れられるようにしました。
(上に貼ってあるGistは修正後です。)

こんな感じに入力できます。

lisp.rb> (+ (* 3
               (+ (* 2 4)
                  (+ 3 5)))
            (+ (- 10 7)
               6))
57
lisp.rb> 

思ったより少ない変更で実現できました。 (Gistの履歴)

2013年7月22日月曜日

Common LispでR (RCLを使ってみる)

前回の続きです。

今回はRCLを使って、Common LispからRの関数を呼んでグラフを描いてみます。

インストール

こちらに書いてある通り、Quicklispでインストールできました。

> (ql:quickload :rcl)

使い方

詳しいマニュアルは無いみたいです。

RCL examples に短いサンプル、 Comparing RCL with RpyRPyとの比較があります。

あとインストール先の examples/ 以下にサンプルコードがいくつかあります。

グラフ描画のコード

サンプルコードを参考にしながら何となく書いてみました。

(asdf:oos 'asdf:load-op :rcl)
(r:r-init)
(r:enable-rcl-syntax)

(defun plot-num-of-animals (animal-state)
  (let ((day (mapcar #'car animal-state))
        (num (mapcar (lambda (x) (length (cadr x))) animal-state)))
    [plot day num :type "o" :pch 20 :ylim [c 0 [max num]]
         :main "動物の個体数" :xlab "経過日数" :ylab "動物の個体数"]))

(defun plot-hist (animal-state)
  (let ((max-x 1.0)
        (max-y 160))
    (mapc (lambda (state)
            (let ((day (car state))
                  (rate (cadr state)))
              [hist rate :breaks [seq 0 max-x 0.02]
                   :ylim [c 0 max-y] :main (format nil "~7d 日後" day)
                   :col "gray" :xlab "直進率" :ylab "動物の個体数"]))
          animal-state)))

(defun animal-state->png (out-dir animal-state)
  (r:with-device ((format nil "~A/%03d" out-dir) :png :width 400 :height 300)
    (plot-num-of-animals animal-state)
    (plot-hist animal-state)))

角括弧[] の部分でRの関数を呼んでいます。 Lispのコードと混ぜて書けるのがおもしろいです。

あと、r:with-device で描画のコードを包むのがLispっぽくて好きです。 (Rだと デバイス開く、描画、デバイス閉じる、と手続き的)

描画

描画を実行します。

まずは前回同様、20万日分シミュレーションを走らせ、1万日間隔で動物の状態をdumpします。

> (defparameter *dump* (dump-animal-state 200000 10000))

グラフを描画します。

> (animal-state->png "Rplot" *dump*)

Rplot/ 以下に連番のpngファイルが生成されます。

前回同様のグラフができました。

感想

今回試した範囲では割と良い感じに使えました。 Rに慣れた人が、Common Lispでデータをいじっていて、ちょっとグラフを描いて確かめたい、 みたいなときに便利かも。

ちなみにCommon Lispでグラフを描く方法は他にもいろいろあるみたいです。 (CLiki: plotting)

2013年6月9日日曜日

Land of Lispの進化シミュレーションの様子をRでグラフ化する

Land of Lisp の10章で進化シミュレーションを実装します。

こんな感じ。

* が草で M が動物です。
ちなみに翻訳者のShiro氏が、遺伝子の傾向の違いを色で表示する版のGaucheのコードを公開しています。 (Island Life - GaucheでもLand of Lisp)

ところでシミュレーションを続けると、行動パターンが異なる2つの種に進化します。

  • ゾウ : ジャングル(画面中央部)の豊富な食物に集中する。ジャングルから離れないように周りをうろうろする。
  • ロバ : 乾燥帯(ジャングル以外の部分)で食物を探し回る。広い範囲を移動する。

面白いですね。

この種の分化がどのように進むのか気になったので、Rを使ってグラフ化してみます。

動物の遺伝子

進化シミュレーションのコードは公式ページで公開されています。(evolution.lisp)

動物の構造体の中身はこんな感じ。

#S(ANIMAL :X 50 :Y 15 :ENERGY 114 :DIR 0 :GENES (5 4 10 9 5 5 9 10))

genes フィードが遺伝子情報で、動物の行動パターンを決定します。 8つの整数は動物が8方向どちらに向きを変えやすいかを表します。 先頭 (スロット 0) の値が大きければその動物は向きを変えずに直進しやすく、 次 (スロット 1) の値が大きければ時計回りに45度向きを変えやすい、といった感じです。

ゾウ? それともロバ?

動物の性質を表すために、 直進率 という指標を考えました。これは動物がまっすぐ進む確率です。 直進率が高い = 長距離を移動する = ロバ、直進率が低い = うろうろする = ゾウ、と考えて良いですかね。

動物の直進率を計算する straight-rate 関数:

(defun straight-rate (animal)
  (let ((genes (animal-genes animal)))
    (/ (car genes)
       (apply #'+ genes))))

動物の状態をdump

進化の過程を知りたいので、一定間隔で動物たちの直進率をdumpします。

(defun dump-animal-state (days by)
  (loop for i
     from 1
     to days
     do (update-world)
     if (zerop (mod i by))
     collect (list i
                   (mapcar #'straight-rate *animals*))))

days で指定した日数分シミュレーションを走らせ、by で指定した間隔で全動物の直進率を計算します。
この関数は (経過日数 全動物の直進率のリスト) のリストを返します。

3日間シミュレーションを走らせて1日間隔でdumpする例:

> (defparameter *dump* (dump-animal-state 3 1))
*DUMP*
> *dump*
((1 (5/23 5/23)) (2 (1/5 1/5 5/23 5/23))
 (3 (5/23 5/23 1/5 1/5 1/5 1/5 5/23 5/23)))

Rにデータを渡す

指定期間の動物の状態を計算できるようになりましたが、 Common LispとRの間でどうやってデータをやり取りしましょう。

今回はCommon Lisp側でS式をRのデータ表現に変換することにしました。 R側はそれをevalするだけなので楽です。evalは危険だけどこういう時は手軽で便利ですね。

animal-state->r 関数は動物の状態のリストをRのデータ表現に変換します。

(defun animal-state->r (animal-state)
  (format nil "list(~{~a~^,~%     ~})~%"
          (mapcar (lambda (x)
                    (format nil
                            "list(day=~a, straightRate=c(~{~f~^, ~}))"
                            (car x)
                            (cadr x)))
                  animal-state)))

先ほどの *dump* を変換してみます。

> (animal-state->r *dump*)
"list(list(day=1, straightRate=c(0.2173913, 0.2173913)),
     list(day=2, straightRate=c(0.2, 0.2, 0.2173913, 0.2173913)),
     list(day=3, straightRate=c(0.2173913, 0.2173913, 0.2, 0.2, 0.2, 0.2, 0.2173913, 0.2173913)))
"

ファイルに書き出す関数:

(defun write-animal-state (fname animal-state)
  (with-open-file (my-stream
                   fname
                   :direction :output
                   :if-exists :supersede)
    (princ (animal-state->r animal-state)
           my-stream)))

では実際にシミュレーションを走らせて、結果をファイルに書き出します。

> (load "evolution")
T
> (write-animal-state "animal-state.R" (dump-animal-state 200000 10000))

20万日分シミュレーションを走らせ、1万日間隔で animal-state.R というファイルに書き出しました。

ちなみにdump前に (load "evolution") を実行しているのは世界をリセットするためです。 オリジナルのコードに手を入れずに手っ取り早く世界をリセットする方法が他に思いつきませんでした...

グラフの描画

Rでグラフを描画します。

plot-evolution.R

stateFile <- "animal-state.R"
stateList <- eval(parse(stateFile))

outDir <- "Rplot"
png(filename=file.path(outDir, "%03d.png"), width=400, height=300)

numTable <- data.frame(t(sapply(stateList,
                                function(x) {
                                    c(day=x[["day"]],
                                      num=length(x[["straightRate"]]))
                                })))
plot(numTable, type="o", pch=20, ylim=c(0, max(numTable[, 2])),
     main="動物の個体数", xlab="経過日数", ylab="動物の個体数")

maxX <- 1.0
maxY <- 160

for (state in stateList) {
    hist(state[["straightRate"]], breaks=seq(0, maxX, 0.02), ylim=c(0, maxY),
         main=sprintf("%7d 日後", state[["day"]]), col="gray",
         xlab="直進率", ylab="動物の個体数")
}

dev.off()

このRスクリプトを実行します。

$ mkdir Rplot
$ Rscript --vanilla plot-evolution.R

Rplot/ 以下に連番のpngファイルが生成されます。

結果

まずは動物の個体数の遷移です。

個体数は160前後で安定しています。これは植物が生える速度が常に一定のため、世界が平衡状態を保っているということでしょうか。

次に動物の直進率と個体数のヒストグラムを見てみます。

1万日後。まだ同じ傾向の動物しかいません。

3万日後。目立った変化無し。

4万日後。種の分化が始まった!

5万日後。完全に分かれました。


14万日後。ロバの直進率がさらに上がっています。


おわりに

  • 種の分化の様子を見ることができて楽しかったです。
  • Common LispとRの連携はやっつけだったので他の手法も試してみたいです。
  • 本の写経以外でCommon Lispのコードを書いたのは初めてなので、おかしなところがあったらどんどん指摘してください。

(追記)

他の手法を試してみました: Common LispでR (RCLを使ってみる)

2013年6月5日水曜日

Land of Lisp

Common Lisp始めました。

いまLand of Lispを読んでいます。
Common Lispでゲームを作って遊べる!楽しい!

読み終わったらまた感想書きます。

2013年5月17日金曜日

((Rで) 書く (Lisp) インタプリタ)

Rに慣れるために、Lispインタプリタを書いてみました。

元ねたはPeter Norvigの (How to Write a (Lisp) Interpreter (in Python)) (日本語訳: ((Pythonで) 書く (Lisp) インタプリタ))です。

コード

遊び方

コードを取ってきます。

$ git clone https://gist.github.com/5598108.git

Rインタプリタを起動してコードを読み込ませます。
repl() を実行するとLispの対話式インタプリタが起動します。

$ cd 5598108/
$ R -q
> source("lisp.R")
> repl()
lisp.R> (+ 1 2)
3
lisp.R> (define l (list 1 2 3))
(1 2 3)
lisp.R> (car l)
1
lisp.R> (cdr l)
(2 3)
lisp.R> (define add2 (lambda (x) (+ x 2)))
#<closure>
lisp.R> (add2 40)
42
lisp.R> (equal? (list 1 2) (list 1 (- 5 3)))
TRUE
lisp.R> (define map (lambda (f l) (if (null? l) (quote ()) (cons (f (car l)) (map f (cdr l))))))
#<closure>
lisp.R> (map add2 (list 10 20 30))
(12 22 32)
lisp.R>

感想

Rの言語機能(ファーストクラスの関数、レキシカルスコープ、クロージャなど)のおかげで、割と楽に書けたと思います。

ただリストの操作はちょっと面倒くさいかな、と思いました。
あと値渡しとか遅延評価とかでところどころはまりました。まだまだ慣れが必要です。

でも今回初めて使った機能がいくつかあって勉強になったし、何より言語処理系を作るのはとても楽しいです。

(追記)

続編を書きました。((Rで) 書く ((もっとRっぽい) Lisp) インタプリタ)

2013年5月10日金曜日

Rとクロージャ その2

前回の続きです。

Rのプロンプトで demo(scoping) を実行すると、クロージャの使用例を見ることができます。
SICPの銀行口座 (3.1.1 Local State Variables) をRで実装したものだと思います。

> demo(scoping)


        demo(scoping)
        ---- ~~~~~~~

Type  <Return>   to start : 

> ## Here is a little example which shows a fundamental difference between
> ## R and S.  It is a little example from Abelson and Sussman which models
> ## the way in which bank accounts work.       It shows how R functions can
> ## encapsulate state information.
> ##
> ## When invoked, "open.account" defines and returns three functions
> ## in a list.  Because the variable "total" exists in the environment
> ## where these functions are defined they have access to its value.
> ## This is even true when "open.account" has returned.  The only way
> ## to access the value of "total" is through the accessor functions
> ## withdraw, deposit and balance.  Separate accounts maintain their
> ## own balances.
> ##
> ## This is a very nifty way of creating "closures" and a little thought
> ## will show you that there are many ways of using this in statistics.
> 
> open.account <- function(total) {
+ 
+     list(
+        deposit = function(amount) {
+            if(amount <= 0)
+                stop("Deposits must be positive!\n")
+            total <<- total + amount
+            cat(amount,"deposited. Your balance is", total, "\n\n")
+        },
+        withdraw = function(amount) {
+            if(amount > total)
+                stop("You don't have that much money!\n")
+            total <<- total - amount
+            cat(amount,"withdrawn.  Your balance is", total, "\n\n")
+        },
+        balance = function() {
+            cat("Your balance is", total, "\n\n")
+        }
+        )
+ }

> ross <- open.account(100)

> robert <- open.account(200)

> ross$withdraw(30)
30 withdrawn.  Your balance is 70 


> ross$balance()
Your balance is 70 


> robert$balance()
Your balance is 200 


> ross$deposit(50)
50 deposited. Your balance is 120 


> ross$balance()
Your balance is 120 


> try(ross$withdraw(500)) # no way..
Error in ross$withdraw(500) : You don't have that much money!

> 

銀行口座の状態(total: 残高)はカプセル化されていて、外部からは直接参照できません。 代わりにアクセサ(deposit: 預入れ、 withdraw: 払出し、 balance: 残高確認) を通してやり取りしています。

オブジェクト指向っぽいものを簡単に実現できていて素敵ですね。

2013年4月29日月曜日

Rとクロージャ

Rでクロージャを使ってみます。

例:アキュムレータを生成する関数

試しにPaul Grahamのエッセイ、Revenge of the Nerds に出てくる「アキュムレータを生成する関数」を書いてみます。

Schemeならこんな感じ。

gosh> (define (make-accumulator n)
        (lambda (i)
          (set! n (+ n i))
          n))
make-accumulator
gosh> (define A (make-accumulator 5))
A
gosh> (A 10)
15
gosh> (A 10)
25

Perl5だと

$ perl -de0

Loading DB routines from perl5db.pl version 1.33
Editor support available.

Enter h or `h h' for help, or `man perldebug' for more help.

main::(-e:1):   0
  DB<1> sub make_accumulator {  \
  cont:   my $n = shift;        \
  cont:   sub { $n += shift }   \
  cont: }

  DB<2> $a = make_accumulator 5

  DB<3> p $a->(10)
15
  DB<4> p $a->(10)
25

(他の言語での実装はこちらを参照。)

これをRで書くと

> makeAccumulator <- function(n) {
+   function(i) {
+     n <<- n + i
+     n
+   }
+ }
> a <- makeAccumulator(5)
> a(10)
[1] 15
> a(10)
[1] 25

短く書くとこんな感じでしょうか。

makeAccumulator <- function(n) function(i) (n <<- n + i)

スーパーアサインメント演算子 <<-

ポイントはスーパーアサインメント演算子 <<- です。

通常の代入演算子 <- は同じ環境の変数に代入します。

<<- は入れ子になった環境を順に外側に検索し、最初に変数の定義が見つかった環境に代入します。
見つからなかった場合はグローバル環境に代入します。
詳細は help(<<-) を実行してみてください。

今回の例では、内側の関数から、外側の関数の環境の n に代入したいので、 <- ではなく <<- を使う必要があります。

2013年4月14日日曜日

Ubuntu 12.04にRをインストール

Ubuntu 12.04にRをインストールします。
$ sudo apt-get install r-base-core

Emacsから使いたい場合はESSも入れると良いです。
$ sudo apt-get install ess
Emacsを起動して M-x R で対話型インタプリタが立ち上がります。