2018年8月17日

std::vector<T> で MATLAB や numpy のように配列計算する

Abstract
  C++ において,クラスのコンストラクタ中でメモリを確保し,ディストラクタ中でメモリを解放することで,メモリ寿命をクラス寿命としてスコープにより制御すること (RAII※) は一般的である.C++ において,単なる配列を RAII に従い実装する場合,新たなクラスを定義することはせず,std::array<T> や std::vector<T> などの標準実装のコンテナ型を用いることが一般的である.これらのコンテナ型は,単なる配列としてだけでなく,数値計算やデータ解析における配列としても多用される.このとき「要素の追加」や「配列全体への演算」といった操作が頻繁に発生する.しかしながら,現状では std::vector<T> に演算子が定義されておらず,MATLAB における a=[1, 2, 3]; a*=4 や numpy における a=np.asarray([1,2,3]); a*=4 といった演算は,std::vector<double> a={1,2,3}; for(int i=0; i<a.size(); i++){ a[i]*=4; } のように,for 文を用いた冗長な記述となる.
  そこで,本投稿では,計算や解析で特に多用される std::vector<T> に対して演算子を定義することで,記述の簡略化を図る.より具体的には,四則演算子と剰余演算子,べき乗演算子を定義した.また,左シフト演算子 (<<) と 左シフト代入演算子 (<<=) をそれぞれ要素の結合と追加 (.push_back()) として割り当てた.この結果,先ほどの演算は,std::vector<double> a={1,2,3}; a*=4; のように簡略化された.

※ RAII は Resource Acquisition Is Initialization の省略形.「資源(リソース)の確保と解放を,クラス型の変数の初期化と破棄処理に結び付けるというプログラミングのテクニック [追記1]」のこと.
Introduction
  C++ において std::vector<T> を用いてメモリを確保し,その寿命をスコープで制御することは一般的である.そのため,数値計算やデータ解析における配列として多用され,要素の追加や配列全体への演算といった操作が頻繁に発生する.また, std::vector<T> で扱うデータ型 T は,double 型や float 型などの浮動小数点数型だけでなく,int 型や MPFR の任意精度の浮動小数点数型を定義した class 等が用いられることもある.例えば,double 型の配列を用いる場合,各配列値を 4 倍するためには,
  Code
#include <iostream>
#include <vector>

template <typename T>
inline void print(const std::vector<T>& rhs){
    std::cout<<" = [ ";
    for(int i=0; i<rhs.size(); i++){ std::cout<<rhs[i]<<' '; }
    std::cout<<"]"<<std::endl;
}
#define printn(var) {printf("%s", #var);print(var);}
#define printn_all(var) {printf("%s(%d): ", __func__, __LINE__);printf("%s", #var);print(var);}

int main(){
    std::vector<double> a={1, 2, 3};
    for(int i=0; i<a.size(); i++){ a[i]*=4; }
    printn(a);

    return 0;
}
  Execution result
a = [ 4 8 12 ]
  Running with Ideone.com
のように,for 文を用いて記述する必要がある.
Previous research
  MATLAB や Python の numpy では,演算子が定義されているため,for 文を用いずに表記することが可能である.下記に,Python における numpy の例を示す.
  Code
import numpy as np

def main():
    a = np.asarray([1,2,3])
    a*=4
    print(a)

main()
  Execution result
[ 4  8 12]
  Running with Ideone.com
Purpose
  本投稿では,演算子 +, +=, -, -=, *, *=, /, /=, %, %=, ^, ^= をオーバーロードし,最初に示したサンプルコードを,
int main(){
    std::vector<double> a={1, 2, 3};
    a*=4;
    printn(a);

    return 0;
}
と記述可能とすることを目的とする.加えて,頻出する要素追加演算 .push_back() を <<= 演算子として割り当て,同様に << 演算子を std::vector の連結演算子として割り当てることを目的とする.また,built-in 型だけでなく,新たなに class 等で定義された演算子に対しても,std::vector<T> に定義する演算子が機能するように,テンプレートを用いた汎用的な実装を目標とする.
Method
  演算子の定義は,演算子をオーバーロード (多重定義) することで達成する.オーバーロード (多重定義) とは,言葉の通り,既に演算が定義されている関数や演算子について,異なる引数をとる定義を追加する,即ち,一つの対象 (関数や演算子) に対して,渡される引数に応じた複数の解釈を多重に定義することである.例えば,built-in 型であれば,はじめから,
