sage_example
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
** はじめに [#w5c1ad9c]
Sageは、Mathematicaのような数式処理を行うオープンソースの...
Sageの特徴を挙げると、次のようなものが挙げられます。
- ブラウザーで使える
- 柔軟な図化機能
- Maxima, R等の既存ツールとの連携
** ブラウザーで気軽に使える [#gcf2b32b]
Sageの最大の特徴は、 FirefoxやInternet Explorer等のブラウ...
ノートブックは、Sageでの一連の計算を記録したノートであり...
*** ノートブックの作成 [#cc63fd06]
Sageのノートブックを体験するには、Sageの開発サイトでアカ...
ログインが完了すると図１のようなノートブック画面になりま...
&ref(fstLogin.png);
図1 初期ログイン時のノートブック画面
*** ワークシートの作成とセルの操作 [#od2d0251]
ノートブック画面でNew Worksheetをクリックすると新しいワー...
ワークシートで式を評価するには、セルと呼ばれるテキストエ...
セルの基本操作は、以下のように行います。
- セルの評価：セルに記述した式を評価するには、シフトキー...
- セルの追加：セルの上下にマウスを移動すると青い帯が表示...
- セルの削除：セル内のテキストをすべて削除し、もう一度バ...
*** ヘルプ [#ucc8af32]
Sageの関数の使い方は、関数名のあとに?を付けてshift-return...
#pre{{
abs?
}}
&ref(abs_help.png);
図２ abs関数のヘルプ画面
*** 補完機能 [#z39d60ad]
関数名やPythonの変数名は、タブキーで補完することができま...
例として、facを入力した後にタブキーを押すと、自動的にfact...
&ref(fac_selection.png);
図３ facに対する補完結果
*** 説明分の挿入 [#a3c2ecd3]
ワークシートには、説明文や図、数式を入力するためのツール...
説明文には、Latex形式で数式が書けるので、計算手法の解説な...
&ref(editor_screen.png);
図４ エディタ画面の編集例
*** Python文法 [#xd09a4e2]
Sageでは、Pythonのインタプリタを使用しているため、Python...
Sage固有の使い方としては、変数への代入結果を表示するため...
Sageは数式処理システムですので、結果は数値ではなく数式を...
数値として表示するには、N関数を使用します（図５の第3セル）
&ref(eval.png);
図５ Sage固有の使い方
*** 変数と関数の定義 [#f7318bde]
Sageの数式で変数として認識させるには、変数をvar関数で宣言...
Sageで関数を定義するには、以下の3つの方法があります。
- シンボリック関数：f(x) = 式の形式で定義します。
- Python関数：Pythonのdefを使って関数を定義します。
- lambda関数：Pythonのlambda式を使って関数を定義します。
例として3乗根をシンボリック関数、Python関数、lamba関数で...
&ref(var_func.png);
図６ Sage変数と関数の定義例
*** よく使われる定数 [#r4c4ec34]
Sageでよく使われる定数を表１に示します。
表１ Sageでよく使われる定数
| 円周率 | \( \pi \) | pi |
| 自然対数の底 | \( e \) | e |
| 虚数単位 | \( I \) | I |
| 無限大 | \( \infty\) | oo |
*** グラフの重ね書きもが簡単 [#dc22b119]
Sageのとても便利な機能にグラフの重ね書きがあります。
例として、Sin曲線＋ノイズのデータと元のSin曲線を重ね書き...
- Sin曲線のブロット結果を変数sin_pltに代入
- データのリストプロット結果を変数data_pltに代入
- sin_pltとdta_pltの和に対して、showメッセージを送って重...
&ref(plot.png);
図７ Sin曲線＋ノイズのデータ（100個）と元のSin曲線の重ね...
*** インタラクティブな処理 [#n226d5a1]
様々な値に対する結果をインタラクティブに提供するために、@...
@interactに続いて、関数定義def _(引数)を記述することで利...
次のようなものの中からユーザが指定することができます。
- 引数にslider(リスト)またはリストを代入するとスライドバー
- 引数にselector(リスト)を代入するとプルダウンメニュー
- 引数に値を代入するとテキストフィールド
図８にSin曲線のフィッティングにおいて、多項式の次数Mをス...
（選択可能な値をリストで指定しています）。
&ref(iteractive.png);
図８ Sin曲線のフィッティング問題で多項式の次数Mをスライ...
*** サンプルワークシートのダウンロードと再利用の方法 [#f2...
本稿で紹介したワークシートや参考文献[1]で使用したワークシ...
各公開ワークシートを表示し、先頭にあるDownloadを押すと、...
ダウンロードしたワークシートは、ユーザのノートブックでUpl...
ダウンロードしたファイルを指定することで、簡単にアップロ...
自分のノートブックでワークシートの計算を再現したり、値を...
** 数式処理システムらしい解法 [#n60e8628]
共役勾配法を例にSageの数式処理システムらしい解法を紹介し...
Sageの数式処理機能とPythonの記述力を合わせると、とてもス...
以下のような２次形式の関数を考えます。
$$
f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - 6 x_1 - 7 x_...
$$
初期値\( x_0 \)が反復計算により極値に達するためには、
勾配▽fからある程度接線方向tにずれた共役勾配\( d_n \)方向...
$$
d_n = - \nabla f(x_n) + \beta_n d_{n-1} , \beta_n = \frac...
$$
\( x_{n+1} \) は刻み値 \( \alpha_n \) による次の反復公式...
$$
x_{n+1} = x_n + \alpha_n d_n, \alpha_n = - \frac{d_n^T \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...
}}
次に、▽ｆを計算する関数nabla_fを定義します。ここでは、あ...
偏微分した結果をdfsに保存した上で、その結果に利用して引数...
返すように定義しています。
#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でそれぞれ偏微分し...
ここでは、変換結果をきれいな行列形式で出力して確認するた...
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for ...
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...
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) resu...
}}
RのグラフをSageで表示するには、一度ファイルに出力する必要...
出力された画像ファイルを表示するために、html関数でpca.pdf...
#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), y...
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']['V...
# 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のレファレンスカードを日本語に翻訳して公...
最後に、本稿の査読にご協力頂きました(株)構造計画研究所の...
** 参考文献 [#k0a21bc9]
- 竹本 浩, Software Design 2010/06, OSS数式処理システム...
** 脚注 [#ped2d672]
- 注１：William Steinの2007年UW CSE コロキアムの資料から...
- 注２：Linux, Mac OSX, Windowsにダウンロードして使うこと...
- 注３：Sage公開サーバのURL：http://www.sagenb.org/、筆者...
- 注４：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...
終了行:
** はじめに [#w5c1ad9c]
Sageは、Mathematicaのような数式処理を行うオープンソースの...
Sageの特徴を挙げると、次のようなものが挙げられます。
- ブラウザーで使える
- 柔軟な図化機能
- Maxima, R等の既存ツールとの連携
** ブラウザーで気軽に使える [#gcf2b32b]
Sageの最大の特徴は、 FirefoxやInternet Explorer等のブラウ...
ノートブックは、Sageでの一連の計算を記録したノートであり...
*** ノートブックの作成 [#cc63fd06]
Sageのノートブックを体験するには、Sageの開発サイトでアカ...
ログインが完了すると図１のようなノートブック画面になりま...
&ref(fstLogin.png);
図1 初期ログイン時のノートブック画面
*** ワークシートの作成とセルの操作 [#od2d0251]
ノートブック画面でNew Worksheetをクリックすると新しいワー...
ワークシートで式を評価するには、セルと呼ばれるテキストエ...
セルの基本操作は、以下のように行います。
- セルの評価：セルに記述した式を評価するには、シフトキー...
- セルの追加：セルの上下にマウスを移動すると青い帯が表示...
- セルの削除：セル内のテキストをすべて削除し、もう一度バ...
*** ヘルプ [#ucc8af32]
Sageの関数の使い方は、関数名のあとに?を付けてshift-return...
#pre{{
abs?
}}
&ref(abs_help.png);
図２ abs関数のヘルプ画面
*** 補完機能 [#z39d60ad]
関数名やPythonの変数名は、タブキーで補完することができま...
例として、facを入力した後にタブキーを押すと、自動的にfact...
&ref(fac_selection.png);
図３ facに対する補完結果
*** 説明分の挿入 [#a3c2ecd3]
ワークシートには、説明文や図、数式を入力するためのツール...
説明文には、Latex形式で数式が書けるので、計算手法の解説な...
&ref(editor_screen.png);
図４ エディタ画面の編集例
*** Python文法 [#xd09a4e2]
Sageでは、Pythonのインタプリタを使用しているため、Python...
Sage固有の使い方としては、変数への代入結果を表示するため...
Sageは数式処理システムですので、結果は数値ではなく数式を...
数値として表示するには、N関数を使用します（図５の第3セル）
&ref(eval.png);
図５ Sage固有の使い方
*** 変数と関数の定義 [#f7318bde]
Sageの数式で変数として認識させるには、変数をvar関数で宣言...
Sageで関数を定義するには、以下の3つの方法があります。
- シンボリック関数：f(x) = 式の形式で定義します。
- Python関数：Pythonのdefを使って関数を定義します。
- lambda関数：Pythonのlambda式を使って関数を定義します。
例として3乗根をシンボリック関数、Python関数、lamba関数で...
&ref(var_func.png);
図６ Sage変数と関数の定義例
*** よく使われる定数 [#r4c4ec34]
Sageでよく使われる定数を表１に示します。
表１ Sageでよく使われる定数
| 円周率 | \( \pi \) | pi |
| 自然対数の底 | \( e \) | e |
| 虚数単位 | \( I \) | I |
| 無限大 | \( \infty\) | oo |
*** グラフの重ね書きもが簡単 [#dc22b119]
Sageのとても便利な機能にグラフの重ね書きがあります。
例として、Sin曲線＋ノイズのデータと元のSin曲線を重ね書き...
- Sin曲線のブロット結果を変数sin_pltに代入
- データのリストプロット結果を変数data_pltに代入
- sin_pltとdta_pltの和に対して、showメッセージを送って重...
&ref(plot.png);
図７ Sin曲線＋ノイズのデータ（100個）と元のSin曲線の重ね...
*** インタラクティブな処理 [#n226d5a1]
様々な値に対する結果をインタラクティブに提供するために、@...
@interactに続いて、関数定義def _(引数)を記述することで利...
次のようなものの中からユーザが指定することができます。
- 引数にslider(リスト)またはリストを代入するとスライドバー
- 引数にselector(リスト)を代入するとプルダウンメニュー
- 引数に値を代入するとテキストフィールド
図８にSin曲線のフィッティングにおいて、多項式の次数Mをス...
（選択可能な値をリストで指定しています）。
&ref(iteractive.png);
図８ Sin曲線のフィッティング問題で多項式の次数Mをスライ...
*** サンプルワークシートのダウンロードと再利用の方法 [#f2...
本稿で紹介したワークシートや参考文献[1]で使用したワークシ...
各公開ワークシートを表示し、先頭にあるDownloadを押すと、...
ダウンロードしたワークシートは、ユーザのノートブックでUpl...
ダウンロードしたファイルを指定することで、簡単にアップロ...
自分のノートブックでワークシートの計算を再現したり、値を...
** 数式処理システムらしい解法 [#n60e8628]
共役勾配法を例にSageの数式処理システムらしい解法を紹介し...
Sageの数式処理機能とPythonの記述力を合わせると、とてもス...
以下のような２次形式の関数を考えます。
$$
f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - 6 x_1 - 7 x_...
$$
初期値\( x_0 \)が反復計算により極値に達するためには、
勾配▽fからある程度接線方向tにずれた共役勾配\( d_n \)方向...
$$
d_n = - \nabla f(x_n) + \beta_n d_{n-1} , \beta_n = \frac...
$$
\( x_{n+1} \) は刻み値 \( \alpha_n \) による次の反復公式...
$$
x_{n+1} = x_n + \alpha_n d_n, \alpha_n = - \frac{d_n^T \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...
}}
次に、▽ｆを計算する関数nabla_fを定義します。ここでは、あ...
偏微分した結果をdfsに保存した上で、その結果に利用して引数...
返すように定義しています。
#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でそれぞれ偏微分し...
ここでは、変換結果をきれいな行列形式で出力して確認するた...
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for ...
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...
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) resu...
}}
RのグラフをSageで表示するには、一度ファイルに出力する必要...
出力された画像ファイルを表示するために、html関数でpca.pdf...
#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), y...
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']['V...
# 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のレファレンスカードを日本語に翻訳して公...
最後に、本稿の査読にご協力頂きました(株)構造計画研究所の...
** 参考文献 [#k0a21bc9]
- 竹本 浩, Software Design 2010/06, OSS数式処理システム...
** 脚注 [#ped2d672]
- 注１：William Steinの2007年UW CSE コロキアムの資料から...
- 注２：Linux, Mac OSX, Windowsにダウンロードして使うこと...
- 注３：Sage公開サーバのURL：http://www.sagenb.org/、筆者...
- 注４：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...
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
