トップページ -> コンピュータ実習 -> Octave を用いた画像処理,信号処理 -> liboctave のベクトル,行列,統計計算に関する機能
[サイトマップへ], [サイト内検索へ],

liboctave のベクトル,行列,統計計算に関する機能

サイト構成 連絡先,業績など コンピュータ実習 データの扱い コンピュータ設定 教材(公開) サポートページ

liboctave とは,Octave の機能(全機能ではないが)が集まった C++ 言語のライブラリファイルとヘッダファイル. 自分で,C++ のプログラムを書き,liboctave の機能を呼び出すことが簡単にできます.

この Web ページでは, ベクトル,行列,統計計算に関する機能のうちごく一部を説明する.

【お断り】 以下,liboctave の機能を呼び出す C++ プログラムの見本と,コマンドラインでのビルドと実行の手順を説明する. 文法エラーは無いことは確認済みです.正しく計算できているかのチェックは済んでいません


 

前準備

  1. Octave のビルドとインストールに関する Web ページ の記述に従って,Octave のインストール

    ※ Octave 3.6 系列とそれ以前で liboctave の使用法に変化があったようで、このWebページでは、現時点での最新版である Octave 3.6 で説明を続けることにします.

  2. (オプション) ソースコードのコピー

    場合によってはソースコードがあった方が便利という場合があります(必ずしも必要というわけではありません)

    cd /tmp
    if [ ! -f octave-3.6.2.tar.gz ]; then 
        wget ftp://ftp.gnu.org/gnu/octave/octave-3.6.2.tar.gz
    fi
    
    sudo rm -rf octave-3.6.2
    tar -xvzof octave-3.6.2.tar.gz
    cd octave-3.6.2
    
    cd /tmp/octave-3.6.2/src/DLD-FUNCTIONS
    sudo cp *.cc /usr/local/include/octave-3.6.2/octave
    

 

コンパイルオブション

最初の見本は、行列のコンストラクタです.

#include <iostream>
#include <stdio.h>
#include <time.h>
#include <octave/config.h>
#include <octave/Matrix.h>

int main( int argc, char **argv )
{
  Matrix X(2, 3, 1.0);
  std::cout << X;

  return 0;
}

Ubuntu でのコンパイルオプションの例

上記のようなソースコードをa.cc のようなファイル名で保存しておく


 

行列(基本操作、行列と行列の積、LU分解、逆行列、行列式、SVD、QR分解、固有値と固有ベクトル、など)


 

Octave の DLD-FUNCTIONS を C++ から使いたい場合

  1. 使いたいライブラリファイルを探す

    ここの説明では,「/usr/local/lib/octave/3.6.2/rand.a」を使いたいとする

  2. 使いたい関数を探す

    ◆ 操作手順の例

    nm -C /usr/local/lib/octave/3.6.2/rand.a 
    

    ◆ 実行結果の例

    ここの説明では, 「Frand(octave_value_list const&, int)」を使いたいとする。

  3. ソースコードを書く

    要点は3つ

    ◆ ソースコードの例

    // liboctave を用いた 乱数
    #include <stdio.h>
    #include <time.h>
    #include <iostream>
    #include <octave/config.h>
    #include <octave/Matrix.h>
    #include <octave/dMatrix.h>
    #include <octave/defun-dld.h>
    
    #define DIM 3
    
    extern octave_value_list Frand(octave_value_list const&, int); 
    
    int main( int argc, char **argv )
    {
      Matrix X(DIM, DIM);
      std::cout << X;
    
      octave_value_list L;
      L.append(DIM);
      L.append(DIM);
      octave_value_list result = Frand(L, L.length());
      X = result(0).float_array_value (); 
      std::cout << X; 
    
      return 0;
    }
    

  4. ビルド

    <要点>

    リンクオプションとして /usr/local/lib/octave/3.6.2/rand.a を付ける

    ※ コンパイルオプションについては、この Web ページの上のほうで説明している.

    ◆ 実行結果の例

  5. 動作確認

    ./a.out
    

    ◆ 実行結果の例

コンボリューション (DLD_FUNCTIONS の Fconv2)

conv2.a の中の Fconv2() を使う。

◆ ソースコードの例

