sage_example の変更点

http://www.pwv.co.jp/~take/TakeWiki/index.php?sage_example

[ トップ ]   [ 編集 | 差分 | バックアップ | 添付 | リロード ]   [ 新規 | 一覧 | 単語検索 | 最終更新 | ヘルプ ]

• 追加された行はこの色です。
• 削除された行はこの色です。
• sage_example へ行く。
• sage_example の差分を削除

---

** はじめに [#w5c1ad9c]
Sageは、Mathematicaのような数式処理を行うオープンソースのソフトウェアです。Sageの歴史はまだ浅く、ウィリアム・スタイン 氏(William Stein)によって2005年2月に開発がスタートし、2006年2月のUCSD SAGE Days 1でSage 1.0が公開され、最新のバージョンは4.7.2です（注１）。

Sageの特徴を挙げると、次のようなものが挙げられます。
- ブラウザーで使える
- 柔軟な図化機能
- Maxima, R等の既存ツールとの連携


** ブラウザーで気軽に使える [#gcf2b32b]
Sageの最大の特徴は、 FirefoxやInternet Explorer等のブラウザーからSage Notebook Serverにアクセスして、気軽に数式処理を実行することが出来ることです（注２）。

ノートブックは、Sageでの一連の計算を記録したノートであり、計算に関する説明文を挿入したり、値を変更して再計算することができます。

*** ノートブックの作成 [#cc63fd06]
Sageのノートブックを体験するには、Sageの開発サイトでアカウントを作成し、ノートブックを作成するのが最も簡単な方法です（注３）。

ログインが完了すると図１のようなノートブック画面になります。

&ref(fstLogin.png);

図1　初期ログイン時のノートブック画面

*** ワークシートの作成とセルの操作 [#od2d0251]
ノートブック画面でNew Worksheetをクリックすると新しいワークシートが作成されます。

ワークシートで式を評価するには、セルと呼ばれるテキストエリアを利用します。

セルの基本操作は、以下のように行います。
- セルの評価：セルに記述した式を評価するには、シフトキーとリターンキーを同時に押す（shift-returnと記す）方法とevaluateをクリックする方法があります。
- セルの追加：セルの上下にマウスを移動すると青い帯が表示されます。この青い帯をクリックするとセルが追加されます。
- セルの削除：セル内のテキストをすべて削除し、もう一度バックスペースキーを押すとセルが削除されます。

*** ヘルプ [#ucc8af32]
Sageの関数の使い方は、関数名のあとに?を付けてshift-returnでヘルプ画面が表示されます（図２参照）。

#pre{{
abs?
}}

&ref(abs_help.png);

図２　abs関数のヘルプ画面

*** 補完機能 [#z39d60ad]
関数名やPythonの変数名は、タブキーで補完することができます。

例として、facを入力した後にタブキーを押すと、自動的にfactorまで入力が補完され、factorとfactorialの２つの候補が表示されています（図３）。

&ref(fac_selection.png);

図３　facに対する補完結果

*** 説明分の挿入 [#a3c2ecd3]
ワークシートには、説明文や図、数式を入力するためのツールが提供されています。セルの上下にマウスを移動すると青い帯が表示されますので、シフトキーとクリックを同時に押すと図４のようなエディタ画面が表示されます。

説明文には、Latex形式で数式が書けるので、計算手法の解説などが簡単に挿入でき、ワークシートの活用に役立ちます。

&ref(editor_screen.png);

図４　エディタ画面の編集例

*** Python文法 [#xd09a4e2]
Sageでは、Pythonのインタプリタを使用しているため、Pythonの文法で記述します。Pythonでは#以下がコメント文となり、セミコロン;で複数の式を1行にまとめて記述することができます。

Sage固有の使い方としては、変数への代入結果を表示するためにセミコロンの後に変数名を追加します（図５の第1セル）。
Sageは数式処理システムですので、結果は数値ではなく数式を評価した値となります（図５の第2セル）。
数値として表示するには、N関数を使用します（図５の第3セル）

&ref(eval.png);

図５　Sage固有の使い方

*** 変数と関数の定義 [#f7318bde]
Sageの数式で変数として認識させるには、変数をvar関数で宣言しなくてはなりません。変数の宣言はvarの引数に文字列で使用する変数名をスペースで区切って行います（図６の第1セル）。

Sageで関数を定義するには、以下の3つの方法があります。
- シンボリック関数：f(x) = 式の形式で定義します。
- Python関数：Pythonのdefを使って関数を定義します。
- lambda関数：Pythonのlambda式を使って関数を定義します。

例として3乗根をシンボリック関数、Python関数、lamba関数で定義し、その結果を表示してみます（図６の第3セル以降）。

&ref(var_func.png);

図６　Sage変数と関数の定義例

*** よく使われる定数 [#r4c4ec34]
Sageでよく使われる定数を表１に示します。

表１　Sageでよく使われる定数
| 円周率 | \( \pi \) | pi |
| 自然対数の底 | \( e \) | e |
| 虚数単位 | \( I \) | I |
| 無限大 | \( \infty\) | oo |

*** グラフの重ね書きもが簡単 [#dc22b119]
Sageのとても便利な機能にグラフの重ね書きがあります。

例として、Sin曲線＋ノイズのデータと元のSin曲線を重ね書きする処理を図７（第1セル）に示します。
- Sin曲線のブロット結果を変数sin_pltに代入
- データのリストプロット結果を変数data_pltに代入
- sin_pltとdta_pltの和に対して、showメッセージを送って重ね書きを実行


&ref(plot.png);

図７　Sin曲線＋ノイズのデータ（100個）と元のSin曲線の重ね書き例

*** インタラクティブな処理 [#n226d5a1]
様々な値に対する結果をインタラクティブに提供するために、@interactコマンドが用意されており、
@interactに続いて、関数定義def _(引数)を記述することで利用可能です。
次のようなものの中からユーザが指定することができます。
- 引数にslider(リスト)またはリストを代入するとスライドバー
- 引数にselector(リスト)を代入するとプルダウンメニュー
- 引数に値を代入するとテキストフィールド


図８にSin曲線のフィッティングにおいて、多項式の次数Mをスライドバーの操作で０から９まで選択可能にした例を示します
（選択可能な値をリストで指定しています）。

&ref(iteractive.png);

図８　Sin曲線のフィッティング問題で多項式の次数Mをスライドバーで与えた例

*** サンプルワークシートのダウンロードと再利用の方法 [#f2d08d73]
本稿で紹介したワークシートや参考文献[1]で使用したワークシート（注4）は、自由にダウンロードできます。

各公開ワークシートを表示し、先頭にあるDownloadを押すと、ダウンロードが開始します。
ダウンロードしたワークシートは、ユーザのノートブックでUploadをクリックし、
ダウンロードしたファイルを指定することで、簡単にアップロードできます。
自分のノートブックでワークシートの計算を再現したり、値を変えて評価することでができますので、是非試してみて下さい。


** 数式処理システムらしい解法 [#n60e8628]
共役勾配法を例にSageの数式処理システムらしい解法を紹介します。
Sageの数式処理機能とPythonの記述力を合わせると、とてもスマートに共役勾配法で2次関数の極値を求めることができます。

以下のような２次形式の関数を考えます。
$$
f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - 6 x_1 - 7 x_2 　（式１）
$$

初期値\( x_0 \)が反復計算により極値に達するためには、
勾配▽fからある程度接線方向tにずれた共役勾配\( d_n \)方向に進む必要があります。ここでは、dの初期値をd0=−∇f(x0) とします。

$$
d_n = - \nabla f(x_n) + \beta_n d_{n-1} , \beta_n = \frac{(\nabla f(x_n))^T \nabla f(x_n)}{(\nabla f(x_{n-1}))^T \nabla f(x_{n-1})} 　　（式２）
$$

\( x_{n+1} \) は刻み値 \( \alpha_n \)  による次の反復公式で更新されます。ここでHは、f(x)のヘッセ行列です。
$$
x_{n+1} = x_n + \alpha_n d_n, \alpha_n = - \frac{d_n^T \nabla f(x_n)}{d_n^T H d_n} 　（式３）
$$


*** Sageでの実装 [#v43079dd]
このように定義される共役勾配法をSageを使って実装していきます。
ベクトルvを変数x1, x2で定義し、ベクトルvを使って関数fを定義します。
#pre{{
# 変数定義
vars = var('x1 x2')
v = vector([x1, x2])
}}

#pre{{
# fを定義
def f(v):
    return 3/2 * v[0]^2 + v[0]*v[1] + v[1]^2 - 6*v[0] - 7*v[1]
}}

次に、▽ｆを計算する関数nabla_fを定義します。ここでは、あらかじめ関数fの各変数で
偏微分した結果をdfsに保存した上で、その結果に利用して引数ベクトルvxの値を代入した結果を
返すように定義しています。
#pre{{
# fを偏微分したリスト
dfs = [diff(f(v), x_i) for x_i in v]
# ▽fを定義　(dfsにvxの要素の値を適応した結果を返す)
def nabla_f(vx):
    # ベクトルvxの各要素の値をvの要素に対応づける
    s = dict(zip(v, vx))
    # ベクトルの各要素の偏微分の結果にsを適応させる
    return vector([df.subs(s) for df in dfs])
}}

ヘッセ行列もSageの数式機能とPythonのリスト内包を使えば、簡単に求めることができます。
最初にdiffコマンドでf(v)を変数x_i, x_jでそれぞれ偏微分した2重リストを作成し、matrix関数で行列型に変換します。
ここでは、変換結果をきれいな行列形式で出力して確認するために、jsmath(H)で画面に表示しています。
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for x_j in v])
print jsmath(H)

}}

&ref(h.png);

\( \alpha_n \) の定義も式３の通りです。

#pre{{
# α_nの定義
def alpha_n(x, d):
    return -d.dot_product(nabla_f(x)) / (d * H * d)

}}

すべての準備が整ったので、共役勾配法を使って極値を計算してみます。
#pre{{
eps = 0.001
x0 = vector([2, 1])
d = - nabla_f(vx=x0)
x = x0
k = 1
while (true):
    o_nabla_f_sqr = nabla_f(x).dot_product(nabla_f(x))
    o_x = x
    x += alpha_n(x, d)*d
    if ((x - o_x).norm() < eps):
        break
    beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f_sqr
    d = -nabla_f(x) + beta*d
    if (d.norm() == 0):    # 0除算への対応
        break
    k += 1
print "x=", x
print "k=", k

}}
結果は、次のようになります。
#pre{{
x= (1, 3)
k= 2
}}
比較のため、Sageの最適化機能で同じ問題を解かせますと、これらの結果は一致します。
#pre{{
# 同様の処理をsageの機能を使って計算してみる
g = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2
minimize(g, [2, 1], algorithm="cg")
}}
#pre{{
       Optimization terminated successfully.
         Current function value: -13.500000
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
(1.0, 3.0)
}}

