- 追加された行はこの色です。
- 削除された行はこの色です。
sageは、mathematicaのような数式処理を行うオープンソースのソフトウェアです。sageの歴史はまだ浅く、ウィリアム・スタイン 氏(William Stein)によって2005年2月に開発がスタートし、2006年2月のUCSD SAGE Days 1でSage 1.0が公開され、最新のバージョンは4.7.2です(注1)。
投稿用原稿のため、掲載先の許可がでるまで削除します。ご了承ください。
Sageの特徴を挙げると、
- ブラウザーで使える
- 柔軟な図化機能
- Maxima, R等の既存ツールとの連携
があります。
** ブラウザーで気軽に使える [#gcf2b32b]
Sageの最大の特徴は、 FirefoxやInternet Explorer等のブラウザーからSage Notebook Serverにアクセスして、気軽に数式処理を実行することが出来ることです(注2)。
Notebookは、Sageでの一連の計算を記録したノートであり、計算に関する説明文を挿入したり、値を変更して再計算することができます。
*** ノートブックの作成 [#cc63fd06]
Sageのノートブックを体験するには、Sageの開発サイトでアカウントを作成し、ワークシートを作成するのが最も簡単な方法です(注3)。
初期画面には、1個のセルが表示され、ここにpython形式で数式を入力し、evaluateボタンで結果が表示されます。
ノートブックには、説明文や図、数式を入力するためのツールが提供されています。セルの上下にマウスを移動すると青い帯が表示されますので、シフトキーとクリックを同時に押すと図1のようなエディタ画面が表示されます。
説明文には、Latex形式で数式が書けるので、計算手法の解説などが簡単に挿入でき、ノートブックの活用に役立ちます。
&re(editor_screen.png);
図1 エディタ画面の編集例
*** グラフの重ね書きもが簡単 [#dc22b119]
Sageのとても便利な機能にグラフの重ね書きがあります(注4)。
例として、Sin曲線のフィッティング用データと元のSin曲線を重ね書きする処理はリスト1のようになります。
- Sin曲線のブロット結果を変数sin_pltに代入
- データのリストプロット結果を変数data_pltに代入
- sin_pltとdta_pltの和に対して、showメッセージを送って重ね書きを実行
#pre{{
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green');
data_plt = list_plot(zip(X, t)); data_plt (data_plt + sin_plt).show()
}}
リスト1 Sin曲線のフィッティング用データと元のSin曲線の重ね書き例
*** インタラクティブな処理 [#n226d5a1]
様々な値に対する結果をインタラクティブに提供するために、@interactコマンドが用意されています。
@interactに続いて、関数定義def _(引数)を記述します。
- 引数にslider(リスト)またはリストを代入するとスライドバー
- 引数にselector(リスト)を代入するとプルダウンメニュー
- 引数に値を代入するとテキストフィールド
からユーザが値を指定することができます。
図2にSin曲線のフィッティング問題で多項式の次数Mに0から9までの値のリストを与えた例を示します。
&ref(iteractive.png);
図2 Sin曲線のフィッティング問題で多項式の次数Mをスライダーバーで与えた例
** 数式処理システムらしい解法 [#n60e8628]
共役勾配法を例にSageの数式処理システムらしい解法を紹介します。
Sageの数式処理機能とPythonの記述力を合わせると、とてもスマートに共役勾配法で2次関数の極値を求めることができます。
以下のような2次形式の関数を考えます。
極値に達するには、勾配▽fからある程度接線方向tにずれた共役勾配dn方向に進みます。
βnは
となり、dの初期値はd0=−∇f(x0) から始めます。
xは刻み値α、
を使って次式で更新します。 ここでHは、f(x)のヘッセ行列です。
ベクトルvを変数x1, x2で定義し、これを使って関数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]
}}
次に、▽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の数式機能を使えば、簡単に求めることができます。
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for x_j in v])
print jsmath(H)
}}
α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)
}}
求まった解をプロットすると、図3のようになります。
#pre{{
# 関数と解をプロット
p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4])
pt = point([1, 3, f(x)], color='red')
(p3d+pt).show()
}}
図3 関数fと求めた解をプロットした図
** 既存ツールとの連携 [#d940443a]
Sageではすべての処理をSage単独で行うのではなく、既存のツールと連携するためのインタフェースを提供しています。
例として、Rでの主成分分析の結果をSageに渡す方法を紹介します(注5)。
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で表示するには、一度ファイルに出力する必要があります(図4参照)。
#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">')
}}
図4 主成分分析の結果をSageで表示した図
*** RのデータをSageに渡す [#c5daceb5]
Rのデータをsageに渡しすためにsageobj関数と_save_メソッドが提供されています。
例として、Rで読み込んだデータセットをSageの3次元プロットを使ってプロットしてみます(図4参照)。
#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()
}}
図5 Rで読み込んだデータをSageでプロットした図
** Sageの今後 [#g16dd496]
このように柔軟で便利な機能をもつSageですが、日本国内ではまだまだ普及していません。本稿を通じて、一人でも多くの方にSageを使って頂ければ幸いです。今後のユーザ会などの支援団体の創設や解説書の出版が待たれます。
沼田泰英氏がSageのレファレンスカードを日本語に翻訳して公開されています(注6)。印刷して傍らに置いておくと便利です。
** 脚注 [#ped2d672]
- 注1:William Steinの2007年UW CSE コロキアムの資料から引用。
- 注2:Linux, Mac OSX, Windowsにダウンロードして使うこともできます。
- 注3:Sage公開サーバのURL:http://www.sagenb.org/、筆者のサーバにも公開サーバ(htttp://www.pwv.co.jp:8000/)を置いています。
- 注4:http://www.pwv.co.jp:8000/home/pub/20/にワークシートを公開しています。
- 注5:http://www.pwv.co.jp:8000/home/pub/14/にワークシートを公開しています。
- 注6:http://www.stat.t.u-tokyo.ac.jp/~numata/nora/sage-doc/