// liboctave を用いた 2次元畳み込み(コンボリューション)
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/defun-dld.h>

#define DIM 3
#define DIM2 2

extern octave_value_list Fconv2(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(DIM, DIM);
  Matrix B(DIM2, DIM2);
  Matrix Z(DIM, DIM);

  X(0,0) = 1; X(0,1) = 4; X(0,2) = 7; 
  X(1,0) = 2; X(1,1) = 5; X(1,2) = 8; 
  X(2,0) = 3; X(2,1) = 6; X(2,2) = 9; 
  B(0,0) = 1; B(0,1) = 1; 
  B(1,0) = 1; B(1,1) = 1; 
  std::cout << X;
  std::cout << B;
  fprintf( stderr, "Z = Fconv2(X, B, \"full\")\n" );


  octave_value_list L;
  L.append(X);
  L.append(B);
  L.append("full");
  octave_value_list result = Fconv2(L, L.length());
  Z = result(0).float_array_value (); 
  std::cout << Z; 

  return 0;
}

◆ 実行結果の例

一様乱数 (DLD_FUNCTIONS の Frand)

rand.a の中の Frand() を使う。

◆ ソースコードの例

// liboctave を用いた 乱数
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/defun-dld.h>

#define H 8
#define W 3


extern octave_value_list Frand(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(H, W);
  std::cout << X;

  octave_value_list L;
  L.append(H);
  L.append(W);
  octave_value_list result = Frand(L, L.length());
  X = result(0).float_array_value (); 
  std::cout << X; 

  return 0;
}

◆ 実行結果の例

正規乱数 (DLD_FUNCTIONS の Frandn)

分布が正規分布になるような乱数. rand.a の中の Frandn() を使う。

◆ ソースコードの例

// liboctave を用いた 乱数
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/defun-dld.h>

#define H 8
#define W 3


extern octave_value_list Frand(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(H, W);
  std::cout << X;

  octave_value_list L;
  L.append(H);
  L.append(W);
  octave_value_list result = Frandn(L, L.length());
  X = result(0).float_array_value (); 
  std::cout << X; 

  return 0;
}

◆ 実行結果の例

1次元のFFT (DLD_FUNCTIONS の fft)

fft.a の中の Ffft() を使う。

◆ ソースコードの例

#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/ov-flt-complex.h>
#include <octave/ov-complex.h>
#include <octave/fCMatrix.h>
#include <octave/CMatrix.h>
#include <octave/defun-dld.h>