求まった解をプロットすると、図９のようになります。
#pre{{
# 関数と解をプロット
p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4])
pt = point([1, 3, f(x)], color='red')
(p3d+pt).show()
}}

&ref(fig3.png);

図９　関数fと求めた解をプロットした図

このように数式処理システムの特徴を活かすようにSageでプログラムを記述すると、
数式での表現と実際の計算手続きがきれいに対応するように記述することが可能です。

** 既存ツールとの連携 [#d940443a]
Sageではすべての処理をSage単独で行うのではなく、既存のツールと連携するためのインタフェースを提供しています。

例として、Rでの主成分分析の結果をSageに渡す方法を紹介します（注５）。

Oil Flowのデータを使って主成分分析をします。データは、あらかじめワークシートにアップロードしておきます。
Rとのインタフェース関数rを使ってRのコマンドを記述します。

データを読み込み、主成分分析をresult変数に代入します。
#pre{{
# Rでデータを読み込みPCAを計算
fileName = DATA + 'DataTrn.txt'
oilflow = r("oilflow <- read.table('%s')" %fileName) result = r("result <- prcomp(oilflow)")
}}

RのグラフをSageで表示するには、一度ファイルに出力する必要があります（図１０参照）。
出力された画像ファイルを表示するために、html関数でpca.pdfを表示するHTMLを記述していいます。

