読者です 読者をやめる 読者になる 読者になる

C++ meets MATLAB

この記事は C++ Advent Calendar 23日目の記事です.
前日の記事はhgodaiさんのenumを便利にするユーティリティでした.

諸注意

記事内のC++に関する言及は間違っているかもしれません.
そしてMATLAB&Simulinkについても完全に初心者なので,間違いがあった場合は@まで指摘お願いします.
また,Fallout 4はSteamで好評発売中です.

f:id:Yonosamochi:20151218130302p:plain
figure.1 うっかりC++ Advent calendarに参加してMATLABの話題を書こうとした私の図

MATLAB

MATLABは、MATrix LABoratoryを略したものであり、行列計算、ベクトル演算、グラフ化や3次元表示などの豊富なライブラリを持った、インタプリタ形式の高性能なテクニカルコンピューティング言語、環境としての機能を持つ。標準で数多くのライブラリを有しているが、それ以上のデータ解析や統計、アプリケーション展開などが必要な場合にはToolboxと呼ばれる拡張パッケージをインストールすることで、MATLABの機能拡張を図ることができる。 - Wikipediaより

要するに統合開発環境であり,プログラミング言語です.
機械学習や画像処理,信号処理などのライブラリ ― MATLABではこれらをToolboxと呼称しています― もかなり充実しています.
そして名前,"Matrix Laboratory"の通りに,特に行列演算を得意とします.

MATLABの基本的な構文と処理

これはC++ Advent Calendarであり,MATLAB Advent Calendarではありません.ですが,なにも言わずに突然MATLABコードを記事に載せても混乱を招くので,MATLABのごく基本的な構文は紹介しておきます.
MATLABは動的型付け言語です.実行時に代入される変数のみが型と値を持ちます.

r = 9;

この変数aは1*1――要するにスカラー型――で,doubleとして宣言されています.
9.0でもないのにdoubleというのは少し奇妙かもしれませんが,MATLABでは倍浮動小数点型が標準です.
でもご安心ください!ちゃんと整数型もMATLABではそろっています.

z = int64(9);

はい.何言ってるんですかね.傍目にはキャストしているようにしか見えません.ここらへんにツッコむとガンガン文字数が増えていくので後は淡々と紹介していきます.

スカラー型の四則演算

r_1 = 114;
r_2 = 514;

add = r_1 + r_2; % 加算 =628
subtract = r_1 - r_2; % 減算 =-400
multiple = r_1 * r_2; % 乗算 =58596
division = r_1 / r_2; % 除算 =0.2218
surplus_1 = rem(r1, r2); % 剰余 =114
surplus_2 = mod(r1, r2); % モジュロ演算による剰余 =114

ベクトルと行列

colVec = [5; 1; 4];
rowVec = [1, 1, 4]; % 実はセミコロンが不要なので,[1 1 4]でもOK.
mat1 = colVec * rowVec;
mat2 = mat_1 * 10;
mat3 = mat_1 * eye(3);

この場合,各変数は

\begin{equation} colVec= \begin{bmatrix} 5 \\ 1 \\ 4 \\ \end{bmatrix} \end{equation}

\begin{equation} rowVec= \begin{bmatrix} 5 &1 &4 \end{bmatrix} \end{equation}

\begin{equation} mat1 \begin{bmatrix} 5 &5 &20 \\ 1 &1 &4 \\ 4 &4 &16 \end{bmatrix} \end{equation}

\begin{equation} mat2= \begin{bmatrix} 50 &50 &200 \\ 10 &10 &40 \\ 40 &40 &160 \end{bmatrix} \end{equation}

\begin{equation} mat3= \begin{bmatrix} 5 &5 &20 \\ 1 &1 &4 \\ 4 &4 &16 \end{bmatrix} \end{equation} となります.
mat3のeye()関数は,単位行列を生成する関数です.この場合は3*3の単位行列,\begin{equation} I= \begin{bmatrix} 1 &0 &0 \\ 0 &1 &0 \\ 0 &0 &1 \end{bmatrix} \end{equation} を生成します.
もしあなたが行列を全然知らない方で,どうして変数mat2, mat3のような行列になるか分からない方がいても,この記事では無視してもらっても構いません.

