(10/25) Mishima.syk #4やります

connpassを用意したので参加登録お願いします。

発表者も絶賛募集中ですので@fmkz___まで。

僕はCytoscape+chemviz2の話かiPython notebookの話でもしようかなと考えています。

この本良さげ。多分買う

ProductName IPython Interactive Computing and Visualization Cookbook
Cyrille Rossant
Packt Publishing / ?円 ( 2014-09-25 )


実践コンピュータビジョンの8章をScikit-learnで

kNN, Naibe Bayesm, SVM, (Random Forrest)をScikit-learnでやってみた。データはiPython NotebookでReSTで出力したものをpandocでmarkdown_strictに変換しなおしてblogに貼り付けた。

描画用のヘルパー関数とデータセットの生成

from matplotlib.colors import ListedColormap
import Image
from numpy import *
from pylab import *
import pickle

def myplot_2D_boundary(X,y):
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02),
                        np.arange(y_min, y_max, 0.02))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

    Z = Z.reshape(xx.shape)
    plt.figure()
    plt.pcolormesh(xx, yy, Z, cmap=cmap_light)

    plt.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())

    plt.show()

with open('points_normal.pkl', 'r') as f:
  class_1 = pickle.load(f)
  class_2 = pickle.load(f)
  labels = pickle.load(f)

X_normal = np.r_[class_1, class_2]
y_normal = labels

with open('points_ring.pkl', 'r') as f:
  class_1 = pickle.load(f)
  class_2 = pickle.load(f)
  labels = pickle.load(f)

X_ring = np.r_[class_1, class_2]
y_ring = labels

with open('points_normal_test.pkl', 'r') as f:
  class_1 = pickle.load(f)
  class_2 = pickle.load(f)
  labels = pickle.load(f)

X_normal_test = np.r_[class_1, class_2]
y_normal_test = labels

with open('points_ring_test.pkl', 'r') as f:
  class_1 = pickle.load(f)
  class_2 = pickle.load(f)
  labels = pickle.load(f)

X_ring_test = np.r_[class_1, class_2]
y_ring_test = labels

cmap_light = ListedColormap(['#FFAAAA', '#AAFFAA'])
cmap_bold = ListedColormap(['#FF0000', '#00FF00'])

kNN

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn import neighbors, datasets

clf = neighbors.KNeighborsClassifier(3)
clf.fit(X_normal, y_normal)

myplot_2D_boundary(X_normal,y_normal)

image

clf = neighbors.KNeighborsClassifier(3)
clf.fit(X_ring, y_ring)

myplot_2D_boundary(X_ring,y_ring)

image

ベイズ

from sklearn.naive_bayes import GaussianNB

clf = GaussianNB()
clf.fit(X_normal, y_normal)
labels_pred = clf.predict(X_normal_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_normal_test.shape[0],(y_normal_test != labels_pred).sum())

myplot_2D_boundary(X_normal,y_normal)

image

clf = GaussianNB()
clf.fit(X_ring, y_ring)
labels_pred = clf.predict(X_ring_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_ring_test.shape[0],(y_ring_test != labels_pred).sum())

myplot_2D_boundary(X_ring,y_ring)

image

SVM

from sklearn import svm
clf = svm.SVC()
clf.fit(X_normal, y_normal)

labels_pred = clf.predict(X_normal_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_normal_test.shape[0],(y_normal_test != labels_pred).sum())

myplot_2D_boundary(X_normal, y_normal)

image

clf = svm.SVC()
clf.fit(X_ring, y_ring)
labels_pred = clf.predict(X_ring_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_ring_test.shape[0],(y_ring_test != labels_pred).sum())

myplot_2D_boundary(X_ring,y_ring)

image

Random Forest

from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier(n_estimators=10)

clf.fit(X_normal, y_normal)
labels_pred = clf.predict(X_normal_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_normal_test.shape[0],(y_normal_test != labels_pred).sum())

myplot_2D_boundary(X_normal,y_normal)

image

from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier(n_estimators=10)

clf.fit(X_ring, y_ring)
labels_pred = clf.predict(X_ring_test)

print "Number of mislabeled points out of a total %d points : %d" % (y_ring_test.shape[0],(y_ring_test != labels_pred).sum())

myplot_2D_boundary(X_ring,y_ring)

image

実践コンピュータビジョンの8章に参加してきました

最近ちょっと忙しくて実践コンピュータビジョンの読書会には初参加なのに発表してきたわけだが。

写経以外に務めたことはipython notebookとscikit-learnを推してきたw。あとディープラーニングの話とかしてた。そして少しまじめにディープラーニングを学ぼうと思った。

懇親会は筋肉居酒屋で。

1411363268

店に入ると鉄アレイ等がお出迎え

1411363265

塩バターラーメン風パスタ(一人前のハーフサイズw)

1411363267

久しぶりに参加して楽しかったですね。主催者がいい感じにバトンタッチしつつ、新しい人も程よく入りながら長いこと続いているいい読書会だなぁと思いました。数えてみたらもうちょっとで5年ですね。

次回は島田でやるそうです。

