ZDD Graphillionを使う 2 - ryamadaのコンピュータ・数学メモ

昨日の記事で、pythonのGraphillionパッケージを使ってみた
大まかなことが見えたので、整理し直して、チュートリアルの例から少し離れて、自分なりにGraphillionを使ってみるつもりで書いてみる
python自体にも不慣れなので、その辺りも確認しながら
まず、準備
- pythonが使える環境。Graphillionが使える環境。Windowsでも行けるはずだが、初期トライアルでこけたので、Macにて行う
- pythonモジュールをコンソールでインストール

sudo easy_install networkx # グラフオブジェクト用ライブラリ
sudo easy_install matplotlib # プロット用ライブラリ
sudo easy_install graphillion # Graphillion本体
sudo easy_install itertools # 列挙・組み合わせ関連のライブラリ

pythonの軌道とライブラリの読み込み

from graphillion import GraphSet # ライブラリのインストール
import graphillion.tutorial as tl
import networkx as nx
import matplotlib.pylab as plt
import itertools as it

世界を作る。Graphillionが扱うのは、全体集合としての無向グラフとそのサブグラフたちであるので、基本は全体集合としてのグラフの定義から始める
- グラフはノード集合とエッジ集合の組でできる。どんなものからスタートしてもよいが、ノード数を与え、その完全グラフ（すべてのノード間にエッジがある）を世界~universeとすれば、そのサブグラフを定義して、それを「小さな世界」とすることで任意のことができるので、まずはそこからやってみる

# 全ノードの全ペアを作るには、itertoolsを使うと便利
# import itertools as it # すでに読み込んである
n = 5 # ノード総数
# 全ペアを発生
# range(1,n) は[1,2,...,n-1]のこと
k = it.combinations(range(1,n),2)
# それをタプル(i,j)を要素とするリスト[(i,j),(i',j'),...]と言うエッジリストにする
universe = [] # エッジリストの格納オブジェクトを空からスタート
for element in k:
	universe.append((element))

universe

>>> universe
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]

このitertoolsというのは組み合わせ・順列などを整理してあるライブラリのようなので、少し横道にそれるが試してみる→こちら
完全グラフを世界にする用意ができたので、graphillionのグラフ世界を作ってみる

# from graphillion import GraphSet
GraphSet.set_universe(universe)

この操作により、GraphSetに「世界」が定まり、そこに定義されたエッジの集合の部分集合を部分グラフとしてオブジェクトとして扱えるようになる。これをしないと、GraphSetオブジェクト（世界グラフの部分グラフの集合）が作れない
グラフセット（部分集合）の作り方はいくつかある
- 部分グラフの集合として作る
- エッジを含める・含めないという条件によって全サブグラフについて選別してサブグラフ集合を作る
- 条件を満足するサブグラフの集合を作る点ではエッジ選別と同じだが、「パス」「クリーク」「木・全域木」「サイクル」などグラフ性質を指定して集合を作る
- 以上は使い勝手がよいように提供された方法だが、カスタマイズするための低レベル方法もある
エッジリストによって部分グラフを表し、その集合を作る

# エッジリストを与え、一つの部分グラフを作り、それを持って、要素数１の部分グラフ集合とする方法
graph1 = [(1,4)]
gs1 = GraphSet([graph1])
gs1 = GraphSet([graph1])
gs1

>>> gs1
GraphSet([[(1, 4)]])

# 部分グラフを複数与えれば、複数の部分グラフを要素とするグラフセットができる
gs2 = [(1,2),(2,3)]
gs2 = GraphSet([graph1,graph2])
gs2

>>> gs2
GraphSet([[(1, 4)], [(1, 2), (2, 3)]])

辺について条件をつけて、その条件に合うすべての部分グラフの集合を作る

# 持つ辺のリスト
edYes = [(1,2),(1,3),(2,3),(2,4)]
# 持たない辺のリスト
edNo = [(3,4)]
gs3 = GraphSet({'include': edYes, 'exclude': edNo})

- 2個しかない

for element in gs3:
	print(element)

>>> for element in gs3:
...     print(element)
... 
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4)]
[(1, 2), (1, 3), (2, 3), (2, 4)]

- {}は条件指定をするが、includeの指定もexcludeの指定もなければ、すべての部分グラフを持つグラフセットができる。今、世界は完全グラフなので、この場合には、ノード集合の冪集合になる