関数

当然ですが,関数も自分で定義可能です.たとえば引数XをZ倍する関数multの場合を考えます.

% filename: mult.m

function R = mult(X, Z)
R = X*Z;

だいたい予想はつくと思いますが,この関数は2つの引数を取り,ひとつの値を返却します.
特に説明はいらないとは思いますのでこの関数では省きます.
もうひとつ関数を考えましょう.今度は値を2つ以上返却する関数です.

% filename: circleinfo.m
% R は円の半径
function[C, S] = circleinfo(R)
C = R*2*pi; % 半径をRとする円の円周
S = R*R*pi; % 半径をRとする円の面積

見ればわかりますね.コメントの通り半径をRとする円の円周と面積を求めています.piはMATLABで用意されている定数です.
この関数circleinfo(R)を使用するには次のように関数を呼び出します

[Circumference, Area] = circleinfo(1);

この場合はCircumference = 6.82, Area = 3.1416になります.

さて,これで基本的な構文については(色々説明不足ですが)終わりです. 次からが今回の本題,Matlab coderの話です.

MATLAB Coder

MATLABにはMATLAB CoderというMATLABコードからC/C++コードを生成するtoolboxが存在します.
これを使用すればMATLABのコードをC++から使えるということですね.べんり.
また,C++のみではなく,ダイナミックライブラリ,スタティックライブラリの生成もサポートしています.そして,マルチコア対応の関数生成もOpenMPによってサポートされます.まぁ一通りの機能は揃っている,ということですね.
細かいことはここらへんに載っています.
今回はどのようなC++コードが吐き出されるかを中心に据えます.

生成

MATLAB CoderでC/C++コードを生成する場合,エントリポイントとなる関数が必要になります.
例として,常に定数1を返却する関数oneをC++に落とし込んでみましょう.

% filename: one.m
function R = one()
R = 1;

この関数oneは次のC++コードを生成します…

one.h

/*
 * Academic License - for use in teaching, academic research, and meeting
 * course requirements at degree granting institutions only.  Not for
 * government, commercial, or other organizational use.
 *
 * one.h
 *
 * Code generation for function 'one'
 *
 */

#ifndef __ONE_H__
#define __ONE_H__

/* Include files */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "tmwtypes.h"
#include "mex.h"
#include "emlrt.h"
#include "blas.h"
#include "rtwtypes.h"
#include "one_types.h"

/* Function Declarations */
extern real_T one();

#ifdef __WATCOMC__

#pragma aux one value [8087];

#endif
#endif

/* End of code generation (one.h) */

one.cpp

/*
 * Academic License - for use in teaching, academic research, and meeting
 * course requirements at degree granting institutions only.  Not for
 * government, commercial, or other organizational use.
 *
 * one.cpp
 *
 * Code generation for function 'one'
 *
 */

/* Include files */
#include "rt_nonfinite.h"
#include "one.h"

/* Function Definitions */
real_T one()
{
  return 1.0;
}

/* End of code generation (one.cpp) */

ぶっちゃけた話が一番上のライセンス条項が怖い
他にもっと生成されるファイルが多いのですが,マジで多いので省略します.

さて,この関数one()を実行実際に使ってみます.mainファイルは次のようになります…

main.cpp

#include <iostream>
#include "one/one.h"
#include "one/one_initialize.h"
#include "one/one_terminate.h"

int main(int argc, char const* argv[])
{
    one_initialize();
    std::cout << one() << std::endl;
    one_terminate();
    return 0;
}

