2013年12月11日水曜日

【PostgreSQL】PL/Pythonで機械学習


PostgreSQL Advent Calendar 2013 12/11 です。

DBに格納したデータに対して機械学習のアルゴリズムを適用したくなったとします。
DBからデータを出力して他のプログラム言語で処理してもいいのですが、PostgreSQLで処理が完結した方が便利なこともあるでしょう。

しかし、PostgreSQLで処理するために機械学習のアルゴリズムを一から自分で実装するのは大変です。

PL/Pythonを使えば、pythonで書かれたライブラリを利用して手軽に機械学習を行うことができます。
また、pythonの機械学習ライブラリはC言語で書かれていることが多いので、速度面でも悪くないです。

インストール

PL/Python

PostgreSQLをソースコードからインストールしている場合は、configureのオプションに--with-pythonをつけてビルドします。
今回はpython2系を使うので、ビルド時にpython2系にパスが通っているか、PYTHON環境変数にpythonのパスを指定するかして下さい。

OSのパッケージ管理でインストールしている場合は、apt-get install postgresql-plpython-9.xなどを試してください。

その後、CREATE EXTENSIONすればPL/Pythonが使えるようになります。
CREATE EXTENSION plpython2u;

scipy

pythonのscipyライブラリをインストールします。
pip install scipy
上記ですんなりインストールできればいいですが、失敗した場合は結構苦労するかもしれません。

実装

k-means法で2次元ベクトルのクラスタリングをやってみます。
scipyの使い方は下記を参考にして下さい。

K-means clustering and vector quantization (scipy.cluster.vq)

kmeans関数の作成

複数のベクトル(x, y)から、xとyそれぞれを配列にして渡し、
各ベクトルがどのクラスタに属すかを配列で取得する関数を作成します。

CREATE FUNCTION kmeans(x float[], y float[]) RETURNS int[] AS $$
 from numpy import array
 from scipy.cluster.vq import vq, kmeans, whiten

 features = array(zip(x, y))
 whitened = whiten(features)
 book = array((whitened[0], whitened[2]))
 codebook, distortion = kmeans(whitened, book)
 code, dist = vq(whitened, codebook)

 return list(code)
$$ LANGUAGE plpython2u;

SELECT kmeans(
 array[1.9,1.5,0.8,0.4,0.1,0.2,2.0,0.3,1.0],
 array[2.3,2.5,0.6,1.8,0.1,1.8,0.5,1.5,1.0])

/*
       kmeans        
---------------------
 {0,0,1,1,1,1,1,1,1}
(1 row)

2番目のベクトル(1.5,2.5)はクラスタ0に属する
*/

クラスタリングする

ベクトルの集合をkmeans関数に渡せるように配列に集約し、kmeansで得られた配列のn番目が集合のn行目のクラスタの値になります。

SELECT
 x, y,
 (
  kmeans(
   array_agg(x) OVER(),
   array_agg(y) OVER()
  )
 )[row_number() OVER()] as kmeans_cluster
FROM (VALUES
 (1.9::float,2.3::float),
 (1.5,2.5),
 (0.8,0.6),
 (0.4,1.8),
 (0.1,0.1),
 (0.2,1.8),
 (2.0,0.5),
 (0.3,1.5),
 (1.0,1.0)
) t(x, y)

/*
  x  |  y  | kmeans_cluster 
-----+-----+----------------
 1.9 | 2.3 |              0
 1.5 | 2.5 |              0
 0.8 | 0.6 |              1
 0.4 | 1.8 |              1
 0.1 | 0.1 |              1
 0.2 | 1.8 |              1
   2 | 0.5 |              1
 0.3 | 1.5 |              1
   1 |   1 |              1
(9 rows)
*/

以上、PL/Pythonでscipyライブラリを使ってk-means法を実装してみました。

おまけ

PostgreSQLのC言語関数で書かれたkmeansのライブラリもあります。
https://github.com/umitanuki/kmeans-postgresql

2012年12月14日金曜日

【PostgreSQL】SQLでライフゲームの実装


PostgreSQL Advent Calendar 2012 12/14 です。

ライフゲームとは「生命の誕生、進化、淘汰などのプロセスを簡易的なモデルで再現したシミュレーションゲーム」のことです。
碁盤目状に配置されたセルが、隣接するセルとの関係により、誕生・生存・死滅していく様子をシミュレーションします。

