ZDDグラフのノード値を行列計算する

[statgenetJimu]

以下をpythonとRとで確かめる。 出来れば、疎行列(sparse Matrix)を使って大規模なグラフでもやってみる。 Rならグラフ用パッケージigraphなどを適宜活用するとよいはず。pythonのグラフ系モジュールに何があるのかはよく知らない。

ZDDのグラフとして、 v1 -(lo)->v2, v1 -(hi)->v3, v2-(lo)->T, v2-(hi)->T, V3-(lo)->B,V3-(hi)->Tとなっているとする。

このとき、v1の値を1とすると、v2=v1_1,v3=v1_1,T= v2_2+V3_1,B=v3*1という多変数の連立一次方程式になっている。 これを線形代数で表すと、vi -> vjのエッジの本数を要素としたノード数 x ノード数の正方行列ができる。

今、この正方行列を(1,0,0,0,0)という、第一成分が1で残りの成分が0のベクトルに何度も左から掛けていくと、最後に、v1,v2,...,T,Bを通るパスの本数が求められる。

ベクトルに何度も掛けるかわりに、行列のべき乗をしてやれば、その第１列が答えともなる。

|python| import numpy as np M = mat([[1,0,0,0,0],[1,0,0,0,0],[1,0,0,0,0],[0,1,1,0,0],[0,1,1,0,0]]) M \ 5 ||< したがって、ZDD情報から、行列を作ってやればよさそう。

[skoyama427]

例に上がっているZDDですが

v1 -(lo)-> v2 v1 -(hi)-> v3 v2 -(lo)-> B <==ここBで良かったですか？ v2 -(hi)-> T V3 -(lo)-> B V3 -(hi)-> T

これならMAT**5の1列目は(1,1,1,2,2)になります。 手で数えても正しいと思います。

[statgenetJimu]

提示例はあくまでも例ですので、気にせずやってください。実際、この行列のべき乗計算で正しいかどうかの確認を取っていないので、その確認もお願いします。

[skoyama427]

しばらく考えてみたのですが、一度作成されたZDDから上記のようにして情報を取り出せるでしょうか？以下、間違いなどありましたらご容赦のうえ、ご指摘ください。

上記のグラフ v1 -(lo)-> v2 v1 -(hi)-> v3 v2 -(lo)-> T v2 -(hi)-> T V3 -(lo)-> B V3 -(hi)-> T で表されるひとつの例として {a, b} = [{0,1}, {0,0}, {1,1}]のハプロタイプ集合が考えられます。 graphillionからこれをZDD化し.dumpsで観察(p01.py)すると 2 2 B T 5 1 T 2 が返ってきました。ノードが2つになっています。これは、今回のグラフセットにおいてaの値が0の時、bはいずれの枝もTを指すため冗長とみなされてZDDのデータ構造に変換する際には削除されることを表しているのではないでしょうか。これから考えるに作成されたZDDはノードを通過するパスの量的情報を失っているように思われます。とするとハプロタイプ集合から作成されたZDDの情報からノードを通過するパスの数を計算するのは難しいように思います・・・。

それとも、全く論点がずれていたりしますでしょうか？？？

[statgenetJimu]

[{0,1}, {0,0}, {1,1}]の例が作るZDDは

print(G1.dumps()) 2 2 B T 3 2 T T 11 1 3 2 .

ではないかと思うのですが。 以下、実行結果です。

from graphillion import GraphSet import pickle from itertools import chain

import networkx as nx import matplotlib.pyplot as plt

haps = [[0,0],[0,1],[1,1]] n = len(haps[1]) ng = len(haps)

universe = []

for i in range(n): universe.append((i,i+1))

GraphSet.set_universe(universe)

g_univ = GraphSet({})

g1 = [] i = 1 for i in range(ng): tmp = [] for j in range(n): if haps[i][j] == 1: tmp.append(universe[j]) if i == 0: g1 = [tmp] else : g1.append(tmp)

G1 = GraphSet(g1) print(G1.dumps())

[skoyama427]

そうでした。。。 小規模ですがRで行列計算のアルゴリズムを作ってみました。。。 出てくる値は各ノードに入射してくるノードの数を表しているようです。。。