double operator+(const double& lhs, const double& rhs){ return lhs+rhs; }
が定義されているのと等価である.この既にある型について定義されている演算子について,引数に std::vector<T> が含まれる場合の動作も定義したい.
  あるテンプレート class に対して,+ 演算子を定義する場合を考える.このとき,class 内部で演算子を宣言するには,
  Code
#include <iostream>
#include <vector>

template <typename T> class example; // in order to extern

template <typename T> inline class example<T> example_add(const class example<T>& lhs, const class example<T>& rhs){ return example<T>(lhs.val()+rhs.val()); }
template <typename T> inline class example<T> example_add(const class example<T>& lhs, const               T & rhs){ return example<T>(lhs.val()+rhs); }

template <typename T>
class example{
private:
    T value;
public:
    example(){ value=(T)0; }
    example(const T& rhs){ value=rhs; }
    ~example(){}
    T val() const { return value; }
    
    class example<T> operator+(const class example<T>& rhs){ return example_add(*this, rhs); }
    class example<T> operator+(const               T & rhs){ return example_add(*this, rhs); }
};

template <typename T>
void print(const class example<T>& rhs){ std::cout<<" = "<<rhs.val()<<std::endl; }

#define printn(var) {printf("%s", #var);print(var);}
#define printn_all(var) {printf("%s(%d): ", __func__, __LINE__);printf("%s", #var);print(var);}

int main(){
    class example<double> v1(1.0), v2(2.0), v3(3.0), buf;
    printn(v1);
    printn(v2);
    printn(v3);
    printn(v1+v2+v3);
    
    buf=v1+v2+v3+4.0;
    printn(buf);
    
    return 0;
}
  Execution result
v1 = 1
v2 = 2
v3 = 3
v1+v2+v3 = 6
buf = 10
  Running with Ideone.com
となる.なお,このように内部で定義する場合,2 項演算子の右辺値しか受け取ることができない.したがって,上記の演算で左辺値が必要となる「buf=4.0+v1+v2+v3;」のような演算を行おうとすると,コンパイルエラーとなる.そのため,今回のように,既にあるクラスに変更を加えないまま演算子を定義する場合や,2 項演算子の左辺値が必要な場合は,
#include <iostream>
#include <vector>

template <typename T>
class example{
private:
    T value;
public:
    example(){ value=(T)0; }
    example(const T& rhs){ value=rhs; }
    ~example(){}
    
    T val() const { return value; }
};

template <typename T> inline class example<T> operator+(const class example<T>& lhs, const class example<T>& rhs){ return example<T>(lhs.val()+rhs.val()); }
template <typename T> inline class example<T> operator+(const class example<T>& lhs, const               T & rhs){ return example<T>(lhs.val()+rhs      ); }
template <typename T> inline class example<T> operator+(const               T & lhs, const class example<T>& rhs){ return example<T>(lhs      +rhs.val()); }

template <typename T>
void print(const class example<T>& rhs){ std::cout<<" = "<<rhs.val()<<std::endl; }

#define printn(var) {printf("%s", #var);print(var);}
#define printn_all(var) {printf("%s(%d): ", __func__, __LINE__);printf("%s", #var);print(var);}

int main(){
    class example<double> v1(1.0), v2(2.0), v3(3.0), buf;
    printn(v1);
    printn(v2);
    printn(v3);
    printn(v1+v2+v3);
    
    buf=0.1+v1+v2+v3+4.0;
    printn(buf);
    
    return 0;
}
  Execution result
v1 = 1
v2 = 2
v3 = 3
v1+v2+v3 = 6
buf = 10.1
  Running with Ideone.com
のように,外部で定義する.
  また,実際の実装においては,それぞれの演算子について,ほとんど同じ実装を行う必要がある.そのまま実装を行うと,ソースコードの見通しが非常に悪くかつ冗長となるため,各演算子における実装を関数形式マクロにより展開することで,簡潔な実装とした.
Implementation
  実際の実装を下記に示す.下記のコードでは,演算子の定義の主要な箇所を,コメントを含め 193 行の実装で達成している.

  ./sstd/src/stdVector_expansion/stdVector_expansion.hpp
#pragma once
#include "../math.hpp"

//-----------------------------------------------------------------------------------------------------------------------------------------------

#define SSTD_DEF_stdVecEx_defInNamespace(Func)                            \
    template <typename T>                   std::vector<T>  Func(const std::vector<T>& lhs, const std::vector<T>& rhs); \
    template <typename T, typename rhsType> std::vector<T>  Func(const std::vector<T>& lhs, const        rhsType& rhs); \
    template <typename T, typename lhsType> std::vector<T>  Func(const        lhsType& lhs, const std::vector<T>& rhs);
