《奇異遞迴模板模式（CRTP）應用--表示式模板(expression template) 2》
奇異遞迴模板模式（CRTP）應用--表示式模板(expression template) 2
DanielWang_2017-07-12 23:30:423869收藏3 分類專欄：c/c++/cpp11 版權
1 表示式模板（expression template）概述
首先分幾個部分介紹下expression template。
1.1 表示式模板（expression template）是什麼？
引用wiki 介紹的Expression templates：
Expression templates are a C++ template metaprogramming technique that builds structures representing a computation at compile time, which structures are evaluated only as needed to produce efficient code for the entire computation. Expression templates thus allow programmers to bypass the normal order of evaluation of the C++ language and achieve optimizations such as loop fusion.
Expression templates were invented independently by Todd Veldhuizen and David Vandevoorde; it was Veldhuizen who gave them their name.They are a popular technique for the implementation of linear algebra software.
結合上面的介紹，簡單總結下expression template，表示式模板（expression-template）是一種C++模板超程式設計技術，它在編譯時刻構建好一些能夠表達某一計算的結構，對於整個計算這些結構僅在
1.2 為什麼引入 expression template
先通過一個例子來說明，假設在機器學習中我們現在需要構造一個類，這個類有size_ 維feature,每個feature的值是個double型別的實數，這個類不同物件間的同一維度的feature可以進行簡單的算術運算（比如+、- maximum…），比較簡單的思想是進行運算子過載，具體見程式碼：
// a naive solution( bad! )
// operator overloading
#include <iostream>
#include <algorithm>
#include <cassert>
using namespace std;
class MyType {
public:
MyType(size_t size, double val = 0.0) : size_(size) {
features_ = new double [size_];
for (size_t i = 0; i < size_; i++) features_[i] = val;
cout << "\tMyType constructor size = " << size_ << "\n";
}
// construct vector using initializer list
MyType(std::initializer_list<double> init) {
size_ = init.size();
features_ = new double[size_];
for (size_t i = 0; i < size_; i++) features_[i] = *(init.begin() + i);
cout << "\tMyType constructor size = " << size_ << "\n";
}
// in c++11 move constructor can be used here
MyType(const MyType& rhs):size_(rhs.size()) {
features_ = new double[size_];
for (size_t i = 0; i < size_; i++) features_[i] = rhs[i];
cout << "\tMyType copy constructor size = " << size_ << "\n";
}
MyType& operator=(const MyType& rhs) {
if (this != &rhs) {
delete[]features_;
size_ = rhs.size();
features_ = new double[size_];
for (size_t i = 0; i < size_; i++) features_[i] = rhs[i];
cout << "\tMyType assignment constructor size = " << size_ << "\n";
}
return *this;
}
~MyType() {
if (nullptr != features_) {
delete [] features_;
size_ = 0;
features_ = nullptr;
}
}
double &operator[](size_t i) { return features_[i]; }
double operator[](size_t i) const { return features_[i]; }
size_t size() const { return size_; }
private:
size_t size_;
double* features_;
};
//This kind of approach is inefficient, because temporal memory is allocated and de-allocated during each operation
MyType operator+(MyType const &u, MyType const &v) {
MyType sum(std::max(u.size(), v.size()));
for (size_t i = 0; i < u.size(); i++) {
sum[i] = u[i] + v[i];
}
cout << "\t in MyType +\n";
return sum;
}
// operator- balabala
int main() {
MyType a(3, 1.1);
MyType b(3, 2.01);
MyType c = {3.01, 3.01, 3.01};
cout << "\t----computing-----\n";
MyType d = a + b + c;
for (size_t i = 0; i < d.size(); i++) {
cout << "\t" << d[i] << "\n";
}
return 0;
}- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
因為我們這個類涉及資源申請，因此需要實現普通建構函式、拷貝構造、賦值構造、析構（其實C++11裡實現移動構造、移動賦值更佳,後面會code說明）。執行結果如下：
對於實現d=a+b+c這樣的多個連加的表示式，我們可以看到它進行多次物件的構造（物件a、b、c進行三次普通物件構造（必需），然後按照+的順序進行運算，按照計算順序，表示式其實是這樣的d=((a+b)+c))，首選a+b相加會構造一個臨時物件，a、b相加的結果(這裡假設a+b=tmp1)返會又要進行拷貝構造，然後表示式實際變成d=tmp1+c, tmp1和c相加又會構造一個臨時物件（這裡假設為tmp2,即有tmp1+c=tmp2），然後tmp2拷貝構造生成最終的結果d，細數下簡單的連加進行了很多物件的構造析構（其實就是很多次記憶體的申請和釋放，這個效率不要太低！！！）。彷彿聽到有童鞋高呼我們可以用C++11的移動語義，來減少臨時物件（記憶體的頻繁申請、釋放）頻繁記憶體操作，好先上程式碼：
// move constructor
MyType(MyType&& rhs) :size_(rhs.size()) {
features_ = rhs.features_;
rhs.size_ = 0;
rhs.features_ = nullptr;
cout << "\tMyType move constructor size = " << size_ << "\n";
}
// move assignment
MyType& operator=(MyType&& rhs) {
if (this != &rhs) {
size_ = rhs.size_;
rhs.size_ = 0;
features_ = rhs.features_;
rhs.features_ = nullptr;
cout << "\tMyType move assignment constructor size = " << size_ << "\n";
}
*this;
}- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
即是利用C++11的移動語義(move constructor 、assignment constructor )來減少記憶體的重複申請與釋放，然而對於d = a + b + c 這樣連加的case至少需要有一次臨時記憶體的申請與釋放和兩次加操作的迴圈（如下圖，標號1的記憶體即為a+b=tmp1構造過程，物件tmp1存放臨時結果，當tmp1與c相加即有tmp1+c=tmp2，再生成tmp2後tmp1便釋放記憶體，tmp2用移動構造出連加結果d），因此，下面不得不引入一個大殺器 expression template。
注意上面幾張圖的結果是在vs2015上編譯的結果，可見沒進行NRVO(named return value optimization)，進行了RVO即Return Value Optimization，是一種編譯器優化技術，可以把通過函式返回建立的臨時物件給”去掉”，然後可以達到少呼叫拷貝構造的操作。號外，這篇有介紹RVO和std::move的原理的文章，可以看下。所以採用不同的編譯器以及優化等級可能結果和上圖不同, g++ 預設是開啟了 NRVO,可以加上這個編譯引數-fno-elide-constructors 去掉NRVO，詳細的測試可以參見，下面給出的github地址。
1.3 expression template 想解決的問題
引用wikiMore C++ Idioms/Expression-template所說的Express Template的意圖：
Intent
- To create a domain-specific embedded language (DSEL) in C++.
- To support lazy evaluation of C++ expressions (e.g., mathematical expressions), which can be executed much later in the program from the point of their definition.
- To pass an expression – not the result of the expression – as a parameter to a function.
總結來說其實Expression Template想要解決的問題是：
1. 在C++中建立一個嵌入式領域專用語言(DSEL)（大概意思就是在一個程式語言裡嵌入一個領域特定語言。參考這個 Domain Specific Languages in C++）。
2. 支援表示式的懶惰計算（延遲計算），相對其定義它可以推遲表示式的計算。
3. 傳遞一個表示式，不是表示式的結果而是把表示式自身作為一個引數傳給一個函式。
1.4 expression template的方法
結合1.2中的naive operator overloading下面舉例說明表示式模板的神奇功效，先show程式碼（對CRTP技術不瞭解的請先看我之前寫的奇異遞迴模板模式（ Curiously Recurring Template Pattern，CRTP）1）：
// expression template
// Copyright (C) 2017 http://ustcdane.github.io/
#include <algorithm>
#include <cstdio>
#include <cassert>
// this is expression, all expressions must inherit it
template<typename A>
struct Expression {
// returns const reference of the actual type of this expression
const A& cast() const { return static_cast<const A&>(*this); }
int size() const { return cast().size(); }// get Expression size
private:
Expression& operator=(const Expression&) {
return *this;
}
Expression() {}
friend A;
};
// binary expression: (binary op,lhs, rhs)
// Func can easily allows user customized theirs binary operators
template<typename Func, typename TLhs, typename TRhs>
struct BinaryOp : public Expression<BinaryOp<Func, TLhs, TRhs> > {
BinaryOp(const TLhs& lhs, const TRhs& rhs):lhs_(lhs.cast()), rhs_(rhs.cast()) {}
// work function, computing this expression at position i, import for lazy computing
double value(int i) const {
return Func::Op(lhs_.value(i), rhs_.value(i));
}
int size() const { return std::max(lhs_.size(), rhs_.size()); }
private:
const TLhs& lhs_;
const TRhs& rhs_;
};
// template binary operation, works for any expressions
template<typename Func, typename TLhs, typename TRhs>
inline BinaryOp<Func, TLhs, TRhs>
expToBinaryOp(const Expression<TLhs>& lhs, const Expression<TRhs>& rhs) {
return BinaryOp<Func, TLhs, TRhs>(lhs.cast(), rhs.cast());
}
// binary operators +
struct Add {
// any function defined inside its class definition is inline
static double Op(double a, double b) {
return a + b;
}
};
// define Minimum
struct Minimum {
static double Op(double a, double b) {
return a < b ? a : b;
}
};
template<typename TLhs, typename TRhs>
inline BinaryOp<Add, TLhs, TRhs>
operator+(const Expression<TLhs>& lhs, const Expression<TRhs>& rhs) {
return expToBinaryOp<Add>(lhs, rhs);
}
template<typename TLhs, typename TRhs>
inline BinaryOp<Minimum, TLhs, TRhs>
min(const Expression<TLhs>& lhs, const Expression<TRhs>& rhs) {
return expToBinaryOp<Minimum>(lhs, rhs);
}
// allocation just by user
// no constructor and destructor to allocate and de-allocate memory
class MyExpType : public Expression<MyExpType> {
public:
MyExpType():size_(0) {}
MyExpType(double *features, int size)
: size_(size), features_(features) {
printf("MyExpType constructor size = %d. No memory allocate.\n", size_);
}
// delete copy constructor,
MyExpType(const MyExpType& src_) = delete;
template<typename ExpType>
MyExpType(const Expression<ExpType>& src_) = delete;
// here is where computing happens,lazy support
template<typename ExpType>
MyExpType& operator=(const Expression<ExpType>& src) {
const ExpType &srcReal = src.cast();
assert(size_ >= srcReal.size()); //
for (int i = 0; i < srcReal.size(); ++i) {
features_[i] = srcReal.value(i); // binary expression value work function
}
printf("MyExpType assignment constructor size = %d\n", size_);
return *this;
}
// computing function
double value(int i) const { return features_[i]; }
int size() const { return size_; }
private:
int size_;
double* features_;
};
void print(const MyExpType& m) {
printf("( ");
for (int i = 0; i < m.size() - 1; ++i) {
printf("%g, ", m.value(i));
}
if (m.size())
printf("%g )\n", m.value(m.size()-1));
else
printf(" )\n");
}
int main() {
const int N = 3;
double sa[N] = { 1.1,1.1, 1.1 };
double sb[N] = { 2.01, 2.01, 2.01 };
double sc[N] = { 3.01, 3.01, 3.01 };
double sd[N] = { 0 };
MyExpType A(sa, N), B(sb, N), C(sc, N), D(sd, N);
printf("\n");
printf(" A = "); print(A);
printf(" B = "); print(B);
printf(" C = "); print(C);
printf("\n\tD = A + B + C\n");
D = A + B + C;
for (int i = 0; i < A.size(); ++i) {
printf("%d:\t%g + %g + %g = %g\n",
i, B.value(i), C.value(i), B.value(i), D.value(i));
}
printf("\n\tD = A + min(B, C)\n");
// D = A + expToBinaryOp<Minimum>(B, C);
D = A + min(B, C);
for (int i = 0; i < A.size(); ++i) {
printf("%d:\t%g + min(%g, %g) = %g\n",
i, A.value(i), B.value(i), C.value(i), D.value(i));
}
return 0;
}- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
程式執行結果如下圖：
如上圖所示，看到expression template的神奇了吧，我們的類這下沒有記憶體申請、釋放（使用者給一個有效地址記憶體及size即可），三連加只有一次賦值構造，無記憶體申請、釋放，而且只有一次迴圈。正像上一小節所說的，Expression template所要解決是延遲計算問題（lazy evaluation ）。它可以在C++中讓operator+返回自定義型別來實現這一功能。當表示式比較長時可以有效的構建出表示式樹如分析樹(parse tree)這樣的形式。對應到本程式的表示式 D=A+B+C ，其對應的表達（型別）樹如下圖所示，表示式的計算只在賦值（operator=）給一個實際物件時才進行。表示式模板通過在編譯期的表示式樹來實現延遲計算。
表示式右側被解析為型別expToBinaryOp<expToBinaryOp<MyExpType, MyExpType>, MyExpType > 這裡記為RHS，如表示式樹所示，表示式樹結點的型別是在編譯期(並且code inlining)從底向上傳遞的。當呼叫型別MyExpType 賦值運算子時，這裡只展開一次迴圈,即在位置i處會呼叫RHS.value(i)，型別RHS會根據operator+ 和value(i)的定義遞迴的把RHS.value(i)展開為 RHS.value(i) = A.value(i) + B.value(i) + C.value(i),可以看到表示式的神奇，無記憶體申請釋放，多次迴圈變成一次迴圈，即如下程式碼：
for (int i = 0; i < D.size(); ++i) {
D.features_[i] = A.value(i) + B.value(i) + C.value(i);
} - 1
- 2
- 3
2. 表示式模板demo
結合表示式模板的相關分析，實現了一個支援一元及二元自定義型別的表示式模板操作庫（當然是用了CRTP技術），來段測試程式碼先：
#include <iostream>
#include <cmath>
#include "valueType.hpp"
using namespace std;
using op::doubleType;
doubleType algorithm(const doubleType x[2]) {
doubleType y = 1;
y += x[1] + x[0] * x[1];
cout << "doubleType y = 1 + x[1] + x[0]*x[1]\n";
return y;
}
doubleType algorithm2(const doubleType x[2]) {
doubleType y = exp(x[0]) + -x[1] + 1;
cout << "doubleType y = exp(x[0]) + -x[1] + 1\n";
return y;
}
doubleType algorithm3(const doubleType x[2]) {
const double PI = 3.141592653589793;
doubleType y = cos(x[0]) + sin(PI/2*log2(x[1]));
cout << "cos(x[0]) + sin(PI/2*log2(x[1]))\n";
return y;
}
doubleType algorithm4(const doubleType x[2]) {
const doubleType m = -0.618;
doubleType y = min(m, min(x[0], x[1]));
cout << "m=-0.618\nmin(m, min(x[0], x[1]))\n";
return y;
}
int main(int argc, char** argv)
{
doubleType x[2];
doubleType y;
x[0] = 0.0;
x[1] = 2.0;
cout << "x[0] = " << x[0] << endl;
cout << "x[1] = " << x[1] << endl;
std::cout << "Begin run algorithm:\n";
y = algorithm(x);
std::cout << "y = " << y.value() << "\n";
y = algorithm2(x);
std::cout << "y = " << value(y) << "\n";
y = algorithm3(x);
std::cout << "y = " << y << "\n";
if (2 == y)
cout << "Yes.\ty = 2\n";
y = algorithm4(x);
std::cout << "y = " << value(y) << "\n";
return 0;
}- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
執行結果：
分析：如程式碼中 algorithm3 y=cos(x[0]) +sin(PI/2*log2(x[1]))，其對應的表示式型別樹如下圖所示：
具體為什麼是這樣的可以參考程式碼，其中ScalarMul是表示式型別和標量型別相乘的表示式型別。本文所有程式碼請到我的github。
3. Expression Template總結
expression template能讓程式碼高效很多，即使現在最新的cpp11移動語義也不能達到我們所說的表示式模板連加的效率，尤其在代數計算領域，表示式模板更是實現相關工具的典範，有機會的話再寫篇關於自動求導(automatic differentiation,AD)的文章吧(好像介紹CRTP技術的時候就說要寫了…)，也是利用了表示式模板的技術。