RHEL6にサイエンス系のPythonパッケージを導入する

最近RHEL6を与えられたのだが、Pythonのバージョンが2.6系なので2.7系を入れつつ以下のパッケージを導入したのでメモ

  • numpy
  • scipy
  • scikit-learn
  • matplotlib
  • ggplot
  • ipython (ipython notebook)

Pythonのインストール

開発環境をrpmで入れておく

yum groupinstall "Development tools"
yum install zlib-devel bzip2-devel openssl-devel ncurses-devel sqlite-devel readline-devel tk-devel gdbm-devel db4-devel libpcap-devel xz-devel

あとはPythonをソースからインストールする

pipのインストール

wget https://bootstrap.pypa.io/get-pip.py
python get-pip.py

blas,lapackのインストール

scipyにはlapack(付きのnumpy)が必要なのだけどyum install blas,lapackは上手くいかないのでソースからインストールした。

wget http://www.netlib.org/lapack/lapack.tgz
tar xzfv lapack.tgz
cd lapack-*/
cp INSTALL/make.inc.gfortran make.inc

meke.incのオプションを修正する -fPICオプションを追加。もし64ビットマシンなら-m64オプションも追加

FORTRAN  = gfortran
OPTS     = -O2 -frecursive -fPIC -m64
DRVOPTS  = $(OPTS)
NOOPT    = -O0 -frecursive -fPIC -m64

書き換えたらmakeする

make blaslib; make lapacklib

出来た*.aを適当なディレクトリに配置して環境変数を設定し、.bashrcとか/etc/profileに追加しておく

export BLAS=/[path]/[to]/librefblas.a
export LAPACK=/[path]/[to]/liblapack.a

numpy,scipyのインストール

pip install numpy

インストール出来たらblas,lapackが使われているかどうかを確認するためimport numpyしてnumpy.show_config()で確認しておく。

OKだったらscipyを入れる。

pip install scipy

matplotlibのインストール

libpngが必要なのでyumで入れる。それからRHEL6のfreetypeは2.3だがmatplotlib1.4.0でも動くので設定ファイルを書き換えてコンパイルする。

yum install libpng libpng-devel

1.4.0のソースをダウンロード

tar xvfz matplotlib-1.4.0.tar.gz

setupext.pyでfreetypeの2.4以上をチェックしているところを2.3に書き換える

python setup.py install

scikit-learn, pandas, ggplot, ipythonのインストール

入れるだけ

pip install scikit-learn
pip install pandas
pip install patsy
pip install ggplot
pip install pyzmq
pip install jinja2
pip install tornado
pip install ipython

PyconJP2014に参加しました

PyconJP2014に参加しました。運営のヒト、発表者のヒト、参加者の皆様、楽しい時間をありがとうございました。

去年のPyconで発表して燃え尽きた感があったんだけど、今年参加したら再生した感があったというか、コードを書く気力が戻ってきたかなw

今回は科学技術計算系のセッションが多くて非常に楽しめた。

参加したセッション(Keynoteを除く)

Deep Learning for Image Recognition in Python (ja)から。

Decaf使ってみようと思ったらDEPRECATEDになってた。後継のcaffeを使えってことらしい

1410783994

リファクタリングツールあれこれ (May the force be with you) (ja)から。

データベーススキーマからSQLAlchemyのモデルを自動生成してくれるsqlacodegenが刺さった。コレ使えばpychembldbもっと簡単にかけたんじゃないの?と。

1410783999

あとは数独を解くソルバの話がとてもおもしろかった。

尚、今週末は実践コンピュータビジョン読書会の画像検索の章なので、参加するといいと思います。

資料作らないトナー

日本酒ラベルの特徴点抽出

今月出番らしいので実践コンピュータビジョンを読み返している。ここ一年くらい興味が読書とか将棋に移ってしまいプログラミングへの関心がちょっと薄れてたけど、写経しだすとやっぱプログラミングは面白いですね☆

ProductName 実践 コンピュータビジョン
Jan Erik Solem
オライリージャパン / 3150円 ( 2013-03-23 )


ちなみに8章はサポートベクターマシンとかkNNを使った画像認識なので2章と7章を事前に読んでおけばいい感じです。

それから久々の富士開催なので、これを機会につけナポリタンとか杉山フルーツの生ゼリーとか食べるといいかと思います。

富士がんもいっちもあったわ。

8章では画像の特徴量としてSIFT記述子を使うので、それを使って最近飲んだ日本酒の画像の特徴抽出をしてみた。

萩の鶴 夕涼み

柔らかい酸と甘さが特徴の萩の鶴

1409619062

風鈴と猫が特徴が出ており、夕涼み感が感じられる

hagi

庭のうぐいす おうから

庭のうぐいすは、結構好きなんだけどおうからはちょっときりっとしすぎな感じでいまいちだった。 燗つければよかったのかなぁ。

1409619064

庭のうぐいすらしくなかったので、うぐいすが特徴としてあまり出てこないのも納得。さすがSIFT

uguisu

澤屋まつもと 守破離