ルールは下記参考
ライフゲーム - Wikipedia

このライフゲームをSQLを使って実装してみます。

データ構造

セルの集合を下記のように表現します。
・1つのセルにつき、1行のレコード
・セルは3つの属性(列)を持つ
 integer x:碁盤の左端から右に何個目か
 integer y:碁盤の上端から下に何個目か
 integer alive:生存状態(生きているセルなら1、死んでいるセルなら0)

データはこれだけです。


初期化

碁盤のサイズ、生存しているセルの割合、を決めてランダムに生成
CREATE TABLE cells AS
SELECT
 x::integer, y::integer
 ,CASE WHEN random() < 0.3 THEN 1 -- 3割のセルが生存
 ELSE 0
 END::integer AS alive
FROM
 generate_series(1,5) t1(x) -- 横のサイズが5
 ,generate_series(1,5) t2(y) -- 縦のサイズが5
;
/*
SELECT * FROM cells;
 x | y | alive 
---+---+-------
 1 | 1 |     1
 1 | 2 |     0
 1 | 3 |     0
 1 | 4 |     0
 1 | 5 |     0
 2 | 1 |     0
 2 | 2 |     0
 2 | 3 |     1
 2 | 4 |     1
 2 | 5 |     0
 3 | 1 |     1
 3 | 2 |     1
 3 | 3 |     1
 3 | 4 |     1
 3 | 5 |     0
 4 | 1 |     0
 4 | 2 |     0
 4 | 3 |     0
 4 | 4 |     0
 4 | 5 |     0
 5 | 1 |     0
 5 | 2 |     0
 5 | 3 |     0
 5 | 4 |     1
 5 | 5 |     0
(25 rows)
*/

描画

初期化したセルから碁盤を描画してみます。
aliveの値をテキストに集約して0,1の文字列で表現しています。
SELECT string_agg(alive_y, E'¥n' ORDER BY y) AS board
FROM(
 SELECT
  y
  ,string_agg(alive::text, '' ORDER BY x) AS alive_y
 FROM cells
 GROUP BY y
)t
;
/*
               board               
-----------------------------------
 10100¥n00100¥n01100¥n01101¥n00000
(1 row)

10100
00100
01100
01101
00000
*/

次世代の生成

まずはSQLを
SELECT
 x, y
 ,(
  SELECT
   CASE WHEN c1.alive = 0 AND c2.alives = 3 THEN 1 -- 誕生
   WHEN c1.alive = 1 AND (c2.alives = 2 OR c2.alives = 3) THEN 1 -- 生存
   ELSE 0 -- 死滅
   END
  FROM(
   SELECT sum(c2.alive) as alives
   FROM cells c2
   WHERE
    NOT(c1.x = c2.x AND c1.y = c2.y) -- 自分自身を除く
    AND
    c1.x BETWEEN c2.x - 1 AND c2.x + 1 -- 横の範囲 
    AND
    c1.y BETWEEN c2.y - 1 AND c2.y + 1 -- 縦の範囲
  )c2
 )::integer AS alive
FROM cells c1
;
/*
 x | y | alive 
---+---+-------
 1 | 1 |     0
 1 | 2 |     0
 1 | 3 |     0
 1 | 4 |     0
 1 | 5 |     0
 2 | 1 |     1
 2 | 2 |     0
 2 | 3 |     0
 2 | 4 |     1
 2 | 5 |     0
 3 | 1 |     0
 3 | 2 |     1
 3 | 3 |     0
 3 | 4 |     1
 3 | 5 |     0
 4 | 1 |     0
 4 | 2 |     1
 4 | 3 |     0
 4 | 4 |     1
 4 | 5 |     0
 5 | 1 |     0
 5 | 2 |     0
 5 | 3 |     0
 5 | 4 |     0
 5 | 5 |     0
(25 rows)


01000
00110
00000
01110
00000
*/

ちょっと解説

