fc2ブログ
ダイクストラ法ですべての最短経路を求める
Pythonの話というより、グラフ理論アルゴリズムの話になりつつあって、ちょっとマニア度が高くなってきましたが、NetworkXのライブラリに少し手を入れて、すべての最短経路を求められるようなメソッドを作って見ました。
重みつき無向グラフの最短経路(shortest path)を求めるアルゴリズムとして、ダイクストラ法(Dijkstra's algorithm)は有名ですが、これを実装したNetworkXのpathsモジュールのメソッドsingle_source_dijkstra(G,source,target=None): は、sourceから、あるtargetまでの最短経路のうち1つしか返してくれません。同じ最短経路長で、いくつかの経路がある場合、それが全部知りたいと思うのは、自然な発想かと。
また、ダイクストラ法は、そのアルゴリズムの本質的な部分で、すべての可能な最短経路を計算しているので、改変はいたって簡単。ソースを貼り付けておきます。

def single_source_dijkstra_all_paths(G,source,target=None):

 if source==target: return (0, [source])
 dist = {} # dictionary of final distances
 paths = {source:[[source]]} # dictionary of paths
 seen = {source:0}
 fringe=[] # use heapq with (distance,label) tuples
 heapq.heappush(fringe,(0,source))
 if not G.is_directed(): G.successors=G.neighbors
 # if unweighted graph, set the weights to 1 on edges by
 # introducing a get_edge method
 # NB: for the weighted graphs (XGraph,XDiGraph), the data
 # on the edge (returned by get_edge) must be numerical
 if not hasattr(G,"get_edge"): G.get_edge=lambda x,y:1

 while fringe:
  (d,v)=heapq.heappop(fringe)
  if v in dist: continue # already searched this node.
  dist[v] = seen[v]
  if v == target: break
   
  for w in G.successors(v):
   vw_dist = dist[v] + G.get_edge(v,w)
   if w in dist:
    if vw_dist < dist[w]:
     raise ValueError,\
       "Contradictory paths found: negative weights?"
   elif w not in seen or vw_dist < seen[w]:
    seen[w] = vw_dist
    heapq.heappush(fringe,(vw_dist,w))
    for _each_path in paths[v]:
     paths.setdefault(w,[]).append( _each_path + [w] )
   elif vw_dist == seen[w]:
    for _each_path in paths[v]:
     paths[w].append( _each_path + [w] )
 return (dist,paths)

pathsの値に、経路のリストを保持するようにして、ループの際に、これまで見つかっている経路長と同じでも、その経路を保持するようにしただけです。

小規模なグラフで動くことは確認したけど、大規模になったときどうかは、確認する必要があるかも。
スポンサーサイト



【2006/10/17 21:03】 | Python | トラックバック(0) | コメント(2) | page top↑
psycoで高速化
Python使っている人には当たり前なのかもしれないけど、なんか実行を高速化してくれるツールがあるらしい。Psycoがそれです。とりあえず、ダウンロードでインストール。windows用のバイナリがあったので、バージョン1.5を選択。psyco-1.5.win32-py2.4.exe

仕組みは少しこれから詳しくなる必要がありそうだけど、使い方は簡単。
import psyco
psyco.full()
と、スクリプトの冒頭に入れるだけ。

Cと同じくらい速いとか、実行速度が10倍になるとか書いてるので、ネットワークの中心性(centrality)を計算するnetworkx.centrality.betweenness_centrality( Graph )を使ってプログラムの実行速度を調べてみました。使ったマシンは、PentiumD 3.2GHz(2G RAM) WindowsXPsp2
XPだと、Microsoftが出しているResource Kitに含まれるtimeitというツールを使える。Windows Server 2003 Resource Kit Tools
C:\Program Files\Windows Resource Kits\ToolsにPATHを通すの忘れずに。

node数5897、edge数15950のネットワークのBetweennessによるCentralityを求める。確かに速度が2倍くらいに。
とりあえず、時間がかかる計算をするときは試してみる価値あるね。ループがないと意味は薄いようだけど。たとえば、巨大なリストからのランダムサンプリングとかは、ダメな模様。
以下は参考までに、timeitの出力結果。

psyco適用前
Version Number: Windows NT 5.1 (Build 2600)
Exit Time: 8:05 pm, Monday, October 16 2006
Elapsed Time: 0:09:15.562
Process Time: 0:09:12.031
System Calls: 1228887
Context Switches: 256903
Page Faults: 815529
Bytes Read: 105270080
Bytes Written: 11654265
Bytes Other: 3072382

psyco適用後
Version Number: Windows NT 5.1 (Build 2600)
Exit Time: 8:15 pm, Monday, October 16 2006
Elapsed Time: 0:05:33.968
Process Time: 0:05:33.062
System Calls: 2064202
Context Switches: 465510
Page Faults: 1811956
Bytes Read: 6464132
Bytes Written: 808205
Bytes Other: 1981950

しかし、Psycoって何の略なんだ?
【2006/10/17 01:02】 | Python | トラックバック(0) | コメント(0) | page top↑
networkxとmatplotlib
複雑ネットワークを解析するプログラムをPythonで書きたいと思っていて、Java言うyFilesのように、グラフとしてのネットワークの特性を扱うライブラリが無かったら、自分で書こうと思っていたけど、ものすごくいいのがありました。それが、networkXで、もしかしたらyFilesより高機能かも。しかも、フリー。yFilesはアカデミアでも20万くらいとるから、Pythonやっぱり万歳か?

利用環境はPython2.4.3で、Windowsのバイナリなら、installはいたって簡単。最新版である0.32を入れました。networkx-0.32.win32.exe
これだけで、ネットワークを扱うプログラムはすぐに書けるようになりますが、描画できる方がいいかと思い、matplotlibも導入。matplotlib-0.87.6.win32-py2.4.exe

>>> import pylab
Traceback (most recent call last):
File "", line 1, in ?
File "C:\Python24\Lib\site-packages\pylab.py", line 1, in ?
from matplotlib.pylab import *
File "C:\Python24\Lib\site-packages\matplotlib\pylab.py", line 197, in ?
import cm
File "C:\Python24\Lib\site-packages\matplotlib\cm.py", line 5, in ?
import colors
File "C:\Python24\Lib\site-packages\matplotlib\colors.py", line 33, in ?
from numerix import array, arange, take, put, Float, Int, where, \
File "C:\Python24\Lib\site-packages\matplotlib\numerix\__init__.py", line 73,
in ?
import numpy
ImportError: No module named numpy

数学ライブラリがが必要なよう。matplotlibのサイトに、

For standard python installations, you will also need to install either Numeric or numarray in addition to the matplotlib installer. matplotlib provides installers for Numeric and numarray users. matplotlib has a numerix setting in the matplotlib rc file (which by default resides in c:\python23\share\matplotlitb\matplotlibrc) and you should make sure this setting corresponds to your preferred array package.

と、書いてるので、先日Rとの連携で入れたNumericをそのまま使うことに。Python界の流れ的には、numarrayが主流になりつつあるようで、matplotlibのデフォルトもそれを使うようになっていますが、
C:\Python24\Lib\site-packages\matplotlib\mpl-data
の、matplotlibrcを編集。
numerix : numpy # numpy, Numeric or numarray
となっているところをコメントアウトし、
numerix : Numeric
に変更することで、無事importできるようになります。

networkxのエライところの1つに、何でもnodeにできるところがあります。
たとえば、以下のようにtupleをnodeにできる。そのまま描画するとわかりやすくていいです。

>>> import networkx as nx
>>> import pylab
>>> g = nx.Graph()
>>> g.add_edge( ('node1',3) , ('node2', 8))
>>> g.add_edge( ('node2',8) , ('node3', 9))
>>> g.add_edge( ('node2',8) , ('node4', 5))
>>> nx.draw( g )
>>> pylab.savefig('net.png')

簡単なネットワークの図

【2006/10/16 15:52】 | Python | トラックバック(0) | コメント(0) | page top↑
PythonからRを呼ぶ
いくら、Pythonに充実した標準ライブラリが備わっているといえども、統計処理をしようと思うと、なかなかすぐにできるものではありません。順位和検定をしたいと思い、Pythonのライブラリを少し探したけど、見つからないのでR(http://www.r-project.org/)との連携を模索したらすぐにいいのが見つかりました。
RPy(http://rpy.sourceforge.net/
現時点で、RC1が発表されていました。2.4系に対応しています。ちなみに試した環境は、2.4.3
RPyは、rpy-1.0-RC1.win32-py2.4.exeを選択。
インストールは、言われるがままに、
Python for Windows extensions(pywin32-210.win32-py2.4.exe
Numerical Python(Numeric-24.2.win32-py2.4.exe
をそれぞれ、インストール。R.dllのある場所を、windowsの環境変数PATHに追加。これで、OK!と思いきや、実行すると
>>> import rpy
RHOME= C:\Program Files\R\R-2.3.1
RVERSION= 2.3.1
RVER= 2031
RUSER= C:\Documents and Settings\tsuji
Loading the R DLL C:\Program Files\R\R-2.3.1\bin\R.dll .. Done.
Loading Rpy version 2031 ..

Traceback (most recent call last):
File "", line 1, in -toplevel-
import rpy
File "C:\Python24\Lib\site-packages\rpy.py", line 112, in -toplevel-
exec("import _rpy%s as _rpy" % RVER)
File "", line 1, in -toplevel-
ImportError: No module named _rpy2031

といったようなメッセージが出る。ここでは、エラーを再現するためにわざとやっているけど、結局、rpyがインストールされたフォルダ(たとえば、C:\Python24\Lib\site-packages)にモジュール(_rpy2031.pyd)がないとダメみたい。今回使ったRC1は、デフォルトでR2.3.1に対応しているので、Rのバージョンをこれにあわせておくと、楽かも。

import rpy

wc_result = rpy.r.wilcox_test( data1 , data2, paired = rpy.r.TRUE );
wc_result.get( "p.value" )

で、計算結果のP値を得られる。TRUE,FALSEは、そのまま文字列では書けないので注意。
【2006/10/14 20:51】 | Python | トラックバック(0) | コメント(0) | page top↑
| ホーム |