酸が強めのお酒が飲みたくて購入

1409619066

円の大きさに酸の強さが現れているはずw

matsumoto

一白水成 特別純米酒

安定した美味しさ

1409619067

サークルの数が美味さを表している☆

suisei

秋鹿

秋鹿、美味いよねー 序盤、中盤、終盤、隙がないと思うよ。だけど…俺は、負けないよ

1409619069

というような決意が特徴として現れていますね

akishika

結論

日本酒のラベルの画像から特徴抽出することで、日本酒の美味さをある程度反映することができる事がわかりました。 興味をもたれたら静岡Developers勉強会 コンピュータビジョン vol.8に参加するといいと思います。

xlsの日付はfloat型

最近はxlsでダウンロードしてきたデータをxlrdで読みだしてSQLAlchemy経由でoracleにimportするようなものを作っている。

Excelで日付が文字列として入っている場合は strptime, fromtimestamp, mktimeを組み合わせてOracleに突っ込めばいい。

普通にExcelのDate型の場合、実態はfloatらしいので読み出すときにちょっと困るが、Date型ではいってる場合はxldate_as_tupleっていうxlrdのユーティリティ関数を使って

publication_date = datetime(*xldate_as_tuple(sheet.cell(r, 21).value, wb.datemode))

みたいにできる。

でもExcelなのでというか、parse出来ないよエラー(floatの値が何故か小さすぎる?)とかValueErrorとか普通に出るというか、元データがダーティ過ぎてそっちの対応に四苦八苦してたw

ggplotっぽい配色

ggplotっぽい配色が欲しかったので探したらsofで見つけたがこれはRのコードなのでPythonで

colorsysというモジュールを使う。

でもこれは0-1の範囲の数字がタプルで返ってきて使いにくいので255をかけた後16進数表記にして返す関数を書いた。

from colorsys import hls_to_rgb

def color_list(n, l=0.65, s=1.0):
    return ['#%02x%02x%02x' % tuple(map(lambda x: int(x * 255), \
        hls_to_rgb((i + 3 - n)/ float(n), l, s))) for i in range(n)]

10個ほど生成させてみるといい感じ

Bioinformatics版のProject Euler

RosalindというProject Eulerっぽいサイトがあったので少し遊んでみた。

楽しくbioinformaticsの問題を解きながらアルゴリズムを覚えていく感じなので、BioPythonを覚えるのにちょうどいいかなと思った。

GC Contentを計算するメソッドくらいあるだろうとは思うんだけど知らないので素朴に計算。リスト内包表記でmaxで取り出すのがいいかなと思ったけど、面倒くさいのでこれまたループで素朴に。

import sys
from Bio import SeqIO

def count_gc(seq):
    return float((seq_record.seq.count("G") + seq_record.seq.count("C"))) * 100 / len(seq_record.seq)

if __name__ == '__main__':
    fasta = sys.argv[1]
    highest_id = ""
    highest_content = 0
    for seq_record in SeqIO.parse(fasta, "fasta"):
        c = count_gc(seq_record)
        if c > highest_content:
            highest_id = seq_record.id
            highest_content = c

    print highest_id
    print "{0:.6f}".format(highest_content)

僕はBioPerlでバイオインフォマティクス関連のプログラミングをしてたので、BioPythonで配列解析をしたことは殆ど無いんですよね。

あらかじめクラスタ数を決めないでクラスタリングする方法(Affinity Propagation)

K-meansのように予めクラスタ数を指定すると、「そのクラスタ数は正しいの?」っていう疑問が浮かぶと思う。

「なんらかの統計値に基づいて適切なクラスタに分割して欲しい」そんな願いを叶えるのがAffinity Propagationというクラスタリングアルゴリズムである

exemplara(セントロイドとかクラスタ中心)になるべきパラメータ(responsibility)とクラスタメンバに属しやすさ(availability)を交互に更新していって収束させる手法なので、K-meansのような初期値依存性がないらしい。

クラスタ数は類似度行列の対角要素(自分との類似度)に依存する(デフォルトはmedian)のでここを変更するとクラスタ数も変わるんだけどね。

Scikit-learnではAffinity Propagationが実装されているのでsykのケミストリースペースを作ってクラスタリングしてみた。ちなみにスライドのPCAの説明は間違っていた(pca.fit(fps).transform(fps)としなければいけなかった)。

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
from sklearn.decomposition import PCA
from sklearn.cluster import AffinityPropagation
from ggplot import *
import numpy as np
import pandas as pd

suppl = Chem.SDMolSupplier('syk.sdf')

fps = []
for mol in suppl:
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    fps.append(arr)

print len(fps)
pca = PCA(n_components=2)
x = pca.fit(fps).transform(fps)
af = AffinityPropagation().fit(x)

d = pd.DataFrame(x)
d.columns = ["PCA1", "PCA2"]
d["labels"] = pd.Series(af.labels_)
g = ggplot(aes(x="PCA1", y="PCA2", color="labels"), data=d) + geom_point() + xlab("PCA1") + ylab("PCA2")
ggsave("ap.png", g)

1404818907