gs4 = GraphSet({})
for element in gs4:
	print(element)

>>> for element in gs4:
...     print(element)
...     
... 
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4)]
[(1, 2), (1, 3), (1, 4), (2, 3), (3, 4)]
[(1, 2), (1, 3), (1, 4), (2, 3)]
[(1, 2), (1, 3), (2, 3), (2, 4), (3, 4)]
[(1, 2), (1, 3), (2, 3), (2, 4)]
[(1, 2), (1, 3), (2, 3), (3, 4)]
[(1, 2), (1, 3), (2, 3)]
[(1, 2), (1, 3), (1, 4), (2, 4), (3, 4)]
[(1, 2), (1, 3), (1, 4), (2, 4)]
[(1, 2), (1, 3), (1, 4), (3, 4)]
[(1, 2), (1, 3), (1, 4)]
[(1, 2), (1, 3), (2, 4), (3, 4)]
[(1, 2), (1, 3), (2, 4)]
[(1, 2), (1, 3), (3, 4)]
[(1, 2), (1, 3)]
[(1, 2), (1, 4), (2, 3), (2, 4), (3, 4)]
[(1, 2), (1, 4), (2, 3), (2, 4)]
[(1, 2), (1, 4), (2, 3), (3, 4)]
[(1, 2), (1, 4), (2, 3)]
[(1, 2), (2, 3), (2, 4), (3, 4)]
[(1, 2), (2, 3), (2, 4)]
[(1, 2), (2, 3), (3, 4)]
[(1, 2), (2, 3)]
[(1, 2), (1, 4), (2, 4), (3, 4)]
[(1, 2), (1, 4), (2, 4)]
[(1, 2), (1, 4), (3, 4)]
[(1, 2), (1, 4)]
[(1, 2), (2, 4), (3, 4)]
[(1, 2), (2, 4)]
[(1, 2), (3, 4)]
[(1, 2)]
[(1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]
[(1, 3), (1, 4), (2, 3), (2, 4)]
[(1, 3), (1, 4), (2, 3), (3, 4)]
[(1, 3), (1, 4), (2, 3)]
[(1, 3), (2, 3), (2, 4), (3, 4)]
[(1, 3), (2, 3), (2, 4)]
[(1, 3), (2, 3), (3, 4)]
[(1, 3), (2, 3)]
[(1, 3), (1, 4), (2, 4), (3, 4)]
[(1, 3), (1, 4), (2, 4)]
[(1, 3), (1, 4), (3, 4)]
[(1, 3), (1, 4)]
[(1, 3), (2, 4), (3, 4)]
[(1, 3), (2, 4)]
[(1, 3), (3, 4)]
[(1, 3)]
[(1, 4), (2, 3), (2, 4), (3, 4)]
[(1, 4), (2, 3), (2, 4)]
[(1, 4), (2, 3), (3, 4)]
[(1, 4), (2, 3)]
[(2, 3), (2, 4), (3, 4)]
[(2, 3), (2, 4)]
[(2, 3), (3, 4)]
[(2, 3)]
[(1, 4), (2, 4), (3, 4)]
[(1, 4), (2, 4)]
[(1, 4), (3, 4)]
[(1, 4)]
[(2, 4), (3, 4)]
[(2, 4)]
[(3, 4)]
[]
>>>

- グラフセットの要素数はlen()関数でわかるので

len(gs4)

>>> len(gs4)
64

木、パス、サイクル、クリークなど
- ある「世界」から条件に合致したサブグラフを抜き出すのに、「世界」が完全グラフだとちょっと退屈なので、次のようにしてみる
- まず、完全グラフという世界から、連結グラフのグラフセットを作り、そのうちの一つを世界に置き換えることにする

# 連結グラフの取り出しは以下の通り
gconn = GraphSet.connected_components(range(1,n))
# ここからランダムに一つを取り出す
# GraphSetはpythonのset(集合)と同じ作りなので、以下のようにすればランダムに要素が取り出せる
# ランダムに取り出されたものはエッジを表すタプルのリストなので、それを指定して要素が１個のグラフセットとしてgrを作る
import random
gr = GraphSet(random.sample(gconn,1))

- どんなグラフなのかを図示するなら、grというグラフセットから要素を取り出し(１個しか要素がないので、取り出したものは、「それ自身」)、チュートリアルのdraw()関数を使う

tl.draw(gr.choice())

- ランダムに色々なグラフができているかを見たければ、以下を繰り返せばよい

tl.draw(GraphSet(random.sample(gconn,1)).choice())

- ちなみに、お絵描きはnetworkxを使っているので、グラフセットの要素をnetworkxのGraphオブジェクトにして、描くこともできる

nxg = nx.Graph(gr.choice()) # networkxのGraphオブジェクトにする
pos = nx.spring_layout(nxg) # 座標を計算し
nx.draw(nxg,pos) # ノードを描き
nx.draw_networkx_labels(nxg,pos) # ラベルを付与し
plt.show() # 描図する

- 「世界」を選ばれた連結グラフに置き換える

GraphSet.set_universe(gr.choice())

- クリークを持つ部分グラフの列挙

gclq = GraphSet.cliques(3)
for element in gclq:
  print(element)

- ルートノードを指定して木グラフの列挙

gtree = GraphSet.trees(1)
for element in gtree:
  print(element)

- - spanning treeオプションもある

gtree2 = GraphSet.trees(1,is_spanning=True)
for element in gtree2:
  print(element)

- Forest(森)は木の集合。is_spanning=Trueなら、複数のルートで全ノードを覆う場合が出せる(１ノードの木は出力しない？？？）

gforest = GraphSet.forests([1,2],is_spanning=True)
for element in gforest:
  print(element)

- サイクル、ハミルトンサイクル

gcycle = GraphSet.cycles()
for element in gcycle:
  print(element)
  
ghamilton = GraphSet.cycles(True)
for element in ghamilton:
  print(element)

- パス、ハミルトンパス

gpath = GraphSet.paths(1,2)
for element in gpath:
  print(element)
  
gpathH = GraphSet.paths(1,2,True)
for element in gpathH:
  print(element)

低レベルでカスタマイズ対応した作成法
- 上記の色々なサブグラフ集合の作成は、GraphSet.graphsという低レベル方法を使って書いたもの。それぞれがどのように書いてあるかを抜き書きしてみることで、GraphSet.graphsについて理解することにする
- graphset.pyに定義されている
- 連結

    @staticmethod
    def connected_components(vertices, graphset=None):
        """Returns a GraphSet of connected components.

        Args:
          vertices: A list of vertices to be connected.
          graphset: Optional.  A GraphSet object.  Components to be
            stored are selected from this object.
        Returns:
          A new GraphSet object.
        """
        return GraphSet.graphs(vertex_groups=[vertices], graphset=graphset)

- - これを見ると、univreseを変更しなくてもオプション指定のgraphsetに適当なグラフセットを与えれば、グラフセットの中で条件を満たすグラフを取り出してくれるようだ。その様式はgraphsの仕様なので、他のすべてのgraphsを用いた方法に共通する
  - 基本は１引数で、引数としてノードのリストを与えることで、graphsが返すのは、指定ノードのすべてが連結であるグラフが選ばれることがわかる
  - 『オプション vertex_groupsは連結成分を指定する』(リストが複数のリストからなっていれば、それぞれのリストが連結を、異なるリストに属するノードは非連結であることを要請する
- クリーク
  - エッジ本数をノード数kの完全グラフの本数としてneで指定している
  - 各ノードの次数は、クリークに含まれていれば、k-1、含まれなければ、0なので、[0,k-1]というリストを持たせる

    def cliques(k, graphset=None):
        """Returns a GraphSet of k-cliques.
        Args:
          k: An integer.  The number of vertices in a clique.
          graphset: Optional.  A GraphSet object.  Cliques to be
            stored are selected from this object.
        Returns:
          A new GraphSet object.
        """
        dc = {}
        for v in GraphSet._vertices:
            dc[v] = xrange(0, k, k - 1)
        ne = range(k * (k - 1) / 2, k * (k - 1) / 2 + 1)
        return GraphSet.graphs(vertex_groups=[[]], degree_constraints=dc,
                               num_edges=ne, graphset=graphset)

- 木
  - ルートノードに指定したノードごとに連結・非連結の分離をし
  - ノードの次数はspanning出なければ制約はなく、spanningならば、1以上、ノード数-1以下
  - ループは許さない→木の条件はこれか！

def trees(root=None, is_spanning=False, graphset=None):
        """Returns a GraphSet of trees.
        Args:
          root:  Optional.  A vertex, at which trees are rooted.
          is_spanning: Optional.  True or False.  If true, trees must
            be composed of all vertices.
          graphset: Optional.  A GraphSet object.  Trees to be stored
            are selected from this object.
        Returns:
          A new GraphSet object.
        """
        vg = [[]] if root is None else [[root]]
        dc = None
        if is_spanning:
            dc = {}
            for v in GraphSet._vertices:
                dc[v] = xrange(1, len(GraphSet._vertices))
        return GraphSet.graphs(vertex_groups=vg, degree_constraints=dc,
                               no_loop=True, graphset=graphset)

- 森：木が複数→ルートが複数で、ルートごとに連結・非連結を制御する
  - 木と違うのは、vertex_groups制約。リストのリストとし、それぞれの要素としてのリストは、ノードID１個を持った「長さ１のリスト」

    def forests(roots, is_spanning=False, graphset=None):
# ここが違う
        vg = [[r] for r in roots]
        dc = None
        if is_spanning:
            dc = {}
            for v in GraphSet._vertices:
                if v not in roots:
                    dc[v] = xrange(1, len(GraphSet._vertices))
        return GraphSet.graphs(vertex_groups=vg, degree_constraints=dc,
                               no_loop=True, graphset=graphset)

- サイクル
  - ハミルトンサイクルの場合にはすべてのノード次数が2
  - ハミルトンではないサイクルの場合には、ノード次数は0か2

    def cycles(is_hamilton=False, graphset=None):
        Args:
          is_hamilton: Optional.  True or False.  If true, cycles must
            be composed of all vertices.
          graphset: Optional.  A GraphSet object.  Cycles to be stored
            are selected from this object.
        Returns:
          A new GraphSet object.
        """
        dc = {}
        for v in GraphSet._vertices:
            dc[v] = 2 if is_hamilton else xrange(0, 3, 2)
        return GraphSet.graphs(vertex_groups=[[]], degree_constraints=dc,
                               graphset=graphset)

- パス
  - サイクルだけれど、始点、終点ノードだけ次数が1

    def paths(terminal1, terminal2, is_hamilton=False, graphset=None):
        Args:
          terminal1 and terminal2: Both end vertices of a paths.
          graphset: Optional.  A GraphSet object.  Paths to be stored
            are selected from this object.
        Returns:
          A new GraphSet object.
        """
        dc = {}
        for v in GraphSet._vertices:
            if v in (terminal1, terminal2):
                dc[v] = 1
            else:
                dc[v] = 2 if is_hamilton else xrange(0, 3, 2)
        return GraphSet.graphs(vertex_groups=[[terminal1, terminal2]],
                               degree_constraints=dc,
                               no_loop=True, graphset=graphset)

- ちなみにgraphsの引数説明は以下の通り

        Args:
          vertex_groups: Optional.  A nested list.  Vertices in an
            inner list are connected while those in different inner
            lists are disconnected.  For `[[1, 5], [3]]`, 1 and 5 are
            connected, while they are not connected with 3.
          degree_constraints: Optional.  A dict with a vertex and a
            range or int.  The degree of a vertex is restricted by the
            range.  For `{1: 2, 6: range(2)}`, the degree of vertex 1
            is 2 and that of 6 is less than 2, while others are not
            cared.
          num_edges: Optional.  A range or int.  This argument
            specifies the number of edges used in graphs to be stored.
            For `range(5)`, less than 5 edges can be used.
          no_loop: Optional.  True or False.  This argument specifies
            if loop is not allowed.
          graphset: Optional.  A GraphSet object.  Graphs to be stored
            are selected from this object.
          linear_constraints: Optional.  A list of linear constraints.
            A linear constraint consists of weighted edges and
            lower/upper bounds.  An edge weight is a positive or
            negative number, which defaults to 1.  Weights of the edges
            that are not included in the constraint are zeros.  For
            instance, `linear_constraints=[([(1, 2, 0.6), (2, 5),
            (3, 6, 1.2)], (1.5, 2.0))]`, feasible graph weights are
            between 1.5 and 2.0, e.g., `[(1, 2), (2, 3), (3, 6)]` or
            `[(1, 2), (2, 5), (5, 6)]`.
            See graphillion/test/graphset.py in detail.