2015-02-06

pythonでDiscrete Exterior Calculus(1)

このパッケージ説明文書を読んでみる
pythonパッケージはこちら
その第１回
冒頭PDFのpp1-11
- 1 Introduction
- 2 Overview of PyDEC
- 3 Simplicial Complex Representation
- 4 Regular Cube Complex Representation
- 5 RipsComplex Representation
- 6 Abstract Simplicial Complex Representation
Simplicial Complex Representation
- Simplicial Complexはこちらにあるように単体的複体
- 三角形とその一般化である単体の集合として表されるもの
- 複体はトポロジー的にはこちらにあるような定義もある
- これをPyDECでどのように表すか
- 点は、0-単体、辺は1-単体、三角形は2-単体、…
- p-単体の集合を取り、複体はその集合と考える
- p-単体集合はを行列表現する
  - 行列の各行は１個のp-単体を表す
  - １個のp-単体は、p+1個の0-単体の順序を気にする並びであるとする(この順序は、辺ならその向き、三角形なら、面の向き付けと関係する)
  - 0-単体の順序は、巡回に関して同一視する
  - 同一視するp-単体には、辞書式順序により順番を決めることで表現の一意性を確保する
Cube Complex Representation
- 正方形とその一般化の集合として表す
- 正方形・立方体・超立方体は、基準点座標と、広がる方向軸とによって表す
- p-Cubeの集合をp=0,1,2,...について集めたものがCube complex
- n次元空間におけるp-cubeは、基準点座標としてn次元空間座標を持ち、広がる方向としてp方向を持つ
- 今、n=3、p=2のp-Cubeは、３次元空間に置かれた正方形であり、それは、たとえば、*1のように表す。(x=0,y=1,z=3)が基準点座標であり、(0,1)は0=x軸、1=y軸方向を表すから、この正方形は、(0,1,3),(1,1,3),(1,2,3),(0,2,3)を頂点とする正方形を表す
- p-Cubeの集合は、このn+pの数値列を1個のp-Cubeに対応させ、それを1行にした、n+p列の行列として表す
Rips Complex Representation
- 空間に座標を持った点集合があったときに、点間ペアワイズ距離が閾値距離以下である点間にエッジを引き、それが作るグラフを単体的複体とみなす
- エッジはエッジリスト表現(２列の行列)にもできるし、正方行列表現もできる、また、各エッジの端点に１を立てた、辺数ｘ点数の行列でも表現できる
- 正方行列表現の行列と変数ｘ点数表現した行列の行列積をとると、そこに2-単体(三角形)がいくつ、どの点の組み合わせ存在するかが線形代数計算で求まる。さらに同様の計算が、再帰的に高次の単体集合の検出にも使える
Abstract Simplicial Complex Representation
- Rips Complexは位置座標を持つ点集合に関するものであったが、点間距離を忘れ、単なる、点が作るグラフ構造(抽象構造)を複体を考えるのがAbstract...なる考え方
- 特に、点を２回描画してもよい(ただし同一の点には、同一の点IDをつける)ことにすれば、ぐるりと巻いた帯と、メビウス帯とが、ほぼ同一の描画表現にできて、唯一の違いは、両端に相当する(実際は一つしかない点だが、便宜上、２回描いてある点)のIDのつけ方の違いとして表現できる。行列の集合としての複体表現では、「描く」必要もないので、両者を表現し分けられている

import numpy

PyDECのExamplesのうちのDancyFlowのdriver.pyをSpyderで実行してみる
- Examplesなどのとり方はこちら
- 処理の冒頭を確認する
  - 点集合(座標)のファイルvertices.txtを読み込んでverticesオブジェクトを作り
  - 三角形集合のファイルtriangles.txtを読み込んでいる。「点ID」なのでint型で読み込み、さらに、ファイルでは点IDが1,2,...になっているが、pyDECでは点は0,1,2,...とID付けするので-1している
  - simplicial_complex()は点の座標情報と三角形を表す行列表現

velocity = array([1.,0])
# Read the mesh
vertices = loadtxt('vertices.txt')
triangles = loadtxt('triangles.txt', dtype='int') - 1
# Make a simplicial complex from it
sc = simplicial_complex((vertices,triangles))

- これをファイル読み込みでなく実施してみると

vertices = array([[1,2],[1,3],[2,4],[1,4]])
triangles = array([[0,1,2],[0,1,3]])
sc = simplicial_complex((vertices,triangles))
vertices
triangles
sc

- vertices,trianglesが２次元ｎアレイであることがわかり、scというオブジェクトは複体を表しており、0単体が4個(点が４つ)、1単体が5個(辺が５個)、2単体が２個(三角形が２個)あることがわかる

vertices
Out[12]: 
array([[1, 2],
       [1, 3],
       [2, 4],
       [1, 4]])
triangles
Out[14]: 
array([[0, 1, 2],
       [0, 1, 3]])

sc
Out[15]: 
simplicial_complex:
  complex:
            2:  2-simplices
            5:  1-simplices
            4:  0-simplices

- Examplesの例にならって、scオブジェクトから0,1,2単体の個数を表示させてみれば

sc[0].num_simplices
sc[1].num_simplices
sc[2].num_simplices

sc[0].num_simplices
Out[18]: 4

sc[1].num_simplices
Out[19]: 5

sc[2].num_simplices
Out[20]: 2

- pythonの情報表示機能を使えば、次のようにもできる

sc?
Type:        simplicial_complex
String form:
simplicial_complex:
  complex:
            2:  2-simplices
            5:  1-simplices
            4:  0-simplices
Length:      3
File:        c:\python27\lib\site-packages\pydec\dec\simplicial_complex.py
Docstring:
simplicial complex

This can be instantiated in several ways:
    - simplicial_complex( (V,S) )
        - where V and S are arrays of vertices and simplex indices
    - simplicial_complex( M )
        - where M is a simplicial_mesh object


Examples
========
>>> from pydec import simplicial_complex, simplicial_mesh
>>> V = [[0,0],[1,0],[1,1],[0,1]]
>>> S = [[0,1,3],[1,2,3]]
>>> M = simplicial_mesh(V, S)
>>> simplicial_complex( (V,S) )
>>> simplicial_complex( M )

*1:0,1,3),(0,1