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++じゃないといやだな。

djangoのせいで

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

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

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

ユークリッドの互除法

ユークリッドの互除法のような簡単なものだと、さすがに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)

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

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

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
なるほどすごい。

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

これは便利。

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

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

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

SELECT * FROM table_name PROCEDURE ANALYSE()

参考: MySQL Presentations: Optimizing MySQL

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

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

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