どのような出力になるかはまぁ自明ですね.1が表示されます.
問題はone_initialize();one_terminate();です.まぁ関数名からわかるとは思いますが,これらはMATLAB Coderによって生成される関数を使用する際に呼び出しておくものです.
one()関数の呼び出す回数にかかわらず,initalize()もterminate()も1回呼び出せば十分です.

次はもう少し複雑な関数――先ほど例として挙げたcircleinfo関数――をC++に落とし込んでみましょう.
今度はmain関数も自動生成してみます.MATLAB Coderは一例としてその関数を使用するmainファイルを生成することもできます.

circleinfo.cpp

//
// Arguments    : double R
//                double *C
//                double *S
// Return Type  : void
//
void circleinfo(double R, double *C, double *S)
{
  //  filename: circleinfo.m
  *C = R * 2.0 * 3.1415926535897931;
  *S = R * R * 3.1415926535897931;
}

main.cpp

// Include Files
#include "rt_nonfinite.h"
#include "circleinfo.h"
#include "main.h"
#include "circleinfo_terminate.h"
#include "circleinfo_initialize.h"

// Function Declarations
static double argInit_real_T();
static void main_circleinfo();

// Function Definitions

//
// Arguments    : void
// Return Type  : double
//
static double argInit_real_T()
{
  return 0.0;
}

//
// Arguments    : void
// Return Type  : void
//
static void main_circleinfo()
{
  double S;
  double C;

  // Initialize function 'circleinfo' input arguments.
  // Call the entry-point 'circleinfo'.
  circleinfo(argInit_real_T(), &C, &S);
}

//
// Arguments    : int argc
//                const char * const argv[]
// Return Type  : int
//
int main(int, const char * const [])
{
  // Initialize the application.
  // You do not need to do this more than one time.
  circleinfo_initialize();

  // Invoke the entry-point functions.
  // You can call entry-point functions multiple times.
  main_circleinfo();

  // Terminate the application.
  // You do not need to do this more than one time.
  circleinfo_terminate();
  return 0;
}

前置文は長いのですべてカットしています.
元のMATLAB関数が返り値を2つ持っていますが,C++ではこれは厳しいので,引数にポインタを渡すようになっていますね.
まぁ,C++でもstd::tupleを使えば多値返却は可能ですが,これはC++'11以降のみのテクニックですから,移植性を売りにしているMATLAB Coderでは論外というところでしょうか.[要出典]

最後に,引数に行列をひとつ渡し,その転置行列との積を計算するプログラムです.MATLABでは次のように表現されます...

fuckthekonami.m

function R = fuckthekonami(A)
R = A*A';

関数名に何か別の意味どころかまるで何かに向けて敵意を感じさせるような名前がついていますが,気にしないでください.仕様です.
さて,今更ですが,MATLABは動的型付け言語です.変数Aが行列型だろうと,スカラー型の値を代入すれば,それはスカラー型になります.しかし,これはC++ではまず許されません.そのため,MATLAB Coderでは引数がどんな型なのか,事前に指定する必要があります.
f:id:Yonosamochi:20151221145455p:plain
GUIMATLAB Coderではこのように型を指定します.singleは所謂単精度浮動小数点型ですね.
さらに,引数がスカラー型なのか,(行/列)ベクトルなのか,行列なのかも指定します.さらについでに言うとベクトルや行列だった場合,その数も指定します...
template?そんなものは甘えだ.

そしてこれがfuckthekonami関数をC++に落とし込むとできる関数です.

fuckthekonami.cpp