次世代のaliveを生成しているサブクエリを具体的なケースに置き換えてみます。
x=2, y=3のセルは次世代で生存しているでしょうか?
サブクエリをc1.x=2, c1.y=3, c1.alive=1で置換したものが下記になります。
  SELECT
   CASE WHEN 1 = 0 AND c2.alives = 3 THEN 1 -- 誕生
   WHEN 1 = 1 AND (c2.alives = 2 OR c2.alives = 3) THEN 1 -- 生存
   ELSE 0 -- 死滅
   END
  FROM(
   -- ここのサブクエリで周囲の生存セル数を計算
   SELECT sum(c2.alive) as alives
   FROM cells c2
   WHERE
    NOT(2 = c2.x AND 3 = c2.y) -- 自分自身を除く
    AND
    2 BETWEEN c2.x - 1 AND c2.x + 1 -- 横の範囲 
    AND
    3 BETWEEN c2.y - 1 AND c2.y + 1 -- 縦の範囲
  )c2
具体的なケースができたら、後は全てのセルに対して同様に処理してやればOKです。

状態を保存する場合はUPDATEするなり、CREATE TABLEするなりしてやりましょう。


まとめ

ライフゲーム?それ、SQLでもできるよ。


Advent Calendar

明日はsakamotomsさんです。



おまけ1

任意の碁盤の状態からcellsを初期化したい場合
CREATE OR REPLACE FUNCTION board_to_cells(board text)
 RETURNS TABLE(x integer, y integer, alive integer)
 AS $$
  SELECT
   row_number() OVER(PARTITION BY y)::integer as x
   ,y::integer
   ,alive::integer
  FROM(
   SELECT
    row_number() OVER() AS y
    ,regexp_split_to_table(alive_y, '') AS alive
   FROM(
    SELECT
     regexp_split_to_table($1, '\n') AS alive_y
   )t
  )t
 $$
LANGUAGE SQL;

/*
SELECT * FROM board_to_cells(E'010\n100\n100');
*/


おまけ2

パフォーマンス実験
手元の環境にて、初期状態が100×100, 生存率3割で試したところ、次世代の計算にかかる時間が20秒強ほど。
もっと高速な実装を考えてみるのも面白いかも。

2012年7月29日日曜日

【JavaScript】【CI】Jasmine + PhantomJS + JenkinsでJavascriptのテスト


Jasmine

JavaScript用のテスティングフレームワーク
使い方は下記等を参考

JasmineによるJavaScriptのテスト その1 « サーバーワークス エンジニアブログ
http://blog.serverworks.co.jp/tech/2010/11/30/jasmine-tutorial-1/

・インストール(というほどのものでもないが)

以下からzipをダウンロード
https://github.com/pivotal/jasmine/downloads

・サンプル

適当なディレクトリに展開してSpecRunner.htmlをブラウザで表示
例)apacheのDocumentRootにjasmineディレクトリを作成
http://localhost/jasmine/SpecRunner.html

テスト結果が表示されればOK


PhantomJS

Headless WebKit with JavaScript API
画面のないwebkitブラウザ

phantomjsを使うことでjsのテストをコマンド実行できる

・インストール

Windowsへのインストールはバイナリがダウンロードできる
http://code.google.com/p/phantomjs/downloads/list

Unix系なら下記参考
http://phantomjs.org/build.html

・サンプル

例)SpecRunner.htmlのテストを実行
SpecRunner.htmlのreporterを変更

実行
/path-to-phantomjs/bin/phantomjs /path-to-phantomjs/examples/run-jasmine.js http://localhost/jasmine/SpecRunner.html
'waitFor()' finished in 200ms.
Player


Jenkins

ソフトウェアのビルドやcronで起動するジョブなどの繰り返しジョブの実行を監視する。

Jenkins - 日本語 - Jenkins Wiki
https://wiki.jenkins-ci.org/display/JA/Jenkins

・インストール

java -jar jenkins.warを実行するか、サーブレットコンテナにデプロイ

例)tomcatにデプロイする場合
/path-to-tomcat/webapps/jenkins.warを配置
http://localhost:8080/jenkinsにアクセス


Phatomjs + Jasmineでテスト結果をJUnit形式のxmlに出力

以下からライブラリをダウンロード
https://github.com/detro/phantomjs-jasminexml-example

phantomjs_jasminexml_runner.jsを/path-to-phantomjs/testに配置
core.jsを/path-to-phantomjs/test/lib/utilsに配置

