2007-11-28

PythonでHarwell-Boeing matrix

C++のような本気言語ではなく、"バッテリー内蔵言語"であるPythonでHarwell-Boeing matrixのようなシビアな世界のデータ構造を書いたらどのぐらいになるかと試した。

このぐらいになる。


#!/usr/bin/python
# _*_ coding: euc_jp _*_
# Copyright (C) 2007 Ayukawa Hiroshi
from array import array
from itertools import izip, islice
from collections import defaultdict

class Hwmatrix(object):
def __init__(self, rownum_type="I", colnum_type="I", val_type="I"):
self.rowindex = array(rownum_type)
self.colindex = array(colnum_type)
self.val = array(val_type)
self.rowindex.append(0)

if val_type in ("f", "d"):
self.val_type = float
else:
self.val_type = int

def addrow(self, row):
row = dict(row)
cols = row.keys()
cols.sort()
vals = [row[k] for k in cols]
self.colindex.extend(cols)
self.val.extend(vals)
self.rowindex.append(self.rowindex[-1] + len(cols))

def getrow_as_pair(self, n):
rownum = self.rowindex[n]
rowlen = self.rowindex[n+1] - rownum
return islice(izip(self.colindex[rownum:], self.val[rownum:]), rowlen)

def __rmul__(self, vector):
if isinstance(vector, dict):
answer = defaultdict(self.val_type)
for k in vector.keys():
weight = vector[k]
for i, v in self.getrow_as_pair(k):
answer[i] += weight * v
return answer

def test():
"""
1 1 0 0 0
0 1 1 0 0
(0 2 3 0 1) * 0 0 1 1 0 = (0 2 5 3 1)
0 0 0 1 1
0 0 0 0 1
"""
mat = Hwmatrix()
mat.addrow({0:1, 1:1})
mat.addrow({1:1, 2:1})
mat.addrow({2:1, 3:1})
mat.addrow({3:1, 4:1})
mat.addrow({4:1})

print {1:2, 2:3, 4:1} * mat

test()


これはひょっとして実用的に使えたりしないだろうか。
会社ではC++で書いたものを使っていますが、自宅の趣味にはこういうのでもいいかも。

追記---
でも、よく考えると、ベクトルとの積__rmul__での処理なんかで、いちいち内部でPyObjectを作って参照カウンタをごにょごにょやっているんだと想像すると、一気に使う気が失せました。
やはりこういうのはC++じゃないといやだな。

2007-11-27

djangoのせいで

S2StrutsのクイックスタートのGreetingのサンプルを見たが、なんかめまいが。

djangoをここ3ヶ月ほど使っていたおかげで、もうJAVAでのWEB開発には戻れない体になってしまいました。

S2Strutsでの開発が自分に回ってきそうになっているが、どうしようか悩み中。

2007-11-17

ユークリッドの互除法

ユークリッドの互除法のような簡単なものだと、さすがにHaskellもPythonもそんなに差はなく、どちらもきれいに書ける。

Haskell

euc m 0 = m
euc m n | m < n = euc n m
| m `mod` n == 0 = n
| otherwise = euc n $ m `mod` n

main = do print $ euc 1071 1029


Python
def euc(m, n):
if m < n: return euc(n, m)
if n == 0: return m
if m % n == 0 : return n
return euc(n, m % n)

print euc(1071, 1029)

2007-11-16

[メモ]数学の記述ではやはりHaskellがPythonよりも上手ですが、変態的Pythonなら肉薄できそう

エラストテネスのふるいという素数生成アルゴリズムがあるが、これを(平方と素数の間に成り立つ不等式を考慮せずシンプルに)Haskellで書くと3行で書ける。

removenum x xs = filter (\a -> a `mod` x /= 0) xs
hurui (x: xs) = x : (hurui $ removenum x xs)

main = do print $ take 100 $ hurui [2..]

これを標準的な作法に従ったPythonで書くと少し行数が増える。(けっこうがんばったのですが。。。)
from itertools import ifilter, count, islice

def hurui(xs):
first = xs.next()
yield first
for x in hurui(ifilter(lambda x: x % first != 0, xs)):
yield x

for i in islice(hurui(count(2)), 100):
print i,

しかも、10000個ぐらい素数を出力させようとすると再帰呼び出しのスタックがあふれてしまう。
Haskellの場合は10000個でも大丈夫だし、きっといくつでも大丈夫なんだろう。(だからこそ関数型なんだろうし)

Pythonで無理矢理遅延評価を実現すれば、以下のようにも書けて、再帰であるのにスタックあふれもなくなる。
from itertools import ifilter, count, islice, takewhile, tee, dropwhile

def hurui(xs):
first = xs.next()
return (first, lambda : hurui(ifilter(lambda x: x % first != 0, xs)))

first, func = hurui(count(2))
for i in range(10000):
print first,
first, func = func()

こういうテクニックって、一般にはなんて言うんだろう。末尾再帰?


追記
コメントで頂いたのだが、以下のサイトによりエレガントな方法が。
http://4topcoder.blogspot.com/2007/07/python-generator.html
なるほどすごい。

