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にアクセスして、気軽に数式処理を実行することが出来ることです（注２）。

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

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

初期画面には、1個のセルが表示され、ここにpython形式で数式を入力し、evaluateボタンで結果が表示されます。

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

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

&re(editor_screen.png);


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

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

例として、Sin曲線のフィッティング用データと元のSin曲線を重ね書きする処理はリスト１のようになります。
- 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() 
}}
リスト１　Sin曲線のフィッティング用データと元のSin曲線の重ね書き例

*** インタラクティブな処理 [#n226d5a1]
様々な値に対する結果をインタラクティブに提供するために、@interactコマンドが用意されています。

@interactに続いて、関数定義def _(引数)を記述します。
- 引数にslider(リスト)またはリストを代入するとスライドバー
- 引数にselector(リスト)を代入するとプルダウンメニュー
- 引数に値を代入するとテキストフィールド

からユーザが値を指定することができます。

図２にSin曲線のフィッティング問題で多項式の次数Mに０から９までの値のリストを与えた例を示します。

&ref(iteractive.png);


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

** 数式処理システムらしい解法 [#n60e8628]
共役勾配法を例にSageの数式処理システムらしい解法を紹介します。
Sageの数式処理機能とPythonの記述力を合わせると、とてもスマートに共役勾配法で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]
}}

次に、▽ｆを計算します。
あらかじめ関数ｆの各変数での偏微分した結果を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)

}}

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

}}

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


** 既存ツールとの連携 [#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で表示するには、一度ファイルに出力する必要があります（図４参照）。
#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">') 

}}


図４　主成分分析の結果を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() 
}}


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

** Sageの今後 [#g16dd496]
このように柔軟で便利な機能をもつSageですが、日本国内ではまだまだ普及していません。本稿を通じて、一人でも多くの方にSageを使って頂ければ幸いです。今後のユーザ会などの支援団体の創設や解説書の出版が待たれます。
沼田泰英氏がSageのレファレンスカードを日本語に翻訳して公開されています（注６）。印刷して傍らに置いておくと便利です。


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