#define SSTD_DEF_stdVecEx_defInNamespace_eq(Func)                        \
    template <typename T>                   std::vector<T>& Func(      std::vector<T>& lhs, const std::vector<T>& rhs); \
    template <typename T, typename rhsType> std::vector<T>& Func(      std::vector<T>& lhs, const        rhsType& rhs);

namespace sstd_stdVecEx{
    // operators for mathematics
    SSTD_DEF_stdVecEx_defInNamespace   (add   ); // +
    SSTD_DEF_stdVecEx_defInNamespace_eq(add_eq); // +=
    SSTD_DEF_stdVecEx_defInNamespace   (sub   ); // -
    SSTD_DEF_stdVecEx_defInNamespace_eq(sub_eq); // -=
    SSTD_DEF_stdVecEx_defInNamespace   (mul   ); // *
    SSTD_DEF_stdVecEx_defInNamespace_eq(mul_eq); // *=
    SSTD_DEF_stdVecEx_defInNamespace   (div   ); // /
    SSTD_DEF_stdVecEx_defInNamespace_eq(div_eq); // /=
    SSTD_DEF_stdVecEx_defInNamespace   (mod   ); // %
    SSTD_DEF_stdVecEx_defInNamespace_eq(mod_eq); // %=
    SSTD_DEF_stdVecEx_defInNamespace   (pow   ); // ^
    SSTD_DEF_stdVecEx_defInNamespace_eq(pow_eq); // ^=
    
    // operators for std::vector
    SSTD_DEF_stdVecEx_defInNamespace   (push_back   ); // <<
    SSTD_DEF_stdVecEx_defInNamespace_eq(push_back_eq); // <<=
}

#undef SSTD_DEF_stdVecEx_defInNamespace    // Deletion of used definition, in order not to pollute the namespace
#undef SSTD_DEF_stdVecEx_defInNamespace_eq // Deletion of used definition, in order not to pollute the namespace

//-----------------------------------------------------------------------------------------------------------------------------------------------

