![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2011/04/22からのアクセス回数 6025
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www.pwv.co.jp:8000/home/pub/21/
また、Sageのサーバを公開しているサイト（http://www.pwv.co.jp:8000/、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
author>Hiroshi TAKEMOTO</author> (<email>take@pwv.co.jp</email>)
インターフェース2011年12号にMATLABの特集があったので、Sageでoctaveを動かしてみようと思います。
sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。
a href="http://www.pwv.co.jp/~take/TakeWiki/index.php?octave%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB">「octaveのインストール」参照</a>
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))
#pre{{
}}
sageからoctaveの機能を呼び出すには、octave, octave.evalインターフェースを使用します。
以下に、octave.eval, octaveの使用例を示します。octave.evalでは結果が出力されますので、これを避けるためにjunk変数に代入しています。
sageへの入力:
# octaveで入力した文字列を評価させる場合
octave.eval('2+2')
'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関数で変換します。
<class 'sage.interfaces.octave.OctaveElement'> <type 'sage.rings.real_mpfr.RealNumber'> 4 12
octaveの使い方をコンパクトに説明したサイトが<a href="http://www.inaba-lab.org/modules/bwiki/index.php?Octave%C6%FE%CC%E7">東海大学 稲葉研究室</a>にあります。
sageへの入力:
# octaveの使い方については、東海大学 稲葉研究室のサイトが簡潔で分かりやすい # http://www.inaba-lab.org/modules/bwiki/index.php?Octave%C6%FE%CC%E7
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
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
4 4 12 15 24 35 49 64 81
octave(MATLAB)の便利な特徴に左演算子\があります。
以下の様なベクトル計算でX=[a; b]としてXを求めるためにちょうど両辺を左からXで割る ような演算（左除算演算子）\を使うと、 $$ B = A X $$ から $$ A \verb|\| 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
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
1.5 -0.666667 3.83333 5.33333 6.83333
先の計算で求められた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()
先の例では、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
[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のワークシートに表示するには、 一端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')
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への入力:
# グレイスケールに変換
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">')
皆様のご意見、ご希望をお待ちしております。