//
// Arguments    : const emxArray_real_T *A
//                emxArray_real_T *R
// Return Type  : void
//
void fuckthekonami(const emxArray_real_T *A, emxArray_real_T *R)
{
  emxArray_real_T *b;
  int i0;
  int br;
  int ar;
  int i1;
  int c;
  int cr;
  int k;
  unsigned int unnamed_idx_0;
  unsigned int unnamed_idx_1;
  int m;
  int ic;
  int ib;
  int ia;
  emxInit_real_T(&b, 2);
  i0 = b->size[0] * b->size[1];
  b->size[0] = A->size[1];
  b->size[1] = A->size[0];
  emxEnsureCapacity((emxArray__common *)b, i0, (int)sizeof(double));
  br = A->size[0];
  for (i0 = 0; i0 < br; i0++) {
    ar = A->size[1];
    for (i1 = 0; i1 < ar; i1++) {
      b->data[i1 + b->size[0] * i0] = A->data[i0 + A->size[0] * i1];
    }
  }

  if ((A->size[1] == 1) || (b->size[0] == 1)) {
    i0 = R->size[0] * R->size[1];
    R->size[0] = A->size[0];
    R->size[1] = b->size[1];
    emxEnsureCapacity((emxArray__common *)R, i0, (int)sizeof(double));
    br = A->size[0];
    for (i0 = 0; i0 < br; i0++) {
      ar = b->size[1];
      for (i1 = 0; i1 < ar; i1++) {
        R->data[i0 + R->size[0] * i1] = 0.0;
        c = A->size[1];
        for (cr = 0; cr < c; cr++) {
          R->data[i0 + R->size[0] * i1] += A->data[i0 + A->size[0] * cr] *
            b->data[cr + b->size[0] * i1];
        }
      }
    }
  } else {
    k = A->size[1];
    unnamed_idx_0 = (unsigned int)A->size[0];
    unnamed_idx_1 = (unsigned int)b->size[1];
    i0 = R->size[0] * R->size[1];
    R->size[0] = (int)unnamed_idx_0;
    R->size[1] = (int)unnamed_idx_1;
    emxEnsureCapacity((emxArray__common *)R, i0, (int)sizeof(double));
    m = A->size[0];
    i0 = R->size[0] * R->size[1];
    emxEnsureCapacity((emxArray__common *)R, i0, (int)sizeof(double));
    br = R->size[1];
    for (i0 = 0; i0 < br; i0++) {
      ar = R->size[0];
      for (i1 = 0; i1 < ar; i1++) {
        R->data[i1 + R->size[0] * i0] = 0.0;
      }
    }

    if ((A->size[0] == 0) || (b->size[1] == 0)) {
    } else {
      c = A->size[0] * (b->size[1] - 1);
      cr = 0;
      while ((m > 0) && (cr <= c)) {
        i0 = cr + m;
        for (ic = cr; ic + 1 <= i0; ic++) {
          R->data[ic] = 0.0;
        }

        cr += m;
      }

      br = 0;
      cr = 0;
      while ((m > 0) && (cr <= c)) {
        ar = -1;
        i0 = br + k;
        for (ib = br; ib + 1 <= i0; ib++) {
          if (b->data[ib] != 0.0) {
            ia = ar;
            i1 = cr + m;
            for (ic = cr; ic + 1 <= i1; ic++) {
              ia++;
              R->data[ic] += b->data[ib] * A->data[ia];
            }
          }

          ar += m;
        }

        br += k;
        cr += m;
      }
    }
  }

  emxFree_real_T(&b);
}

f:id:Yonosamochi:20151221150344p:plain
figure 2: C++が嫌い

読みたくないですね.眼が潰れそうです.
C++がわからない私にはなにをやってるのかサッパリですが,とにかくこいつは引数の行列とその転置行列の行列積を計算して返す処理をしています.少なくとも,MATLABだった時はそういう処理をしていました.
ご多分に漏れず,今回もmainまで自動生成しています.

main.cpp

// Include Files
#include "rt_nonfinite.h"
#include "fuckthekonami.h"
#include "main.h"
#include "fuckthekonami_terminate.h"
#include "fuckthekonami_emxAPI.h"
#include "fuckthekonami_initialize.h"

// Function Declarations
static double argInit_real_T();
static emxArray_real_T *c_argInit_UnboundedxUnbounded_r();
static void main_fuckthekonami();

