サイトマップ サイト内検索 まとめ

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

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

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

【お断り】 以下，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. ソースコードを書く

    要点は３つ

    • 「extern octave_value_list Frand(octave_value_list const&, int)」を入れる
    • DLD-FUNCTIONS の関数の引数については、octave_value_list を使うのが便利．
    • 関数の返り値は octave_value_list 型の変数（例えば result）に入れ、result(0) のようにして取り出す

    ◆ ソースコードの例 

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

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

    • プロジェクト名： Hello
    • プロジェクト・タイプ： 実行可能

    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;
  /* 省略 */
};