#pre{{
#　ラベルの読み込み
fileName = DATA + 'DataTrnLbls.txt'
labels = r("oilflow.labels <- read.table('%s')" %fileName)
# プロットファイル名の設定
filename = DATA+'pca.pdf'
r.pdf(file='"%s"' %filename)
# 結果のプロット
r("col <- colSums(t(oilflow.labels) * c(4,3,2))")
r("pch <- colSums(t(oilflow.labels) * c(3,1,4))")
r("plot(result$x[,1:2], col=col, pch=pch, xlim=c(-3,3), ylim=c(-3,3))")
r.dev_off()
# 式を変えたときにはブラウザーで再読込必要
html('<img src="pca.pdf">')
}}

&ref(fig4.png);

図１０　主成分分析の結果をSageで表示した図

*** RのデータをSageに渡す [#c5daceb5]
RのデータをSageに渡しすためにsageobj関数と_save_メソッドが提供されています。

例として、Rで読み込んだデータセットをSageの3次元プロットを使ってプロットしてみます（図１１参照）。
#pre{{
# RのデータセットをSageの形式に変換すると
# 'DATA'ディクショナリに列単位でV1, V2のようにセットされる 
lb = sageobj(labels)
# これをzipでまとめて使う 
lbs = zip(lb['DATA']['V1'],lb['DATA']['V2'],lb['DATA']['V3'])
# point3dでプロット
N = len(lbs)
plt = Graphics()
for n in range(N):
    [x, y, z] = rs[n]
    if lbs[n][0] == 1:
        plt += point3d([x, y, z], rgbcolor='blue')
    elif lbs[n][1] == 1:
        plt += point3d([x, y, z], rgbcolor='green')
    else:
        plt += point3d([x, y, z], rgbcolor='red')    
plt.show()
}}