// Function Definitions

//
// Arguments    : void
// Return Type  : double
//
static double argInit_real_T()
{
  return 0.0;
}

//
// Arguments    : void
// Return Type  : emxArray_real_T *
//
static emxArray_real_T *c_argInit_UnboundedxUnbounded_r()
{
  emxArray_real_T *result;
  static int iv0[2] = { 2, 2 };

  int idx0;
  int idx1;

  // Set the size of the array.
  // Change this size to the value that the application requires.
  result = emxCreateND_real_T(2, *(int (*)[2])&iv0[0]);

  // Loop over the array to initialize each element.
  for (idx0 = 0; idx0 < result->size[0U]; idx0++) {
    for (idx1 = 0; idx1 < result->size[1U]; idx1++) {
      // Set the value of the array element.
      // Change this value to the value that the application requires.
      result->data[idx0 + result->size[0] * idx1] = argInit_real_T();
    }
  }

  return result;
}

//
// Arguments    : void
// Return Type  : void
//
static void main_fuckthekonami()
{
  emxArray_real_T *R;
  emxArray_real_T *A;
  emxInitArray_real_T(&R, 2);

  // Initialize function 'fuckthekonami' input arguments.
  // Initialize function input argument 'A'.
  A = c_argInit_UnboundedxUnbounded_r();

  // Call the entry-point 'fuckthekonami'.
  fuckthekonami(A, R);
  emxDestroyArray_real_T(R);
  emxDestroyArray_real_T(A);
}

//
// Arguments    : int argc
//                const char * const argv[]
// Return Type  : int
//
int main(int, const char * const [])
{
  // Initialize the application.
  // You do not need to do this more than one time.
  fuckthekonami_initialize();

  // Invoke the entry-point functions.
  // You can call entry-point functions multiple times.
  main_fuckthekonami();

  // Terminate the application.
  // You do not need to do this more than one time.
  fuckthekonami_terminate();
  return 0;
}

f:id:Yonosamochi:20151221150344p:plain
Figure.3 C++がわからない時に使う画像

私には何をやっているのかはさっぱりわかりませんが,とにかく関数を実行しているんでしょう.
最後に,入力された値に複素数iの積を返す関数です.

returnComp.m

% filename: returnComp.m
function R = returnComp(V)
R = V*1i;

一応説明しておくとiは二乗すると-1になる数です\begin{equation} \sqrt{-1} \end{equation} といったところですね. C++ではstd::complexがありますから,きっとMATLAB Coderはstd::complexを使ったコードを生成してくれるでしょう.

returnComp.cpp

// Function Definitions

//
// Arguments    : double V
// Return Type  : creal_T
//
creal_T returnComp(double V)
{
  creal_T R;

  //  filename: returnComp.m
  R.re = V * 0.0;
  R.im = V;
  return R;
}

一体なぜstd::complexは生まれてきたのでしょうか.このcreal_TはMATLAB Coderで生成される独自の型です.
そもそも,行列を返す関数の時もemxArray_real_TというC++にはない型を生成していました.しかし,そもそもC++には行列型はありませんから,当然ですが.

なぜこのようなコードが…

流石に自動生成となると,心のこもっていないコードが生成されますね.人の手で書かれてこそのC++だと改めて実感しました.(冗談です)

注:ここからは少々曖昧な文言になっています.

MATLAB CoderはC++の機能を忌避しているような節がそこかしこに見受けられます.
templateは生成されるコードに含まれませんし,STLも一切使用されていません.
恐らくは移植性の問題です--(ここからは筆者の推測です)
C++は,古株ながら,現在も成長し続けている言語です.C++'11や'14ですね.clangやgccは'11や'14の機能を早めに取り込みました[要出典].しかし他のコンパイラは違います.msvcとか,intel c++ compiler等がそれですね.
一例として,C++'14の機能の一つ,constexpr関数の制限緩和をclangとintel C++ Compilerで比較しますと,clangでは3.4で実装されているのに対し,intel C++ CompilerではVersion 15でもまだ実装されていません. extended clangは単なる一例ですが,もし商用コンパイラまでサポートだけでなく,古いコンパイラのまま開発しているデベロッパーにまで合わせるとなると,templateが使用されていないのも納得です.

