Octave_example

4837 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

SageでOctaveを使う

インターフェース2011年12号にMATLABの特集があったので、Sageでoctaveを動かしてみようと思います。

sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。 「octaveのインストール」参照

事前設定項目

sageからoctaveを使用する場合、あらかじめ以下の設定をすると便利です。

  • GNUTERMをpngにセットする
  • グラフ保存と表示を行うsaveAndShowPlot関数を定義する

sageとoctaveがローカルのマシンにインストールされている場合には、GNUTERMの設定は不要です。 その場合には、sageでoctaveのプロットを行うとローカルの画面にプロット結果が表示されますので、 逐次確認を行いながらワークシートを作成することができます。

import time # octaveのグラフ出力ターミナルをpngにセット #junk = octave.eval('putenv("GNUTERM", "png")') # octave3.6.1の場合コメントアウトしてください。 # sageとoctaveがインストールされているマシンで、octaveのグラフを確認しながら # 実行する場合には、上記の設定をコメントアウトしてください。 # octaveのグラフをファイルに保存し、表示する関数を定義（画像のキャッシュをしない版） def saveAndShowPlot(filename, fac=0.75): output = DATA + filename octave.eval("print -dpng '%s'"%output) width = int(640*fac) html('<img src="%s?%s" width="%spx">'%(filename, time.time(), width)) 
       

sage-octaveインターフェース

sageからoctaveの機能を呼び出すには、octave, octave.evalインターフェースを使用します。

  • octave.evalは、文字列で指定した式をoctaveに評価される場合に使用します
  • octaveは、octaveの結果をsageの変数に代入し、octaveの評価を再利用するために使用します
以下に、octave.eval, octaveの使用例を示します。octave.evalでは結果が出力されますので、これを避けるためにjunk変数に代入しています。

# octaveで入力した文字列を評価させる場合 octave.eval('2+2') 
        
'ans = 4'
'ans = 4'

octaveの戻り値はOctaveElemntとして返されます。OctaveElementに対しては、 算術計算ができるためsageでのループ処理に利用できます。

また、sageで数値と扱うには、N関数を使って一度数値に変換して使用します。

# 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
<class 'sage.interfaces.octave.OctaveElement'> <type 'sage.rings.real_mpfr.RealNumber'>  4  12

octaveの参考サイト

octaveの使い方をコンパクトに説明したサイトが東海大学　稲葉研究室にあります。

# 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] $$

# 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
 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 $$ の結果を以下に示します。

# 要素毎への計算には.を付ける C = octave("A.*B"); C 
        
 4 4 12
 15 24 35
 49 64 81
 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] $$

# 左除算演算子\を使うと方程式や最小自乗法を一つの表現で処理できる # 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
 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]を求めると以下の様になります。

# 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
 1.5
 -0.666667

 

 3.83333
 5.33333
 6.83333

求まった結果をsageでグラフに表示

先の計算で求められたX=[a; b]の値を使って、解の直線y = a*x +bとデータの点をプロットしてみます。

a, bに代入する場合には、N関数を使ってoctave要素から数値に変換します。そして、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からマトリックスをoctaveに渡す場合

先の例では、octaveの結果（ベクトル）から要素を取り出し、sageで結果をプロットしましたが、 逆にsageのマトリックスをoctaveに渡すときには、octave.sage2octave_matrix_string を使用して、一端octaveのマトリックス表現文字列に変換します。

変換されたマトリックス表現文字列をoctaveに渡すことでsageからoctaveへのマトリックス 変換ができます。

# 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
[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関数を 呼び出すことで、簡単に行うことができます。

# octaveのプロット結果を表示する octave.eval('x = 0:0.1:10;') octave.eval('y = sin(x);') octave.eval('plot(x, y);') saveAndShowPlot('test.png') 
       

octaveは信号処理や画像処理が得意

octaveには、信号処理や画像処理を行うためのライブラリが充実しており、 それらをsageから使うことによってsageの処理の幅が広がります。

画像処理の例として、lena_std.jpgの画像をグレイスケールに変換してみます。

# octaveは信号処理や画像処理が得意 html('<img src="lena_std.jpg" width="256px">') # 画像を読み込む input = DATA+"lena_std.jpg" junk = octave.eval('rgb = imread("%s");'%input) 
       
# グレイスケールに変換 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">')