#define LEN 256

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Ffft(octave_value_list const&, int); 
extern octave_value_list Fifft(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  ComplexMatrix X(1, LEN);
  ComplexMatrix Y(1, LEN);

  octave_value_list L;
  octave_value_list result;

  L.append(1);
  L.append(LEN);
  result = Frandn(L, L.length());
  X = result(0).complex_array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  L.append(LEN);
  result = Ffft(L, L.length());
  Y = result(0).complex_array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例

2次元のFFT (DLD_FUNCTIONS の fft)

fft.a の中の Ffft2() を使う。

◆ ソースコードの例

#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/ov-flt-complex.h>
#include <octave/ov-complex.h>
#include <octave/fCMatrix.h>
#include <octave/CMatrix.h>
#include <octave/defun-dld.h>

#define SZ 16

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Ffft2(octave_value_list const&, int); 
extern octave_value_list Fifft2(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  ComplexMatrix X(SZ, SZ);
  ComplexMatrix Y(SZ, SZ);

  octave_value_list L;
  octave_value_list result;

  L.append(SZ);
  L.append(SZ);
  result = Frandn(L, L.length());
  X = result(0).complex_array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  result = Ffft2(L, L.length());
  Y = result(0).complex_array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例

Convex Hull (DLD_FUNCTIONS の convhulln)

convhull.a の中の Fconvhull() を使う。

◆ ソースコードの例

#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/ov-flt-complex.h>
#include <octave/ov-complex.h>
#include <octave/fCMatrix.h>
#include <octave/CMatrix.h>
#include <octave/defun-dld.h>

#define DIM 3
#define NPOINTS 40

extern octave_value_list Frandn(octave_value_list const&, int); 
extern octave_value_list Fconvhulln(octave_value_list const&, int); 

int main( int argc, char **argv )
{
  Matrix X(NPOINTS, DIM);
  Matrix Y(NPOINTS, DIM);

  octave_value_list L;
  octave_value_list result;

  L.append(NPOINTS);
  L.append(DIM);
  result = Frandn(L, L.length());
  X = result(0).array_value ();  
  std::cout << "X\n"; 
  std::cout << X; 

  L.clear();
  L.append(X);
  result = Fconvhulln(L, L.length());
  Y = result(0).array_value ();  
  std::cout << "Y\n"; 
  std::cout << Y; 

  return 0;
}

◆ 実行結果の例


経過時間計測の例

ここでは、Octave の 関数 Fconv2() を用いて, 2次元の畳み込み(コンボリューション)を行う.

【要点】

// liboctave を用いた 2次元畳み込み(コンボリューション)
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <octave/config.h>
#include <octave/Matrix.h>
#include <octave/dMatrix.h>
#include <octave/defun-dld.h>

#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_statistics.h>


extern octave_value_list Fconv2(octave_value_list const&, int); 


// 正方行列のサイズを指定
#define DIM 3
#define DIM2 2

void matrix_rng_uniform( Matrix& M, const int N )
{
  // ここからは,GSL の機能を使って,乱数を配列 X, B に格納
  gsl_rng_env_setup();
  gsl_rng_type *T = (gsl_rng_type *)gsl_rng_default;
  /* 乱数発生器 */
  gsl_rng *r = gsl_rng_alloc(T);
  /* システムクロックを使って乱数の初期値を設定 */
  gsl_rng_set (r, time (NULL));

  int i, j;
  for(i = 0; i < N; i++) {
    for(j = 0; j < N; j++) {
      M(i, j) =  gsl_rng_uniform(r); 
    }
  }

  return;
}

int main( int argc, char **argv )
{

  Matrix X(DIM, DIM);
  Matrix B(DIM2, DIM2);
  Matrix Z(DIM, DIM);

  std::cout << X;
  std::cout << B;

  matrix_rng_uniform( X, DIM );
  matrix_rng_uniform( B, DIM2 );

  fprintf( stderr, "start, \n" );
  fprintf( stderr, "Z = Fconv2(X, B, \"full\"\n" );
  clock_t c = clock();

  octave_value_list L;
  L.append(X);
  L.append(B);
  L.append("full");
  octave_value_list result = Fconv2(L, L.length());
  Z = result(0).float_array_value (); 
  std::cout << Z; 

  fprintf( stderr, "done, elapsed time = %f [sec]\n", ( (double)clock() - (double)c ) / CLOCKS_PER_SEC );

  return 0;
}


ビルドの手順

◆ 実行結果の例


OpenMP を使うように DLD-FUNCTIONS/conv2.cc を書き換える.2 次元の畳み込みを高速

【要点】


Eclipse を使う場合

  1. Eclipse のインストールが済んでいること.

  2. 「Eclipse に C/C++ 開発ツール (CDT) をインストール」の Web ページの記述に従って,CDT プラグインのインストールが済んでいること.

  3. 「Octave の機能を呼び出す C++ プログラムの見本と実行手順」の Web ページの記述に従って,プロジェクトの作成と,プロジェクトのプロパティーの設定が済んでいること.

あらかじめ決めておく事項

Eclipse でのビルドと実行の手順

  1. (前準備)「Octave の機能を呼び出す C++ プログラムの見本と実行手順」の Web ページの記述に従って,プロジェクトの作成と,プロジェクトのプロパティーの設定が済んでいること.

  2. すべてビルド

    CTRL + B(コントロールキーを押しながら「B」) または,「プロジェクト」メニューで「すべてビルド」

    「コンソール」にコンパイルコマンド列が表示されるので,エラーが無いことを確認しておく.

  3. 実行

    デフォルトでは,実行ファイルは, C:\workspace\Hello\Debug\Hello.exe に出来る.(「Hello」はプロジェクト名)

    Eclipse のプロジェクト・エクスプローラで 「Debug」を展開.Hello.exe を右クリックし,「実行」から「ローカル C/C++ アプリケーション」を選ぶ. 数秒ほど待つと,コンソールに結果が出る.

    ※ Windows の Cygwin のコンソールを使うこともできる.。 C:\workspace\Hello\Debug\Hello.exe を実行. 数秒ほど待つと,コンソールに結果が出てくる.

    ※ Cygwin で実行がうまくいかない場合, C:\cygwin\bin\cygwin1.dll を C:\Windows\System32 にコピーすることを忘れている可能性がある.


(参考)liboctave のヘッダファイルのうち関係部分の抜粋(Octave バージョン3.2.0)

ソースコードの引用に先立つ著作権表示

/*

Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
              John W. Eaton
Copyright (C) 2008, 2009 Jaroslav Hajek

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
.

*/

Array.h

Array クラスは配列. MArray2 クラスの親クラスが Array2 クラスで,Array2 クラスの親クラスが Array

template <class T>
class
Array
{
  /* 省略 */
  void fill (const T& val);
  /* 省略 */
  octave_idx_type rows (void) const { return dim1 (); }
  octave_idx_type cols (void) const { return dim2 (); }
  /* 省略 */
  const dim_vector& dims (void) const { return dimensions; }
  /* 省略 */
  T& elem (octave_idx_type i, octave_idx_type j) { return elem (dim1()*j+i); }
  /* 省略 */
  T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); }
  /* 省略 */
  void print_info (std::ostream& os, const std::string& prefix) const;
};

