![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
|
2012/03/25からのアクセス回数 5770 ここで紹介したSageワークシートは、以下のURLからダウンロードできます。 http://www.pwv.co.jp:8000/home/pub/21/ また、Sageのサーバを公開しているサイト（http://www15191ue.sakura.ne.jp:8000/、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。 SageでOctaveを使う †インターフェース2011年12号にMATLABの特集があったので、Sageでoctaveを動かしてみようと思います。 sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。 octaveのインストール参照 事前設定項目 †sageからoctaveを使用する場合、あらかじめ以下の設定をすると便利です。
sageとoctaveがローカルのマシンにインストールされている場合には、GNUTERMの設定は不要です。 その場合には、sageでoctaveのプロットを行うとローカルの画面にプロット結果が表示されますので、 逐次確認を行いながらワークシートを作成することができます。 sageへの入力: # octaveのグラフ出力ターミナルをpngにセット
junk = octave.eval('putenv("GNUTERM", "png")')
# sageとoctaveがインストールされてマシンで、octaveのグラフを確認しながら
# 実行する場合には、上記の設定をコメントアウトしてください。
sageへの入力: # octaveのグラフをファイルに保存し、表示する関数を定義
def saveAndShowPlot(filename, fac=0.75):
output = DATA + filename
octave.eval("print -dpng '%s'"%output)
width = int(640*fac)
html('<img src="%s" width="%spx">'%(filename, width))
sage-octaveインターフェース †sageからoctaveの機能を呼び出すには、octave, octave.evalインターフェースを使用します。
以下に、octave.eval, octaveの使用例を示します。octave.evalでは結果が出力されますので、これを避けるためにjunk変数に代入しています。 sageへの入力: # octaveで入力した文字列を評価させる場合
octave.eval('2+2')
sageの出力: 'ans = 4' octaveの戻り値はOctaveElemntとして返されます。OctaveElementに対しては、 算術計算ができるためsageでのループ処理に利用できます。 また、sageで数値と扱うには、N関数を使って一度数値に変換して使用します。 sageへの入力: # octaveの計算結果をsageの変数に代入する場合
sum = octave('2+2')
print type(sum), type(N(sum)), sum, sum*3
# 通常の計算ではそのまま使えますが、数値として扱う場合にはN関数で変換します。
sageの出力: <class 'sage.interfaces.octave.OctaveElement'> <type 'sage.rings.real_mpfr.RealNumber'> 4 12 octaveの参考サイト †octaveの使い方をコンパクトに説明したサイトが 東海大学 稲葉研究室 にあります。 sageへの入力: # octaveの使い方については、東海大学 稲葉研究室のサイトが簡潔で分かりやすい # http://www.inaba-lab.org/modules/bwiki/index.php?Octave%C6%FE%CC%E7 MATLABクローンのoctaveの特徴 †octaveはMATLABのフリーのクローンであり、不完全ながらMATLABの関数と互換性があります。 sageの不得意とするデータ処理はoctaveの得意とするところであり、その豊富なライブラリと 行列計算の特徴を活かして、信号処理、画像処理、制御に多く使われています。 以下ではマトリックスAとBの演算をoctaveで行う場合の例を示します。 $$ A = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 3 & 4 & 5 \\ 7 & 8 & 9 \end{array} \right] $$ $$ B = \left[ \begin{array}{ccc} 4 & 2 & 4 \\ 5 & 6 & 7 \\ 7 & 8 & 9 \end{array} \right] $$ sageへの入力: # octaveは、MATLABのクローンの一つであり、マトリックス計算が得意
# マトリックスの計算
A = octave("A=[1, 2, 3; 3, 4, 5; 7, 8, 9]")
# 区切り記号はカンマの他にスペースも使える
B = octave("B=[4 2 4;5 6 7;7 8 9]")
print A, B, A+B, A*B
sageの出力: 1 2 3 3 4 5 7 8 9 4 2 4 5 6 7 7 8 9 5 4 7 8 10 12 14 16 18 35 38 45 67 70 85 131 134 165 要素演算子 †octaveのマトリックス計算では通常の計算の他に、行列の要素単位での演算を 行うために.（ピリオド）の後に演算子を指定する方法が提供されています。 かけ算*、左除算\、べき乗が^等を使うことができます。 $$ A .* B $$ の結果を以下に示します。 sageへの入力: # 要素毎への計算には.を付ける
C = octave("A.*B"); C
sageの出力: 4 4 12 15 24 35 49 64 81 左除算演算子 †octave(MATLAB)の便利な特徴に左演算子\があります。 以下の様なベクトル計算でX=[a; b]としてXを求めるためにちょうど両辺を左からXで割る ような演算（左除算演算子）\を使うと、 $$ B = A X $$ から $$ A ￥ B = X $$ となり、Xを求めることができます。 例として以下の計算をoctaveで試してみましょう。結果は、X=[1; 1], A*X=[4, 5]となり, Bの値と一致します。 $$ \left[ \begin{array}{c} 4 \\ 5 \end{array} \right] = \left[ \begin{array}{cc} 3 & 1\\ 4 & 1 \end{array} \right] \left[ \begin{array}{c} a \\ b \end{array} \right] $$ sageへの入力: # 左除算演算子\を使うと方程式や最小自乗法を一つの表現で処理できる
# B = A*X となるXを求めるために両辺を左からAで割るイメージ
A = octave("A=[3,1;4,1]")
B = octave("B=[4;5]")
X = octave("X=A\B")
print X, A*X
sageの出力: 1 1 4 5 最小自乗法（フィッチング曲線）を求める †左除算演算子\のすごいところは、そのまま最小自乗法で誤差を最小とするような解を計算してくれることです。 実験の結果、xの値[3; 4; 5]に対して観測値[4; 5; 6]が求まったとする。これを行列で表現すると、 $$ \left[ \begin{array}{c} 4 \\ 5 \\ 7 \end{array} \right] = \left[ \begin{array}{cc} 3 & 1\\ 4 & 1\\ 5 & 7 \end{array} \right] \left[ \begin{array}{c} a \\ b \end{array} \right] $$ となります。これからX= [a; b]を求めると以下の様になります。 sageへの入力: # 4 = a*3 + b
# 5 = a*4 + b
# 7 = a*5 + b
# のデータを得られた場合の誤差をもっとも小さくするようなa, bを求めると
A = octave("A=[3,1;4,1;5,1]")
B = octave("B=[4;5;7]")
X = octave("X=A\B")
print X, A*X
sageの出力: 1.5 -0.666667 3.83333 5.33333 6.83333 求まった結果をsageでグラフに表示 †先の計算で求められたX=[a; b]の値を使って、解の直線y = a*x +bとデータの点をプロットしてみます。 a, bに代入する場合には、N関数を使ってoctave要素から数値に変換します。そして、sageのプロット関数 を使えば簡単にベストフィッティングの解が求まったことが確認できます。 sageへの入力: a = N(X(1)) b = N(X(2)) pts = [[3,4], [4,5], [5,7]] pt_plt = list_plot(pts) line_plt = plot(a*x+b, [x, 1, 5]) (pt_plt+line_plt).show() sageの出力:
sageからマトリックスをoctaveに渡す場合 †先の例では、octaveの結果（ベクトル）から要素を取り出し、sageで結果をプロットしましたが、 逆にsageのマトリックスをoctaveに渡すときには、octave.sage2octave_matrix_string を使用して、一端octaveのマトリックス表現文字列に変換します。 変換されたマトリックス表現文字列をoctaveに渡すことでsageからoctaveへのマトリックス 変換ができます。 sageへの入力: # sageのマトリックスを作成 m = matrix([[1,2,3],[4,5,6],[7,8,9]]) # sageのマトリックスからoctaveのマトリックス指定文字列に変換 om = octave.sage2octave_matrix_string(m) print om # sageのmをoctaveに生成する mm = octave(om) print type(mm), mm sageの出力: [1, 2, 3; 4, 5, 6; 7, 8, 9] <class 'sage.interfaces.octave.OctaveElement'> 1 2 3 4 5 6 7 8 9 octaveのプロットをsageで表示するには †octaveのプロット結果をsageのワークシートに表示するには、 一端png形式の画像ファイルに保存し、HTMLのimgタグを使って 画像を表示します。 こられの処理は冒頭ににありました、saveAndShowPlot関数を 呼び出すことで、簡単に行うことができます。 sageへの入力: # octaveのプロット結果を表示する
octave.eval('x = 0:0.1:10;')
octave.eval('y = sin(x);')
octave.eval('plot(x, y);')
saveAndShowPlot('test.png')
sageの出力:
octaveは信号処理や画像処理が得意 †octaveには、信号処理や画像処理を行うためのライブラリが充実しており、 それらをsageから使うことによってsageの処理の幅が広がります。 画像処理の例として、lena_std.jpgの画像をグレイスケールに変換してみます。 sageへの入力: # octaveは信号処理や画像処理が得意
html('<img src="lena_std.jpg" width="256px">')
# 画像を読み込む
input = DATA+"lena_std.jpg"
junk = octave.eval('rgb = imread("%s");'%input)
sageの出力:
sageへの入力: # グレイスケールに変換
octave.eval('mono = rgb2gray(rgb);')
output = DATA+"lena_gray.png"
octave.eval("imwrite(mono, '%s', 'png')"%output)
html('<img src="lena_gray.png" width="256px">')
sageの出力:
コメント †皆様のご意見、ご希望をお待ちしております。 Tweet |
lena_std.jpg
test.png
sage0.png
lena_gray.png
