1// Copyright (C) 2016-2018 T. Zachary Laine
2//
3// Distributed under the Boost Software License, Version 1.0. (See
4// accompanying file LICENSE_1_0.txt or copy at
5// http://www.boost.org/LICENSE_1_0.txt)
6//[ self_evaluation
7#include <boost/yap/expression.hpp>
8
9#include <boost/optional.hpp>
10#include <boost/hana/fold.hpp>
11#include <boost/hana/maximum.hpp>
12
13#include <algorithm>
14#include <cassert>
15#include <iostream>
16#include <vector>
17
18
19// A super-basic matrix type, and a few associated operations.
20struct matrix
21{
22 matrix() : values_(), rows_(0), cols_(0) {}
23
24 matrix(int rows, int cols) : values_(rows * cols), rows_(rows), cols_(cols)
25 {
26 assert(0 < rows);
27 assert(0 < cols);
28 }
29
30 int rows() const { return rows_; }
31 int cols() const { return cols_; }
32
33 double operator()(int r, int c) const
34 { return values_[r * cols_ + c]; }
35 double & operator()(int r, int c)
36 { return values_[r * cols_ + c]; }
37
38private:
39 std::vector<double> values_;
40 int rows_;
41 int cols_;
42};
43
44matrix operator*(matrix const & lhs, double x)
45{
46 matrix retval = lhs;
47 for (int i = 0; i < retval.rows(); ++i) {
48 for (int j = 0; j < retval.cols(); ++j) {
49 retval(i, j) *= x;
50 }
51 }
52 return retval;
53}
54matrix operator*(double x, matrix const & lhs) { return lhs * x; }
55
56matrix operator+(matrix const & lhs, matrix const & rhs)
57{
58 assert(lhs.rows() == rhs.rows());
59 assert(lhs.cols() == rhs.cols());
60 matrix retval = lhs;
61 for (int i = 0; i < retval.rows(); ++i) {
62 for (int j = 0; j < retval.cols(); ++j) {
63 retval(i, j) += rhs(i, j);
64 }
65 }
66 return retval;
67}
68
69// daxpy() means Double-precision AX Plus Y. This crazy name comes from BLAS.
70// It is more efficient than a naive implementation, because it does not
71// create temporaries. The covnention of using Y as an out-parameter comes
72// from FORTRAN BLAS.
73matrix & daxpy(double a, matrix const & x, matrix & y)
74{
75 assert(x.rows() == y.rows());
76 assert(x.cols() == y.cols());
77 for (int i = 0; i < y.rows(); ++i) {
78 for (int j = 0; j < y.cols(); ++j) {
79 y(i, j) += a * x(i, j);
80 }
81 }
82 return y;
83}
84
85template <boost::yap::expr_kind Kind, typename Tuple>
86struct self_evaluating_expr;
87
88template <boost::yap::expr_kind Kind, typename Tuple>
89auto evaluate_matrix_expr(self_evaluating_expr<Kind, Tuple> const & expr);
90
91// This is the primary template for our expression template. If you assign a
92// self_evaluating_expr to a matrix, its conversion operator transforms and
93// evaluates the expression with a call to evaluate_matrix_expr().
94template <boost::yap::expr_kind Kind, typename Tuple>
95struct self_evaluating_expr
96{
97 operator auto() const;
98
99 static const boost::yap::expr_kind kind = Kind;
100
101 Tuple elements;
102};
103
104// This is a specialization of our expression template for assignment
105// expressions. The destructor transforms and evaluates via a call to
106// evaluate_matrix_expr(), and then assigns the result to the variable on the
107// left side of the assignment.
108//
109// In a production implementation, you'd need to have specializations for
110// plus_assign, minus_assign, etc.
111template <typename Tuple>
112struct self_evaluating_expr<boost::yap::expr_kind::assign, Tuple>
113{
114 ~self_evaluating_expr();
115
116 static const boost::yap::expr_kind kind = boost::yap::expr_kind::assign;
117
118 Tuple elements;
119};
120
121struct use_daxpy
122{
123 // A plus-expression, which may be of the form double * matrix + matrix,
124 // or may be something else. Since our daxpy() above requires a mutable
125 // "y", we only need to match a mutable lvalue matrix reference here.
126 template <typename Tuple>
127 auto operator()(
128 boost::yap::expr_tag<boost::yap::expr_kind::plus>,
129 self_evaluating_expr<boost::yap::expr_kind::multiplies, Tuple> const & expr,
130 matrix & m)
131 {
132 // Here, we transform the left-hand side into a pair if it's the
133 // double * matrix operation we're looking for. Otherwise, we just
134 // get a copy of the left side expression.
135 //
136 // Note that this is a bit of a cheat, done for clarity. If we pass a
137 // larger expression that happens to contain a double * matrix
138 // subexpression, that subexpression will be transformed into a tuple!
139 // In production code, this transform should probably only be
140 // performed on an expression with all terminal members.
141 auto lhs = boost::yap::transform(
142 expr,
143 [](boost::yap::expr_tag<boost::yap::expr_kind::multiplies>,
144 double d,
145 matrix const & m) {
146 return std::pair<double, matrix const &>(d, m);
147 });
148
149 // If we got back a copy of expr above, just re-construct the
150 // expression this function mathes; in other words, do not effectively
151 // transform anything. Otherwise, replace the expression matched by
152 // this function with a call to daxpy().
153 if constexpr (boost::yap::is_expr<decltype(lhs)>::value) {
154 return expr + m;
155 } else {
156 return boost::yap::make_terminal(t&: daxpy)(lhs.first, lhs.second, m);
157 }
158 }
159};
160
161
162// This is the heart of what self_evaluating_expr does. If we had other
163// optimizations/transformations we wanted to do, we'd put them in this
164// function, either before or after the use_daxpy transformation.
165template <boost::yap::expr_kind Kind, typename Tuple>
166auto evaluate_matrix_expr(self_evaluating_expr<Kind, Tuple> const & expr)
167{
168 auto daxpy_form = boost::yap::transform(expr, use_daxpy{});
169 return boost::yap::evaluate(daxpy_form);
170}
171
172template<boost::yap::expr_kind Kind, typename Tuple>
173self_evaluating_expr<Kind, Tuple>::operator auto() const
174{
175 return evaluate_matrix_expr(*this);
176}
177
178template<typename Tuple>
179self_evaluating_expr<boost::yap::expr_kind::assign, Tuple>::
180 ~self_evaluating_expr()
181{
182 using namespace boost::hana::literals;
183 boost::yap::evaluate(elements[0_c]) = evaluate_matrix_expr(elements[1_c]);
184}
185
186// In order to define the = operator with the semantics we want, it's
187// convenient to derive a terminal type from a terminal instantiation of
188// self_evaluating_expr. Note that we could have written a template
189// specialization here instead -- either one would work. That would of course
190// have required more typing.
191struct self_evaluating :
192 self_evaluating_expr<
193 boost::yap::expr_kind::terminal,
194 boost::hana::tuple<matrix>
195 >
196{
197 self_evaluating() {}
198
199 explicit self_evaluating(matrix m)
200 { elements = boost::hana::tuple<matrix>(std::move(m)); }
201
202 BOOST_YAP_USER_ASSIGN_OPERATOR(self_evaluating_expr, ::self_evaluating_expr);
203};
204
205BOOST_YAP_USER_BINARY_OPERATOR(plus, self_evaluating_expr, self_evaluating_expr)
206BOOST_YAP_USER_BINARY_OPERATOR(minus, self_evaluating_expr, self_evaluating_expr)
207BOOST_YAP_USER_BINARY_OPERATOR(multiplies, self_evaluating_expr, self_evaluating_expr)
208
209
210int main()
211{
212 matrix identity(2, 2);
213 identity(0, 0) = 1.0;
214 identity(1, 1) = 1.0;
215
216 // These are YAP-ified terminal expressions.
217 self_evaluating m1(identity);
218 self_evaluating m2(identity);
219 self_evaluating m3(identity);
220
221 // This transforms the YAP expression to use daxpy(), so it creates no
222 // temporaries. The transform happens in the destructor of the
223 // assignment-expression specialization of self_evaluating_expr.
224 m1 = 3.0 * m2 + m3;
225
226 // Same as above, except that it uses the matrix conversion operator on
227 // the self_evaluating_expr primary template, because here we're assigning
228 // a YAP expression to a non-YAP-ified matrix.
229 matrix m_result_1 = 3.0 * m2 + m3;
230
231 // Creates temporaries and does not use daxpy(), because the A * X + Y
232 // pattern does not occur within the expression.
233 matrix m_result_2 = 3.0 * m2;
234
235 return 0;
236}
237//]
238

source code of boost/libs/yap/example/self_evaluation.cpp