グラフの経営への応用
グラフの諸定理はすこ～しずつ理解できてきました。問題はこれの応用です。工学や交通量などへの応用は多いようなのですが、経営問題に応用例がみつかりません。工程管理で応用している例は経営数学の入門書でみたことがあるのですが、その程度です。小林みどり先生という方が数学を専門としつつも経営への応用を試みているようです。
応用例で参考にならないかと思って書籍を取り寄せることにしました。2010/12/08 19:33:20 PST
グラフの視覚化
グラフを頂点、辺の集合としてリスト表記する方法をとっていて、必要に応じて行列との相互変換をしています。コンピューター内部ではそれでいいのですが、人間にとっては画面で視覚化されていた方が扱いやすいです。自分でそんなのをつくるととても手間がかかるしどうしようかと思っていたら、すでにあるらしいです、そういうもの。
Graphviz だそうです。Gaucheからの変換関数を書けば視覚化も容易そうです。 2010/12/06 22:16:21 PST
簡単なのでdot言語、試してみました。
(define (graphviz g)
(let ((e (cadr g)))
(call-with-output-file "graph.dot"
(lambda (out)
(display "digraph sample {" out)
(newline out)
(display "rankdir = LR;" out)
(newline out)
(let loop ((edge (car e)) (rest (cdr e)))
(if (null? rest)
(display "}" out)
(begin
(display (symbol->string (cadr edge)) out)
(display " -> " out)
(display (symbol->string (caddr edge)) out)
(display ";" out)
(newline out)
(loop (car rest) (cdr rest)))))))))
これでgraph.dotをつくっておいて
gosh> (run-process "dot -Tgif graph.dot -o graph.gif")
*
** SYSTEM-ERROR: cannot find program 'dot -Tgif graph.dot -o graph.gif': 指定されたファイルが見つかりません。
Stack Trace:
_______________________________________
0 (sys-fork-and-exec (car argv) argv :iomap iomap :directory dir :si ...
At line 116 of "C:\\Gauche\\share\\gauche\\0.9\\lib/gauche/process.scm"
gosh>
Windowsだからでしょうか。cmd.exeから直接、コマンドを送るとｏｋで、無事グラフは表示されました。2010/12/07 01:21:05 PST
(use gauche.process)
(load "./graph-test.scm")
(define (graphviz g)
(let ((e (cadr g)))
(call-with-output-file "graph.dot"
(lambda (out)
(display "digraph sample {" out)
(newline out)
(display "rankdir = LR;" out)
(newline out)
(let loop ((edge (car e)) (rest (cdr e)))
(display (symbol->string (cadr edge)) out)
(display " -> " out)
(display (symbol->string (caddr edge)) out)
(display ";" out)
(newline out)
(cond ((null? rest)
(display "}" out))
(else
(loop (car rest) (cdr rest)))))))))
(define (make-gif)
(run-process '("dot" "-Tgif" "graph.dot" "-o" "graph.gif")))
(define (make-dot g)
(graphviz g)
(make-gif))
gosh> (make-dot g8)
#<process #<win:handle process 3040 @0x117c670> "dot" active>
gosh>
キャッシュフローとグラフ理論
仕事上、応用したいのはキャッシュフロー。マトリクス会計とか言いつつもその方面の 人たちはほとんどが文系出身者で線形代数もグラフも知らない人たちで占められています。 財務諸表で採用されているキャッシュフロー計算書は原始時代の石斧のレベル。 なんとか応用を見出し、科学的な分析ツールにしたいものです。2010/12/06 12:13:42 PST
最大フロー問題
「情報科学のためのグラフ理論」加納幹雄 著 朝倉書店、 ここにFord,Fulkersonのアルゴリズムの詳細説明があります。
ライブラリはVer0.17で正常動作しているはずです。 以前に書いたものは誤りを含んでいて検索阻害要因になるので削除します。
wikipediaの↑の解説は間違ってると思います。アルゴリズムの記述の部分、例のステップ数の箇所。
;;最大フロー問題有向グラフ
;;有向グラフでは弧(arc)をa1,a2...で表す
;;「グラフ理論序説」プレアデス出版、ｐ88
(define g6 '((v1 v2 v3 v4) ((a1 v1 v2 4)(a2 v2 v3 1)(a3 v1 v3 2)(a4 v2 v4 3)(a5 v3 v4 4))))
;;流れない経路を考慮しないと正しく計算されない例
(define g7 '((v1 v2 v3 v4) ((a1 v1 v2 4)(a2 v2 v3 1)(a3 v1 v3 4)(a4 v2 v4 4)(a5 v3 v4 4))))
;;p89の例,最大フロー１０
(define g8 '((v1 v2 v3 v4 v5 v6 v7 v8)
((a1 v1 v3 4)
(a2 v1 v5 4)
(a3 v1 v7 4)
(a4 v2 v3 4)
(a5 v2 v6 3)
(a6 v3 v4 3)
(a7 v3 v6 2)
(a8 v5 v6 3)
(a9 v5 v8 2)
(a10 v6 v4 3)
(a11 v7 v5 2)
(a12 v7 v8 2)
(a13 v8 v6 2)
(a14 v8 v4 4))))
;;ネット上の解説資料の例題
(define g9 '((v1 v2 v3 v4 v5 v6)
((a1 v5 v1 60)
(a2 v1 v2 30)
(a3 v1 v3 20)
(a4 v1 v4 10)
(a5 v3 v4 50)
(a6 v4 v2 40)
(a7 v2 v6 20)
(a8 v4 v6 50)
(a9 v5 v3 10)
(a10 v3 v2 10))))
;;情報科学のためのグラフ理論p122
(define g10 '((v1 v2 v3 v4 v5 v6)
((a1 v5 v1 4)
(a2 v5 v4 4)
(a3 v4 v1 1)
(a4 v4 v2 2)
(a5 v1 v2 3)
(a6 v4 v3 2)
(a7 v3 v2 1)
(a8 v3 v6 3)
(a9 v2 v6 4))))
gosh> (maximum-flow g6 1 4)
6
gosh> (maximum-flow g7 1 4)
8
gosh> (maximum-flow g8 1 4)
10
gosh> (maximum-flow g9 5 6)
60
gosh> (maximum-flow g10 5 6)
6
gosh>
2010/12/06 00:49:43 PST
ダイクストラ、最短経路問題アルゴリズム
ライブラリに追加しました。シンプルだけれど強力なアルゴリズムみたいですね。 こないだ車につけたナビゲーターもあれで最短経路を見つけてるのかいなぁ。2010/11/30 07:31:02 PST
グラフのスペクトル分析
エール大学の先生が研究されてるらしいです。 http://www.cs.yale.edu/homes/spielman/sgta/SpectTut.pdf 本で読んだのは序説なのですが、上記はかなりつっこんだ研究なようです。 理解できるかどうかはわかりませんが、なんだかワクワクします。2010/11/29 14:55:51 PST
グラフ理論ライブラリ
ものすご～く不完全ですが、グラフを扱う関数を書き始めました。主にグラフの固有値について興味がありそのあたりを書き進めたいと思っています。
http://homepage1.nifty.com/~skz/Entry/graph.html
2010/11/29 08:26:17 PST
グラフにはまる
グラフ理論の本にある諸定理をGaucheで試していてその不思議さに魅了されはじめました。グラフは頂点と辺の集合ですからこれをリストで表して、計算はコンピューターにとっては行列の方が楽そうなので隣接行列にして。な～んてことでコードを書きま始めました。グラフ理論には美しい定理があれこれと転がっていそうです。RubyにはグラフのライブラリがあるようなのでGaucheにだってあってもいいはず。まあ、アマチュアの私が書くものだから品質には疑問符がつきますが。できあがったらご報告します。マトリクス会計のコード書きはさらに遅延（^^;
2010/11/26 04:50:23 PST
ますます不思議
Σ(i=1 to p) λi^2 = tr(A^2) = 2|E(G)|
(define k (list->matrix '((0 1 1 1)
(1 0 1 1)
(1 1 0 1)
(1 1 1 0))))
gosh> (matrix-jacobi k)
(-1.0 -1.0 -1.0 3.0)
gosh> (matrix-tr (m-expt k 2))
12
gosh>
(-1)^2+(-1)^2+(-1)^2+3^2 = 12
４つの頂点の完全グラフなので辺は６
へ～～～～。2010/11/25 07:15:12 PST
閉路グラフの固有値
へ～、不思議なものです。えっ！複素数のオイラーの公式もひょっとしたら関係あるんですか。
(define g (list->matrix '((0 1 0 0 0 1)
(1 0 1 0 0 0)
(0 1 0 1 0 0)
(0 0 1 0 1 0)
(0 0 0 1 0 1)
(1 0 0 0 1 0))))
(define h (list->matrix '((0 1 0 0 0 0 0 1)
(1 0 1 0 0 0 0 0)
(0 1 0 1 0 0 0 0)
(0 0 1 0 1 0 0 0)
(0 0 0 1 0 1 0 0)
(0 0 0 0 1 0 1 0)
(0 0 0 0 0 1 0 1)
(1 0 0 0 0 0 1 0))))
gosh> (matrix-jacobi g)
(-2.0 -1.0 -1.0 1.0 1.0 2.0)
gosh> (matrix-jacobi h)
(-2.0 -1.414213562373095 -1.4142135623730905 0 0 1.4142135623730745 1.4142135623730936 2.0)
gosh>
2010/11/25 02:29:26 PST
ヤコビ法固有値計算 バグ修正
気になったのでバグ修正。Wilkinsonの公式を理解するのに久々に数式をいじくりまわしたり、行列の計算をしたりで良い勉強になりました。お昼休みにやっつけるつもりが1時間オーバー。遅れた仕事は夜、残って残業して取り返そう。久々に数学、数学って楽しいなぁ。2010/11/24 21:27:24 PST
(define e (list->matrix '((0 1 0)
(1 0 1)
(0 1 0))))
gosh> (matrix-jacobi e)
(-1.414213562373095 0 1.414213562373095)
gosh>
不思議
グラフの固有値、隣接行列のべき乗を計算して楽しんでいます。なんとも不思議です。
完全グラフの固有値はわかりやすいパターンがあるようです。
(define b (list->matrix '((0 1)
(1 0))))
(define k (list->matrix '((0 1 1 1)
(1 0 1 1)
(1 1 0 1)
(1 1 1 0))))
(define j (list->matrix '((0 1 1)
(1 0 1)
(1 1 0))))
gosh> (matrix-jacobi b)
(-1.0 1.0)
gosh> (matrix-jacobi k)
(-1.0 -1.0 -1.0 3.0)
gosh> (matrix-jacobi j)
(-1.0 -1.0 2.0)
gosh>
閉路グラフだともうちょっと複雑。
(define l (list->matrix '((0 1 0)
(1 0 1)
(0 1 0))))
gosh> (matrix-jacobi l)
(-1.41642665362641 -0.18144066652157043 1.59786732014798)
gosh>
scilabでの計算と違ってる。まだバグってる。収束判定がおかしいのかも。
隣接行列のべき乗をとると面白いパターンが
gosh> (mat-print (m-expt k 2)) (3 2 2 2) (2 3 2 2) (2 2 3 2) (2 2 2 3) #<undef> gosh> (mat-print (m-expt k 3)) (6 7 7 7) (7 6 7 7) (7 7 6 7) (7 7 7 6) #<undef> gosh> (mat-print (m-expt l 2)) (1 0 1) (0 2 0) (1 0 1) #<undef> gosh> (mat-print (m-expt l 3)) (0 2 0) (2 0 2) (0 2 0) #<undef> gosh>
2010/11/24 06:29:42 PST
ヤコビ法による固有値計算 バグ
隣接行列の固有値を計算しようとずっと前に書いた自作ライブラリを起動してみました。 対角要素が０の行列の固有値を計算してみると、あれ？おかしい。手計算してみると やはり間違っている。コードのコメントをみるとウィルキンソンの方法により対角要素を 計算しているとあるのだけど、何をしていたのか思い出せない。う～ん、グラフの前にまず この関数を正しく直さないと。まあ、趣味のプログラムだから誰に迷惑をかけるでもなし、いいのだけど。2010/11/22 18:24:02 PST
一応、直して今度はいいみたい。
(define a (list->matrix '((0 0 1)
(0 0 1)
(1 1 0))))
gosh> (matrix-jacobi a)
(-1.414213562373095 0 1.414213562373095)
gosh>
(define k (list->matrix '((0 1 1 1)
(1 0 1 1)
(1 1 0 1)
(1 1 1 0))))
gosh> (matrix-jacobi k)
(-1.0 -0.9999999999999999 -0.9999999999999998 2.9999999999999996)
gosh>
計算誤差を含むけれど、こんなもんかなぁ。
誤差はイプシロンの範囲内なら補正して整数へ変換。これでグラフの固有値の和が0になる定理を確認できる。
gosh> (matrix-jacobi k) (-1.0 -1.0 -1.0 3.0) gosh>
探索のための実験コード
マトリクス会計では資金の流れを追うのに勘定科目から勘定科目に至る経路を知る機能を入れようとしています。会計の場合には単純で入金なら 売上げ->売掛金->受取手形->現金 程度です。複数のルートといっても限られています。が、理解のために実験コードを書いてみました。ウインストンの古いLisp本を参考にしています。待ち行列を使ったものです。2010/11/22 15:27:42 PST コードには無駄な部分、冗長な部分がありますが、ウインストンもそれは承知の上で理解のために敢えてやっているようですのでツッコミはご容赦を。
;;Gaucheで属性リストを扱えるように以下を追加。
(define prop (make-hash-table 'equal?))
(define (putprop name key val)
(hash-table-put! prop (list name key) val))
(define (getprop name key)
(hash-table-get prop (list name key) '()))
(define-macro (defprop name val prop)
`(putprop ',name ',prop ',val))
;;ウインストンのLisp本の例題より「探索を含む例題」
(putprop 's 'neighbors '(l o))
(putprop 'l 'neighbors '(m f))
(putprop 'm 'neighbors '(n))
(putprop 'n 'neighbors '(f))
(putprop 'o 'neighbors '(p q))
(putprop 'p 'neighbors '(f))
(putprop 'q 'neighbors '(f))
(putprop 'o 'neighbors '(s))
;;グラフの連結情報を得るのに隣接行列のget-neighborsを使ってるが
;;(getprop (car path) 'neighbors)にすれば属性リストのグラフを使える。
(define (old-expand path)
(map (lambda (e) (cons e path))
(get-neighbors (car path))))
;;ネットワークである場合にループを回避するためのもの。
(define (new-expand new-paths)
(cond ((null? new-paths) '())
((member (caar new-paths)
(cdar new-paths))
(new-expand (cdr new-paths)))
(else (cons (car new-paths)
(new-expand (cdr new-paths))))))
;;縦型探索
;;オリジナルはprogで書かれている。再帰に書き直した。
;;ウインストンのオリジナルは解がひとつ見つかった時点で探索を終了しているが
;;こちらは強制バックトラックしており全ルートを探索している。
(define (depth start finish)
(define (iter queue)
(cond ((null? queue) #f)
((equal? finish (caar queue))
(print (reverse (car queue)))
(iter (cdr queue)))
(else (iter (append (new-expand (old-expand (car queue))) (cdr queue))))))
(iter (list (list start))))
;;データを隣接行列にしてみる
;;有向グラフなので行からみた列が接続先
;; s l m n o p q f
;;s 0 1 0 0 1 0 0 0
;;l 0 0 1 0 0 0 0 1
;;m 0 0 0 1 0 0 0 0
;;n 0 0 1 0 0 0 0 1
;;o 0 0 0 0 0 1 1 0
;;p 0 0 0 0 0 0 0 1
;;q 0 0 0 0 0 0 0 1
;;f 0 0 0 0 0 0 0 0
;;自作行列ライブラリ
(use math.matrix)
(define mat
(list->matrix
'((0 1 0 0 1 0 0 0)
(0 0 1 0 0 0 0 1)
(0 0 0 1 0 0 0 0)
(0 0 1 0 0 0 0 1)
(0 0 0 0 0 1 1 0)
(0 0 0 0 0 0 0 1)
(0 0 0 0 0 0 0 1)
(0 0 0 0 0 0 0 0))))
;;隣接行列をシンボルで扱うためのもの。
;;頂点を自然数で表すならば以下は不用。
(define (num->sym n)
(case n
((1) 's)
((2) 'l)
((3) 'm)
((4) 'n)
((5) 'o)
((6) 'p)
((7) 'q)
((8) 'f)))
(define (sym->num sym)
(- (+ (matrix-row mat) 1) (length (member sym '(s l m n o p q f)))))
(define (get-neighbors sym)
(define i (sym->num sym))
(define (iter j)
(cond ((<= j 1) '())
((not (zero? (matrix-ref mat i j))) (cons j (iter (- j 1))))
(else (iter (- j 1)))))
(map num->sym (iter (matrix-row mat))))
gosh> (depth 's 'f)
(s o q f)
(s o p f)
(s l f)
(s l m n f)
#f
gosh>
最近、やっていることは
あちこちと気が多いというか、根気がないというか、食い散らかしてます。 化学の構造式をやっていたら飽きてきて、プロジェクトオイラーを解いていたり。 中途半端にしていたマトリクス会計からみでグラフ理論、隣接行列の固有値を勉強してみたり。 そうそう、対称行列に使えるヤコビ法の固有値、固有ベクトル計算は既にライブラリにしてたことを 思い出しました。隣接行列は対称行列なのでこれが使えるじゃん。ああ、数式処理は相変わらず中途状態。2010/11/22 07:35:26 PST
やっているのはこんなこと
まだよく考えやコードがまとまったわけではないのですが、今やっていること（というか遊んでいること）は高校レベルの化学で、こんなことです。
;;ブタン
(define butane '(C H H H (C H H (C H H (C H H H)))))
;;イソブタン
(define isobutane '(C H (C H H H) (C H H H) (C H H H)))
;;カルボキシル基
(define carboxyl-group '(C (=2 O) (O H)))
;;メチル基
(define methyl-group '(C H H H))
;;ヒドロキシ基
(define hydroxy-group '(O H))
;;アセトアルデヒド
(define acetaldehyde '(C (C H H H) (=2 O) H))
;;メタノール
(define methanol '(C H H H (O H)))
(define (atomic-weight element)
(cond ((eq? element 'H) 1.0079)
((eq? element 'C) 12.011)
((eq? element 'N) 14.0067)
((eq? element 'O) 15.9994)
((eq? element 'Na) 22.9898)
((eq? element 'S) 32.06)
((eq? element 'Cl) 35.453)
(else (error "This element not handle:" element))))
;;分子量を構造式から求める
(define (molecular-weight fmla)
(cond ((null? fmla) 0)
((and (symbol? fmla)(or (eq? fmla '=2)(eq? fmla '=3))) 0)
((symbol? fmla) (atomic-weight fmla))
(else (+ (molecular-weight (car fmla))
(molecular-weight (cdr fmla))))))
分子構造には再帰的な要素があると思います。異性体の判定をしたり、アルデヒド基をリストで表された構造式から探したりするコードを模索しています。構造式に再帰があってリストでうまく表せたり再帰で計算できたらハッピーじゃん、というノリです。2010/10/31 07:07:34 PDT
初等整数論から飛んでイスタンブールで、有機化学の問題に取り組んでいます。といっても高校レベル。異性体、同一物質の判定という初歩の初歩のところです。構造式は隣接行列じゃないとダメかと思っていたらリストでもいけますね。二重結合、三重結合を工夫しつつリストで表された有機物質をGaucheを使って弄繰り回してイタズラをしています。化学の復習（何十年ぶり）,お勉強にもなるので一石二鳥です。
エタノールからアセトアルデヒド、酢酸への過程を理解して二日酔いを回避したい。（笑） 2010/10/25 07:31:51 PDT
どうも調子に乗ってあれこれと書き過ぎたかもしれません。書いたものはすべて過去ログに移転しました。コード、数学の誤った記述がインターネット上のごみとして検索阻害要因になるといけません。以後の思いつきのメモはmixiにて公開いたします。私は和製チャールズ・バベージ 2010/10/13 17:01:35 PDT