jasmine.phantomjs-reporter.jsをjasmineのテストhtmlで読み込み、reporterを変更

例)SpecRunner.htmlのテスト結果をxmlに出力
SpecRunner.htmlのreporterを変更


実行
/path-to-phantomjs/bin/phantomjs /path-to-phantomjs/test/phantomjs_jasminexml_runner.js http://localhost/jasmine/SpecRunner.html [path-to-output]
[path-to-output]で指定したディレクトリにTEST-Player.xmlが出力される。


Jenkinsでテストをジョブ登録

Jenkinsの画面から
①新規ジョブの作成を押す
②ジョブ名を入力し、フリースタイル・プロジェクトのビルドを選択してOKを押す
③ビルド→ビルド手順の追加→シェルの実行でテスト実行コマンドを登録

例)SpecRunner.htmlの場合は下記を指定
/path-to-phantomjs/bin/phantomjs /path-to-phantomjs/test/phantomjs_jasminexml_runner.js http://localhost/jasmine/SpecRunner.html $WORKSPACE
$WORKSPACEはjenkinsの環境変数。通常は$JENKINS_HOME/workspace/[job-name]
tomcatにデプロイした場合、$JENKINS_HOMEは$TOMCAT_USER_HOME/.jenkinsがデフォルト

④ビルド後の処理→JUnitテスト結果の集計を選択してテスト結果XMLを指定
  $WORKSPACEに出力している場合は*.xml
⑤設定を保存し、ビルドを実行して動作確認。ビルド履歴からテスト結果等が閲覧できる。
⑥その他、ビルド・トリガで自動実行、ビルド後の処理でメール通知、等を設定

2012年3月11日日曜日

PostgreSQLで統計解析 PL/R 基礎編3 : CRAN(Rライブラリ)利用


CRAN

perlにCPANがあるように、CRANというRのモジュールを公開するネットワークがあり、多くのプログラムが公開されている。
The Comprehensive R Archive Network

下記記事を参考に、PL/RでCRANから取得したモジュールを利用してみる。
zooパッケージを使って行列の欠損値を補間する(R Advent Calendar 2011) - My Life as a Mock Quant


Install R Package

パッケージのinstallはRから
r #R起動
> install.packages("zoo")


Use R Package

パッケージはlibrary(zoo)で読み込み。
行列の欠損を補完する関数を作成
CREATE OR REPLACE FUNCTION na_locf(matrix int4[][]) RETURNS int4[][] AS '
    library(zoo);
    return(na.locf(matrix))
' LANGUAGE 'plr';
SELECT na_locf(ARRAY[ARRAY[1,2,3],ARRAY[4,null,6],ARRAY[7,8,9]]);
-- nullが2で補完される。

CREATE OR REPLACE FUNCTION na_approx(matrix int4[][]) returns int4[][] as ' 
    library(zoo);
    return(0.5 * (na.approx(matrix) + t(na.approx(t(matrix)))))
' LANGUAGE 'plr';
SELECT na_approx(ARRAY[ARRAY[1,2,3],ARRAY[4,null,6],ARRAY[7,8,9]]);
--nullが5で補完される。
PostgreSQLとRのデータ構造の違いに気をつければライブラリの使用は比較的簡単

PostgreSQLで統計解析 PL/R 基礎編2 : R用データ読み込み


前回 に続き、
今回はR用のサンプルデータを読み込んでPostgreSQLのデータとして受け取ってみる。


Test Data

データ集めは下記を参考に
統計を学びたい人へ贈る、統計解析に使えるデータセットまとめ - ほくそ笑む

今回は、世界の48都市の経済状態を記録したデータのcsvを読み込んでみる。


Composite Type Sample

読み込むcsvは文字列と数値が含まれており、Rではデータフレームという形式で読み込まれる。
データフレーム形式のデータはPostgreSQLでは複合型として受け取れる。
PL/Rのマニュアル記載のサンプルは下記の通り。
CREATE TYPE emp AS (
    name text
    ,age int2
    ,salary int4
);
CREATE OR REPLACE FUNCTION get_emps() RETURNS SETOF emp AS '
   names <- c("Joe","Jim","Jon")
   ages <- c(41,25,35)
   salaries <- c(250000,120000,50000)
   df <- data.frame(name = names, age = ages, salary = salaries)
   return(df)