また,STLについては恐らくこれだけではなく,SIMULINKの問題があるのでしょう.MATLABにはSIMULINKという環境があり,これらはC/C++の関数をSIMULINKに読み込ませる上で,STLが吐く例外がSIMULINKの動作に影響を与えるそうです.

そして,生成されたC++コードはとてもじゃないですがreadableとは言いがたいですね.(この辺りは各人のリーディングスキルも関わってきますが,少なくとも私には人が読むコードには見えません).
しかし,人の手を加えないかぎりは生成されたコードにバグが入る可能性は(MATLABのコードにバグがある場合を除いて)ほとんどないでしょうから,気にする必要はないでしょう.つまり本当は実装なんて見る必要はなく,ブラックボックスなライブラリとしてなにも考えずに利用できるメリットがあります.
つまりどんな魔術で動いていたとしても,我関せずを貫き通せるということです.Boostみたいにね.

MATLAB Coderに付き纏う制限

MATLABコードなら何でもC++に落とし込めるわけではありません.C++コードに変換できないMATLABコードも勿論あります.
あるtoolboxのある関数がC++への生成を許していないというものです.例えば,Image Processing Toolboxのimread関数はC++コード生成をサポートしません.

X = imread('./大川隆法.jpg'); % コード生成不可

そして,疎行列(スパース行列),ラムダ式(無名関数),セル配列,入れ子関数,再帰,例外(try~catch ステートメント)はサポートされません.
これらの機能が使用されたコードは生成前にエラーとなります.

さらに,MATLABは動的型付け言語です(本日3回目).MATLABコードで動的な型付けの特徴を悪用したコードは,C++コード生成では控えるべきです(もっとも,処理速度の低下に繋がるので,そもそもMATLABコードでもやるなという指南もありますし,そもそもコードを検査する時点で弾かれるでしょうが)
例えば同じ変数に大きさの違う行列を代入する行為ですね.

ちなみに

MATLABコード内に定数式が存在した場合は,可能な限り生成時に処理を済ませてしまいます.実行時には処理されません.

wkaru.m

function R = wkaru()
R = [1, 2, 3; 4, 5, 6];
R = R * R';

このコードはC++コード生成時に計算されます.

wkaru.cpp

void wkaru(double R[4])
{
  int i0;
  static const signed char iv0[4] = { 14, 32, 32, 77 };

  for (i0 = 0; i0 < 4; i0++) {
    R[i0] = iv0[i0];
  }
}

お前は行列を返すんじゃないのか!?

f:id:Yonosamochi:20151223101355j:plain
figure 4: 行列を返す関数なのに配列を返されて激高するJOJO

おわりに

C++ Advent Calendar 23日目の記事,如何でしたでしょうか.
今回は紹介しませんでしたが,MATLABでは他にSimulink CoderというSimulinkコードをC/C++に落とし込めたり,Embedded Coderで生成されるC/C++コードを細かく制御したり,Raspberry PiArduinoの様な組み込みハード向けにC/C++や,HDLコードに変換したりするtoolboxがそろっています.

もし,この記事がポエム Advent Calendarの記事になるようでしたら,自殺 Advent Calendarの記事になるかと思います.
冗談です.
ご指摘があるようでしたら,@に殴りつけてください.C++初心者なので,もっと勉強する必要がありますから,どうかご教授をお願いします.

Figure 2, 3で使用したC++がわからないかわいい女の子の画像は@氏のご厚意によるものです.ありがとうございました.

明日24日めの記事はwordijpさんの「boost.RangeのOven拡張を自作してみる」です.

広告を非表示にする