MArray2.h

Marray2クラスは 2 次元配列.Matrix クラスの親クラス

template <class T>
class
MArray2 : public Array2<T>
{
  /* 省略 */
  MArray2 (octave_idx_type n, octave_idx_type m) : Array2<T> (n, m) { }
  MArray2 (octave_idx_type n, octave_idx_type m, const T& val) : Array2<T> (n, m, val) { }
  /* 省略 */
  MArray2<T>& insert (const Array2<T>& a, octave_idx_type r, octave_idx_type c)
  {
    Array2<T>::insert (a, r, c);
    return *this;
  }
};

dMatrix.h

Matrix クラスはdouble 型の 2 次元配列

class
OCTAVE_API
Matrix : public MArray2<double>
{
  /* 省略 */
  Matrix (octave_idx_type r, octave_idx_type c) : MArray2<double> (r, c) { }
  Matrix (octave_idx_type r, octave_idx_type c, double val) : MArray2<double> (r, c, val) { }
  /* 省略 */
  Matrix& insert (const Matrix& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const RowVector& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c);
  Matrix& insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c);
  /* 省略 */
  friend OCTAVE_API Matrix real (const ComplexMatrix& a);
  friend OCTAVE_API Matrix imag (const ComplexMatrix& a);
  /* 省略 */
  Matrix transpose (void) const { return MArray2<double>::transpose (); }
  /* 省略 */
  RowVector row (octave_idx_type i) const;
  ColumnVector column (octave_idx_type i) const;
  /* 省略 */
  Matrix inverse (void) const;
  /* 省略 */
  ComplexMatrix fourier2d (void) const;
  ComplexMatrix ifourier2d (void) const;
  /* 省略 */
  // Generic interface to solver with probing of type
  Matrix solve (const Matrix& b) const;
  ComplexMatrix solve (const ComplexMatrix& b) const;
  ColumnVector solve (const ColumnVector& b) const;

EIG.h

class
OCTAVE_API
EIG
{
  /* 省略 */
  EIG (const Matrix& a, bool calc_eigenvectors = true)
    { init (a, calc_eigenvectors); }
  /* 省略 */
  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
    { init (a, calc_eigenvectors); }
  /* 省略 */
  ComplexColumnVector eigenvalues (void) const { return lambda; }
  ComplexMatrix eigenvectors (void) const { return v; }
  /* 省略 */
};

dbleSVD.h

class
OCTAVE_API
SVD
{
  /* 省略 */
  enum type
    {
      std,
      economy,
      sigma_only
    };
  /* 省略 */
  SVD (const Matrix& a, type svd_type = SVD::std) { init (a, svd_type); }
  /* 省略 */
  DiagMatrix singular_values (void) const { return sigma; }
  Matrix left_singular_matrix (void) const;
  Matrix right_singular_matrix (void) const;
  /* 省略 */
};