' LANGUAGE 'plr';

SELECT * FROM get_emps();
/*
name | age | salary
text | smallint | integer
------+-----+-----------
Joe | 41 | 250000
Jim | 25 | 120000
Jon | 35 | 50000
*/
複合型を作成するのが手間なので、下記のようにしてみたらうまくいった。
CREATE OR REPLACE FUNCTION get_emps(out name text, out age int4, out salary int4) RETURNS SETOF record AS '
   names <- c("Joe","Jim","Jon")
   ages <- c(41,25,35)
   salaries <- c(250000,120000,50000)
   df <- data.frame(name = names, age = ages, salary = salaries)
   return(df)
' LANGUAGE 'plr';
SELECT * from get_emps();
が、引数を与えた場合はうまく動作しない模様…。
マニュアルにもデータフレームをsetof recordで返せるとは明示されていないので、対応範囲外なのかもしれない。


Get CityEcon

では、本題のcityecon.csvを読み込んでみる。
CREATE TYPE cityecon AS (
    city text
    ,work int4
    ,price float8
    ,salary float8
);
CREATE OR REPLACE FUNCTION get_cityecon(csvpath text) RETURNS setof cityecon AS '
    x <- read.csv(csvpath)
    return(x)
' LANGUAGE 'plr';
SELECT * FROM get_cityecon('/path-to-csv/cityecon.csv');
-- ERROR:  invalid input syntax for integer: "-"
元データにnull値として"-"が入っているためエラー。
Rで"-"をNAに置き換えておく。
CREATE OR REPLACE FUNCTION get_cityecon(csvpath text) RETURNS setof cityecon AS '
    x <- read.csv(csvpath)
    x$Work <- ifelse(x$Work=="-", NA, x$Work) 
    x$Salary <- ifelse(x$Salary=="-", NA, x$Salary) 
    return(x)
' LANGUAGE 'plr';
SELECT * FROM get_cityecon('/path-to-csv/cityecon.csv');
-- OK
さらに、CSV以外のファイル形式にも対応できるようにしてみる。
CREATE OR REPLACE FUNCTION get_cityecon(path text, header boolean, separator text) RETURNS setof cityecon AS '
    if(is.null(separator)){
        sep <- ""
    } else sep <- separator    
    if(header){
        x <- read.table(path, header = T, sep = sep)
    }else{
        x <- read.table(path, header = F, sep = sep)
    }
    x$Work <- ifelse(x$Work=="-", NA, x$Work) 
    x$Salary <- ifelse(x$Salary=="-", NA, x$Salary) 
    return(x)
' LANGUAGE 'plr';
SELECT * FROM get_cityecon('/path-to-csv/cityecon.csv', true, ',')
-- OK
データの区切り文字とヘッダーが含まれるかどうかを指定可能。

読み込むデータに合わせて複合型を作成したり、NULL値の変換をしたりと結構めんどくさい。
もうちょっと簡単に読み込みたいなあという感想。

2011年12月6日火曜日

PostgreSQLで統計解析 PL/R 基礎編 : PostgreSQL Advent Calendar #6


PostgreSQL Advent Calendar 12/6 です。

ネタを決めるにあたり、analyticsをテーマに何かないものかと探してみたら、
PL/R(R手続き型言語)なるものがあるそうで。

Rは統計解析分野で人気のプログラミング言語。
オープンソースで、コミュニティの活動も活発なようです。
参考:R Advent Calendar 2011 : ATND

では、PL/Rを導入して、ちょこっと触ってみます。


Install
今回構築した環境は下記。
・Mac OS X 10.6 Snow Leopard
・gcc 4.2.1
・postgresql 9.0.4
・R 2.14.0
・plr-8.3.0.13

Windows,Linuxへのインストールは下記をご参考
R と PostgreSQL - RjpWiki

まずRをインストール。
(macのhomebrewを使用)
brew install r
Error: This formula requires a fortran compiler
失敗。fortranコンパイラをインストールして再チャレンジ。
brew install gfortran
brew install r
OK。インストール時に一般的なインストール先に
リンクさせた方がいいよとメッセージが出るので、
sudo ln -s "/usr/local/Cellar/r/2.14.0/R.framework" /Library/Frameworks
これでRのインストールは完了。
なお、パッケージインストールも可能。
参考:R for Mac OS X