// operators for mathematics
#define SSTD_DEF_stdVecEx_o(Func, Ope)                                \
    template <typename T>                                                \
    inline std::vector<T> Func(const std::vector<T>& lhs, const std::vector<T>& rhs){ \
        std::vector<T> ret(lhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=lhs[p] Ope rhs[p]; }    \
        return ret;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T> Func(const std::vector<T>& lhs, const rhsType& rhs){ \
        std::vector<T> ret(lhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=lhs[p] Ope rhs; }        \
        return ret;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T> Func(const rhsType& lhs, const std::vector<T>& rhs){ \
        std::vector<T> ret(rhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=lhs Ope rhs[p]; }        \
        return ret;                                                        \
    }
#define SSTD_DEF_stdVecEx_o_eq(Func, Ope)                                \
    template <typename T>                                                \
    inline std::vector<T>& Func(std::vector<T>& lhs, const std::vector<T>& rhs){ \
        for(uint p=0; p<lhs.size(); p++){ lhs[p] Ope rhs[p]; }            \
        return lhs;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T>& Func(std::vector<T>& lhs, const rhsType& rhs){ \
        for(uint p=0; p<lhs.size(); p++){ lhs[p] Ope rhs; }                \
        return lhs;                                                        \
    }
#define SSTD_DEF_stdVecEx_f(Func, Func2)                                \
    template <typename T>                                                \
    inline std::vector<T> Func(const std::vector<T>& lhs, const std::vector<T>& rhs){ \
        std::vector<T> ret(lhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=Func2(lhs[p], rhs[p]); } \
        return ret;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T> Func(const std::vector<T>& lhs, const rhsType& rhs){ \
        std::vector<T> ret(lhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=Func2(lhs[p], rhs); }    \
        return ret;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T> Func(const rhsType& lhs, const std::vector<T>& rhs){ \
        std::vector<T> ret(rhs.size());                                    \
        for(uint p=0; p<ret.size(); p++){ ret[p]=Func2(lhs, rhs[p]); }    \
        return ret;                                                        \
    }
#define SSTD_DEF_stdVecEx_f_eq(Func, Func2)                                \
    template <typename T>                                                \
    inline std::vector<T>& Func(std::vector<T>& lhs, const std::vector<T>& rhs){ \
        for(uint p=0; p<lhs.size(); p++){ lhs[p]=Func2(lhs[p], rhs[p]); } \
        return lhs;                                                        \
    }                                                                    \
    template <typename T, typename rhsType>                                \
    inline std::vector<T>& Func(std::vector<T>& lhs, const rhsType& rhs){ \
        for(uint p=0; p<lhs.size(); p++){ lhs[p]=Func2(lhs[p], rhs); }    \
        return lhs;                                                        \
    }
SSTD_DEF_stdVecEx_o   (sstd_stdVecEx::add   , + );
SSTD_DEF_stdVecEx_o_eq(sstd_stdVecEx::add_eq, +=);
SSTD_DEF_stdVecEx_o   (sstd_stdVecEx::sub   , - );
SSTD_DEF_stdVecEx_o_eq(sstd_stdVecEx::sub_eq, -=);
SSTD_DEF_stdVecEx_o   (sstd_stdVecEx::mul   , * );
SSTD_DEF_stdVecEx_o_eq(sstd_stdVecEx::mul_eq, *=);
SSTD_DEF_stdVecEx_o   (sstd_stdVecEx::div   , / );
SSTD_DEF_stdVecEx_o_eq(sstd_stdVecEx::div_eq, /=);
SSTD_DEF_stdVecEx_o   (sstd_stdVecEx::mod   , % );
SSTD_DEF_stdVecEx_o_eq(sstd_stdVecEx::mod_eq, %=);
SSTD_DEF_stdVecEx_f   (sstd_stdVecEx::pow   , sstd::pow); // ^
SSTD_DEF_stdVecEx_f_eq(sstd_stdVecEx::pow_eq, sstd::pow); // ^=
#undef SSTD_DEF_stdVecEx_o    // Deletion of used definition, in order not to pollute the namespace
#undef SSTD_DEF_stdVecEx_o_eq // Deletion of used definition, in order not to pollute the namespace
#undef SSTD_DEF_stdVecEx_f    // Deletion of used definition, in order not to pollute the namespace
#undef SSTD_DEF_stdVecEx_f_eq // Deletion of used definition, in order not to pollute the namespace

//---

// operators for std::vector
template <typename T>
inline std::vector<T> sstd_stdVecEx::push_back(const std::vector<T>& lhs, const std::vector<T>& rhs){
    std::vector<T> ret(lhs.size()+rhs.size());
    uint i=0;
    for(uint p=0; p<lhs.size(); p++){ ret[i]=lhs[p]; i++; }
    for(uint p=0; p<rhs.size(); p++){ ret[i]=lhs[p]; i++; }
    return ret;
}
template <typename T, typename rhsType>
inline std::vector<T> sstd_stdVecEx::push_back(const std::vector<T>& lhs, const rhsType& rhs){
    std::vector<T> ret(lhs.size()+1);
    for(uint p=0; p<lhs.size(); p++){ ret[p]=lhs[p]; }
    ret[lhs.size()]=rhs;
    return ret;
}
template <typename T, typename lhsType>
inline std::vector<T> sstd_stdVecEx::push_back(const lhsType& lhs, const std::vector<T>& rhs){
    std::vector<T> ret(rhs.size()+1);
    ret[0]=lhs;
    for(uint p=0; p<rhs.size(); p++){ ret[p+1]=rhs[p]; }
    return ret;
}

template <typename T>
inline std::vector<T>& sstd_stdVecEx::push_back_eq(std::vector<T>& lhs, const std::vector<T>& rhs){
    lhs.insert(lhs.end(), rhs.begin(), rhs.end());
    return lhs;
}
template <typename T, typename rhsType>
inline std::vector<T>& sstd_stdVecEx::push_back_eq(std::vector<T>& lhs, const rhsType& rhs){
    lhs.push_back(rhs);
    return lhs;
}

//-----------------------------------------------------------------------------------------------------------------------------------------------

#define SSTD_DEF_stdVecEx_Operator(Func, Ope)                            \
    template <typename T>                   inline std::vector<T> operator Ope(const std::vector<T>& lhs, const std::vector<T>& rhs){ return Func<T>         (lhs, rhs); } \
    template <typename T, typename rhsType> inline std::vector<T> operator Ope(const std::vector<T>& lhs, const        rhsType& rhs){ return Func<T, rhsType>(lhs, rhs); } \
    template <typename T, typename lhsType> inline std::vector<T> operator Ope(const        lhsType& lhs, const std::vector<T>& rhs){ return Func<T, lhsType>(lhs, rhs); }
#define SSTD_DEF_stdVecEx_Operator_eq(Func, Ope)                        \
    template <typename T>                   inline std::vector<T>& operator Ope(std::vector<T>& lhs, const std::vector<T>& rhs){ return Func<T>         (lhs, rhs); } \
    template <typename T, typename rhsType> inline std::vector<T>& operator Ope(std::vector<T>& lhs, const        rhsType& rhs){ return Func<T, rhsType>(lhs, rhs); }

// operators for mathematics
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::add   , + );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::add_eq, +=);
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::sub   , - );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::sub_eq, -=);
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::mul   , * );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::mul_eq, *=);
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::div   , / );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::div_eq, /=);
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::mod   , % );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::mod_eq, %=);
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::pow   , ^ );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::pow_eq, ^=);