2007-11-15

[メモ]MySQLのプロセス監視

これは便利。

watch -n 15 'mysql -u*** -p*** -e "show processlist"'

[メモ]MySQLの最適なテーブル定義を知る方法

既存のテーブルの無駄に長いフィールドなどを調整したい場合は以下のクエリーで最適なフィールドの長さなどの情報を表示できる。

SELECT * FROM table_name PROCEDURE ANALYSE()


参考: MySQL Presentations: Optimizing MySQL

2007-11-13

[日記]まだまだアッカーマン関数

いまだにアッカーマン関数が原始的でないことの証明をやっている。
あと1パターンの不等式を証明できればできあがりなんだけど。
ここの日記に最初にアッカーマン関数の事を書いてからもう8日がすぎています。
もしこれが大学の授業の宿題とかだったりしたら、自分は
「1週間やってたけどできませんでした。」
という、まじめはまじめなんだけど、できない子みたいだなあ。

すべて完成したら、次にもう一度この証明を見直すときのために証明の要点をここの日記に整理して書こうと思います。
しかし、まずは残りをやらなければ。

2007-11-09

[雑談]「これから一からプログラミングを覚えようと考えています。」

これから一からプログラミングを覚えようと考えています。
様々な言語がありますが、どれを覚えるべきでしょうか?


PHPすごい人気ですね。
初心者に最適という評価はあるようだ。

でも初心者でなくなったときにPHPに感じる落胆もでかいんだよなあ。

数年PHPを使い込んだ過去のある自分としてはおすすめできない。
LL等のイベントでPHPの重鎮の発表の自虐的な雰囲気を見てから判断してもいいと思った。

(このブログにトラックバック機能があれば、トラックバックしたかったが、残念ながらないみたい。。)

Pythonはこういう時に自信を持って勧めていいものかどうか迷います。
他の言語にあらかた疲れた時に、ふと使ってみるとはまるという性質もあるからねえ。

2007-11-08

[日記]まだアッカーマン関数

毎日少ししか勉強時間が取れないので、アッカーマン関数の演習問題から抜け出せていない。

2日悩んだ問題「アッカーマン関数をNプログラムで計算する方法を考えよ」は、スタックをNプログラムで書ければうまくいきそうだとわかったが、それには次の章以降の知識がいるようだ。。。スタックを書く方法を一日考えていたがやはり自力では思いつかないので素直に次章へ進むしかないか。
しかし、「考えよ」ということだから、(「書け」というわけではないから、)ここまでできれば完了なんだろうか。なんともすっきりしない。

次の問題はアッカーマン関数が原始的関数でないことの証明問題なので、着々と進んでいるが、長いので時間がかかりそうだ。

2007-11-06

アッカーマン関数やばい

アッカーマン関数を計算論の教科書の演習問題で見て、問題がうまく解けないのでためしにアッカーマン関数の計算結果の一部を一覧にして眺めようと思った。
でも、計算量が多すぎてすぐにとんでもないことになるんだね。お願いだから教科書に書いておいてくださいと思いました。

f 0 y = y + 1
f x 0 = f (x - 1) 1
f x y = f (x - 1) (f x (y - 1))

main = do print $ f 4 2

これをHaskellで実行すると死ぬ。
f 3 2 ぐらいはまだ大丈夫だが。

どのぐらい計算量がやばいことになるかを、f 2 2 の場合で置換モデルでノートに書き下していったら、なんと27行もノートを使うということが判明。
f 4 2なんか、そりゃとんでもないことになってるんだろうなあ。

[メモ]正規表現モジュール re のscanner機能

Pythonのドキュメントには書かれていないが、正規表現モジュール(re)にはscannerという機能がある。
一旦作成した正規表現オブジェクトで、長い文章中をスキャンしていくような感じで、全てのマッチする部分を検出していく機能だ。
便利なんだが、使いたい時に毎回使い方を忘れているのでこの際ここにメモしておこう。

>>> a="aaabbbcccdddaaa111fffsss"
>>> import re
>>> pat = re.compile(r"aaa(...)")
>>> dir(pat)
['__copy__', '__deepcopy__', 'findall', 'finditer', 'match', 'scanner',
'search', 'split', 'sub', 'subn']
>>> scanner = pat.scanner(a)
>>> dir(scanner)
['match', 'search']
>>> m = scanner.search()
>>> m
<_sre.SRE_Match object at 0x2a98006210>
>>> m.groups()
('bbb',)
>>> m = scanner.search()
>>> m.groups()
('111',)
>>> m = scanner.search()
>>> m
>>> #もう候補がないので m = None となっている。

2007-11-01

[メモ]有機EL

EUCで「有機EL」という語の「機」の後半バイトと「E」の前半バイトはくっつけると「。」と同じコードになる。

>>> s = lambda x: map( lambda y: hex(ord(y)) , unicode(x, "utf-8").encode("euc_jp"))
>>> s("有機EL")
['0xcd', '0xad', '0xb5', '0xa1', '0xa3', '0xc5', '0xa3', '0xcc']
>>> s("。")
['0xa1', '0xa3']