金子邦彦研究室プログラミングOctave の活用liboctave のベクトル,行列,統計計算に関する機能

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」を使いたいとする

    [image]
  2. 使いたい関数を探す

    ◆ 操作手順の例

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

    ◆ 実行結果の例

    [image]

    ここの説明では, 「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 ページの上のほうで説明している.

    ◆ 実行結果の例

    [image]
  5. 動作確認

    ./a.out
    

    ◆ 実行結果の例

    [image]

    コンボリューション (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;
    }
    

    ◆ 実行結果の例

    [image]

    一様乱数 (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;
    }
    

    ◆ 実行結果の例

    [image]

    正規乱数 (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;
    }
    

    ◆ 実行結果の例

    [image]

    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;
    }
    

    ◆ 実行結果の例

    [image]

    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;
    }
    

    ◆ 実行結果の例

    [image]

    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;
    }
    

    ◆ 実行結果の例

    [image]

    経過時間計測の例

    ここでは、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;
    }
    
    
    

    ビルドの手順

    ◆ 実行結果の例

    [image]

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

    【要点】

    [image]

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