&ref(fig5.png);

図１１　Rで読み込んだデータをSageでプロットした図

** Sageの今後 [#g16dd496]
このように柔軟で便利な機能をもつSageですが、日本国内ではまだまだ普及していません。本稿を通じて、一人でも多くの方にSageを使って頂ければ幸いです。今後のユーザ会などの支援団体の創設や解説書の出版が待たれます。

沼田泰英氏がSageのレファレンスカードを日本語に翻訳して公開されています（注６）。印刷して傍らに置いておくと便利です。

最後に、本稿の査読にご協力頂きました(株)構造計画研究所の鈴木 由宇氏、本稿のLaTexへの変換にご協力頂きました神戸大学の長坂 耕作准教授に感謝申し上げます。

** 参考文献 [#k0a21bc9]
- 竹本　浩, Software Design 2010/06, OSS数式処理システム「Sage」, p89
- 竹本　浩, Software Design 2010/06, OSS数式処理システム「Sage」, p89-94

** 脚注 [#ped2d672]
- 注１：William Steinの2007年UW CSE コロキアムの資料から引用。
- 注２：Linux, Mac OSX, Windowsにダウンロードして使うこともできます。
- 注３：Sage公開サーバのURL：http://www.sagenb.org/、筆者のサーバにも公開サーバhttp://www.pwv.co.jp:8000/を置いています。
- 注４：http://www.pwv.co.jp:8000/home/pub/16/にワークシートを公開しています。
- 注５：http://www.pwv.co.jp:8000/home/pub/14/にワークシートを公開しています。
- 注６：http://www.stat.t.u-tokyo.ac.jp/~numata/nora/sage-doc/

       

Site admin: anonymous

PukiWiki 1.4.7 Copyright © 2001-2006 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji. Powered by PHP 5.1.6. HTML convert time: 0.003 sec.