次にPL/Rをインストール。
ソースは→ http://joeconway.com/plr/ から入手。
su - postgres
cd /path-to-postgres_source/contrib #ここにPL/Rのソースを配置
tar zxf plr-8.3.0.13.tar.gz
cd plr
export R_HOME=/Library/Frameworks/R.framework/Resources
make
make install
PL/Rのインストール完了。テスト用のDBを作成して言語を登録。
createdb r_test
cd /path-to-pgsql/contrib
psql r_test < plr.sql
psql -c "SELECT * FROM pg_language WHERE lanName = 'plr'" r_test #確認
これで導入完了。


Try working

PL/Rのドキュメントにある集約関数をテストしてみます。
(中央値を返す集約関数medianを作成)
create or replace function r_median(_float8) returns float as '
  median(arg1)
' language 'plr';
CREATE AGGREGATE median (
  sfunc = plr_array_accum,
  basetype = float8,
  stype = _float8,
  finalfunc = r_median
);

SELECT median(v)
FROM (values(1.0),(2.5),(5.7),(3.4),(-0.9)) AS t(v);
-- => select 2.5
集約関数はfloat8を配列に集約して最後にr_median(_float8)を呼び出し、
Rの関数medianを計算しています。

ポイントはデータ型の違い。
R言語はベクトル処理言語であり、データはベクトルで持ちます。
例えば、上記のmedianの計算は、Rでは下記のように行います。
x <- c(1.0, 2.5, 5.7, 3.4, -0.9) #ベクトルを作成してxに代入
median(x) # => [1] 2.5
postgresからRへの引数の受け渡しについては、
・スカラ ⇒ 単一要素のベクトル(例外あり)
・ 一次元配列 ⇒ 複数要素のベクトル
・二次元配列 ⇒ 行列
・三次元配列 ⇒ 三次元配列
・三次元以上の配列 ⇒ サポートせず
・ 複合型 ⇒ データフレーム
となります。
参考:Passing Data Values

では最後に、統計解析の基礎として相関係数を計算する関数を作成してみます。
参考:R-Source
CREATE OR REPLACE FUNCTION correlation(_float8, _float8) RETURNS float8 AS '
 cor(arg1, arg2)
' LANGUAGE 'plr';
列xと列yの相関係数を計算する場合は下記で。
SELECT correlation(array_agg(x),array_agg(y))
FROM (values(0.7,1.9),(-1.6,0.8),(-0.2,1.1),(-1.2,0.1),(-0.1,-0.1)
      ,(3.4,4.4),(3.7,5.5),(0.8,1.6),(0.0,4.6),(2.0,3.4)) AS t(x,y);
-- => select 0.795102...

以上、簡単にさわってみました。
それほど敷居は高くないように感じますが、どういう用途に使ったものか...


Be continued
応用編に続く ...のか?


Next
明日はs87さんです。よろしくお願いします。

2011年9月30日金曜日

【PostgreSQL】2つの配列に共通する要素のインデックスを取得する関数

やりたいことは、配列A{1,3,4,5,8,10,11,12}と配列B{1,2,5,11,13}から共通する要素(1,5,11)を探し、
配列Aの要素番号(1,4,7)を取得すること。(配列Aの要素はユニークである仮定)


配列Bが数値1つでよければ、contribにidx関数がある。

idx(int[], int item)
戻り値の型:int
説明:itemに一致する要素番号(存在しなければ0)
例:idx(array[11,22,33,22,11], 22)
結果:2

intarray
http://www.postgresql.jp/document/current/html/intarray.html


これを利用して関数を作成(sort()とuniq()もintarrayモジュールに含まれる)

CREATE OR REPLACE FUNCTION idx_array(_int4, _int4) RETURNS _int4 AS
'SELECT
    uniq(sort(array_agg(idx($1, t.v))))
FROM
    unnest($2) AS t(v)
WHERE
    idx($1, t.v) > 0'
LANGUAGE SQL IMMUTABLE STRICT;

SELECT idx_array(ARRAY[1,3,4,5,8,10,11,12], ARRAY[1,2,5,11,13])
結果:{1,4,7}