// operators for std::vector
SSTD_DEF_stdVecEx_Operator   (sstd_stdVecEx::push_back   , << );
SSTD_DEF_stdVecEx_Operator_eq(sstd_stdVecEx::push_back_eq, <<=);

#undef SSTD_DEF_stdVecEx_Operator    // Deletion of used definition, in order not to pollute the namespace
#undef SSTD_DEF_stdVecEx_Operator_eq // Deletion of used definition, in order not to pollute the namespace

//---

template <typename T> inline std::vector<T>& operator++(std::vector<T>& rhs)     { for(uint p=0; p<rhs.size(); p++){ rhs[p]++; } return rhs; } // ++rhs
template <typename T> inline std::vector<T>& operator++(std::vector<T>& rhs, int){ for(uint p=0; p<rhs.size(); p++){ rhs[p]++; } return rhs; } //   rhs++
template <typename T> inline std::vector<T>& operator--(std::vector<T>& rhs)     { for(uint p=0; p<rhs.size(); p++){ rhs[p]--; } return rhs; } // --rhs
template <typename T> inline std::vector<T>& operator--(std::vector<T>& rhs, int){ for(uint p=0; p<rhs.size(); p++){ rhs[p]--; } return rhs; } //   rhs--

//-----------------------------------------------------------------------------------------------------------------------------------------------

Results
  実際に,上記で示した実装を実行すると,
  ./main.cpp
#include <sstd/sstd.hpp>

int main(){
    std::vector<double> a={1, 2, 3}, b={4, 5, 6};
    sstd::printn(a);
    a*=4;                 // same as "for(int i=0; i<a.size(); i++){ a[i]*=4; }"
    sstd::printn(a);
    printf("\n");
    
    a<<=b;                // same as "a.insert(a.end(), b.begin(), b.end());".
    sstd::printn(a);
    a<<=7.0;              // same as "a.push_back(7.0);".
    sstd::printn(a);
    sstd::printn(a<<8.0); // same as "{ std::vector<double> tmp(a.size()+1); for(uint p=0;p<a.size();p++){tmp[p]=a[p];} tmp[a.size()]=8.0; sstd::printn(tmp); }"
    
    return 0;
}
  Execution result
$ ./exe 
a[3] = [ 1.000000 2.000000 3.000000 ]
a[3] = [ 4.000000 8.000000 12.000000 ]

a[6] = [ 4.000000 8.000000 12.000000 4.000000 5.000000 6.000000 ]
a[7] = [ 4.000000 8.000000 12.000000 4.000000 5.000000 6.000000 7.000000 ]
a<<8.0[8] = [ 4.000000 8.000000 12.000000 4.000000 5.000000 6.000000 7.000000 8.000000 ]
  Running with Ideone.com
となる.なお,本ページに表示されている実行結果は,拙ライブラリ sstd を呼び出した実行結果であり,リンク先の Ideone.com にベタ打ちされた実装の結果と,ほぼ同様である.例には示していないが,当然,a<<b; や,a*b; などの演算はすべて実装されている.気になる演算があれば,Ideone.com のリンクから fork し,動作を確認してもらいたい.
  使用方法としては,sstd を利用するか,上記の Ideone.com の内容をヘッダにコピーし,インクルードするとよい.
Conclusion
  テンプレートクラス std::vector<T> に,演算子 +, +=, -, -=, *, *=, /, /=, %, %=, ^, ^=, <<, <<= を定義し,MATLAB や numpy のような配列計算を可能とした.
References
  [追記1] RAII - Wikipedia - 2019年01月05日閲覧
Future work
  拙ライブラリ sstd では,同様の演算子を,列優先 (col-major) の行列型 sstd::mat<T> に対して実装している.同様の実装を,行優先 (row-major) の行列型 sstd::mat_r<T> についても行うことが今後の予定となる.

2018.12.06 追記:
  ・push_front(): >>= の実装.
  ・std::move() などにより,fn(T&& rhs) を受け取った場合,swap() によりコピーを削減する.
2019.08.09 追記:
  ・push_front(): >>= の実装完了.
Source code
  本投稿で使用したソースコードを下記に示す.


また,今後変更があった場合は,上記ではなく,下記に示す sstd 中の ./sstd/src/stdVector_expansion/stdVector_expansion.hpp に反映される.


0 件のコメント:

コメントを投稿