vector_expression.hpp
Go to the documentation of this file.
1/*!
2 * \brief expression templates for vector valued math
3 *
4 * \author O. Krause
5 * \date 2013
6 *
7 *
8 * \par Copyright 1995-2015 Shark Development Team
9 *
10 * <BR><HR>
11 * This file is part of Shark.
12 * <http://image.diku.dk/shark/>
13 *
14 * Shark is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License as published
16 * by the Free Software Foundation, either version 3 of the License, or
17 * (at your option) any later version.
18 *
19 * Shark is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Lesser General Public License for more details.
23 *
24 * You should have received a copy of the GNU Lesser General Public License
25 * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26 *
27 */
28#ifndef REMORA_VECTOR_EXPRESSION_HPP
29#define REMORA_VECTOR_EXPRESSION_HPP
30
31#include "detail/expression_optimizers.hpp"
32#include "kernels/dot.hpp"
33#include "kernels/vector_fold.hpp"
34#include "kernels/vector_max.hpp"
35
36namespace remora{
37
38template<class T, class VecV, class Device>
39typename std::enable_if<
40 std::is_convertible<T, typename VecV::value_type >::value,
41 typename detail::vector_scalar_multiply_optimizer<VecV>::type
42>::type
43operator* (vector_expression<VecV, Device> const& v, T scalar){
44 typedef typename VecV::value_type value_type;
45 return detail::vector_scalar_multiply_optimizer<VecV>::create(v(),value_type(scalar));
46}
47template<class T, class VecV, class Device>
48typename std::enable_if<
49 std::is_convertible<T, typename VecV::value_type >::value,
50 typename detail::vector_scalar_multiply_optimizer<VecV>::type
51>::type
52operator* (T scalar, vector_expression<VecV, Device> const& v){
53 typedef typename VecV::value_type value_type;
54 return detail::vector_scalar_multiply_optimizer<VecV>::create(v(),value_type(scalar));
55}
56
57template<class VecV, class Device>
58typename detail::vector_scalar_multiply_optimizer<VecV>::type
59operator-(vector_expression<VecV, Device> const& v){
60 typedef typename VecV::value_type value_type;
61 return detail::vector_scalar_multiply_optimizer<VecV>::create(v(),value_type(-1));
62}
63
64///\brief Creates a vector having a constant value.
65///
66///@param scalar the value which is repeated
67///@param elements the size of the resulting vector
68template<class T>
69typename std::enable_if<std::is_arithmetic<T>::value, scalar_vector<T, cpu_tag> >::type
70repeat(T scalar, std::size_t elements){
71 return scalar_vector<T, cpu_tag>(elements,scalar);
72}
73
74
75#define REMORA_UNARY_VECTOR_TRANSFORMATION(name, F)\
76template<class VecV, class Device>\
77typename detail::vector_unary_optimizer<VecV,typename device_traits<Device>:: template F<typename VecV::value_type> >::type \
78name(vector_expression<VecV, Device> const& v){\
79 typedef typename device_traits<Device>:: template F<typename VecV::value_type> functor_type;\
80 return detail::vector_unary_optimizer<VecV, functor_type >::create(v(), functor_type());\
81}
82REMORA_UNARY_VECTOR_TRANSFORMATION(abs, abs)
83REMORA_UNARY_VECTOR_TRANSFORMATION(log, log)
84REMORA_UNARY_VECTOR_TRANSFORMATION(exp, exp)
85REMORA_UNARY_VECTOR_TRANSFORMATION(tanh,tanh)
86REMORA_UNARY_VECTOR_TRANSFORMATION(sin,sin)
87REMORA_UNARY_VECTOR_TRANSFORMATION(cos,cos)
88REMORA_UNARY_VECTOR_TRANSFORMATION(tan,tan)
89REMORA_UNARY_VECTOR_TRANSFORMATION(asin,asin)
90REMORA_UNARY_VECTOR_TRANSFORMATION(acos,acos)
91REMORA_UNARY_VECTOR_TRANSFORMATION(atan,atan)
92REMORA_UNARY_VECTOR_TRANSFORMATION(erf,erf)
93REMORA_UNARY_VECTOR_TRANSFORMATION(erfc,erfc)
94REMORA_UNARY_VECTOR_TRANSFORMATION(sqr, sqr)
95REMORA_UNARY_VECTOR_TRANSFORMATION(sqrt, sqrt)
96REMORA_UNARY_VECTOR_TRANSFORMATION(cbrt, cbrt)
97REMORA_UNARY_VECTOR_TRANSFORMATION(sigmoid, sigmoid)
98REMORA_UNARY_VECTOR_TRANSFORMATION(softPlus, soft_plus)
99REMORA_UNARY_VECTOR_TRANSFORMATION(elem_inv, inv)
100#undef REMORA_UNARY_VECTOR_TRANSFORMATION
101
102///\brief Adds two vectors
103template<class VecV1, class VecV2, class Device>
104vector_addition<VecV1, VecV2 > operator+ (
105 vector_expression<VecV1, Device> const& v1,
106 vector_expression<VecV2, Device> const& v2
107){
108 REMORA_SIZE_CHECK(v1().size() == v2().size());
109 return vector_addition<VecV1, VecV2>(v1(),v2());
110}
111///\brief Subtracts two vectors
112template<class VecV1, class VecV2, class Device>
113auto operator- (
114 vector_expression<VecV1, Device> const& v1,
115 vector_expression<VecV2, Device> const& v2
116) -> decltype (v1 + (-v2)){
117 REMORA_SIZE_CHECK(v1().size() == v2().size());
118 return v1 + (-v2);
119}
120
121///\brief Adds a vector plus a scalar which is interpreted as a constant vector
122template<class VecV, class T, class Device>
123typename std::enable_if<
124 std::is_convertible<T, typename VecV::value_type>::value,
125 vector_addition<VecV, scalar_vector<T, Device> >
126>::type operator+ (
127 vector_expression<VecV, Device> const& v,
128 T t
129){
130 return v + scalar_vector<T, Device>(v().size(),t);
131}
132
133///\brief Adds a vector plus a scalar which is interpreted as a constant vector
134template<class T, class VecV, class Device>
135typename std::enable_if<
136 std::is_convertible<T, typename VecV::value_type>::value,
137 vector_addition<VecV, scalar_vector<T, Device> >
138>::type operator+ (
139 T t,
140 vector_expression<VecV, Device> const& v
141){
142 return v + scalar_vector<T, Device>(v().size(),t);
143}
144
145///\brief Subtracts a scalar which is interpreted as a constant vector.
146template<class VecV, class T, class Device>
147typename std::enable_if<
148 std::is_convertible<T, typename VecV::value_type>::value ,
149 decltype(std::declval<VecV const&>() + T())
150>::type operator- (
151 vector_expression<VecV, Device> const& v,
152 T t
153){
154 return v + (-t);
155}
156
157///\brief Subtracts a vector from a scalar which is interpreted as a constant vector
158template<class VecV, class T, class Device>
159typename std::enable_if<
160 std::is_convertible<T, typename VecV::value_type>::value,
161 decltype(T() + (-std::declval<VecV const&>()))
162>::type operator- (
163 T t,
164 vector_expression<VecV, Device> const& v
165){
166 return t + (-v);
167}
168
169#define REMORA_BINARY_VECTOR_EXPRESSION(name, F)\
170template<class VecV1, class VecV2, class Device>\
171vector_binary<VecV1, VecV2, typename device_traits<Device>:: template F<typename common_value_type<VecV1,VecV2>::type> >\
172name(vector_expression<VecV1, Device> const& v1, vector_expression<VecV2, Device> const& v2){\
173 REMORA_SIZE_CHECK(v1().size() == v2().size());\
174 typedef typename common_value_type<VecV1,VecV2>::type type;\
175 typedef typename device_traits<Device>:: template F<type> functor_type;\
176 return vector_binary<VecV1, VecV2, functor_type >(v1(),v2(), functor_type());\
177}
178REMORA_BINARY_VECTOR_EXPRESSION(operator*, multiply)
179REMORA_BINARY_VECTOR_EXPRESSION(element_prod, multiply)
180REMORA_BINARY_VECTOR_EXPRESSION(operator/, divide)
181REMORA_BINARY_VECTOR_EXPRESSION(element_div, divide)
182REMORA_BINARY_VECTOR_EXPRESSION(pow, pow)
183REMORA_BINARY_VECTOR_EXPRESSION(min, min)
184REMORA_BINARY_VECTOR_EXPRESSION(max, max)
185#undef REMORA_BINARY_VECTOR_EXPRESSION
186
187
188//operations of the form op(v,t)[i] = op(v[i],t)
189#define REMORA_VECTOR_SCALAR_TRANSFORMATION(name, F)\
190template<class T, class VecV, class Device> \
191typename std::enable_if< \
192 std::is_convertible<T, typename VecV::value_type >::value,\
193 vector_binary<VecV, scalar_vector<typename VecV::value_type, Device>, \
194 typename device_traits<Device>:: template F<typename VecV::value_type> > \
195>::type \
196name (vector_expression<VecV, Device> const& v, T t){ \
197 typedef typename VecV::value_type type;\
198 typedef typename device_traits<Device>:: template F<type> functor_type;\
199 return vector_binary<VecV, scalar_vector<type, Device>, functor_type >(v(), scalar_vector<type, Device>(v().size(),(type)t), functor_type()); \
200}
201REMORA_VECTOR_SCALAR_TRANSFORMATION(operator/, divide)
202REMORA_VECTOR_SCALAR_TRANSFORMATION(operator<, less)
203REMORA_VECTOR_SCALAR_TRANSFORMATION(operator<=, less_equal)
204REMORA_VECTOR_SCALAR_TRANSFORMATION(operator>, greater)
205REMORA_VECTOR_SCALAR_TRANSFORMATION(operator>=, greater_equal)
206REMORA_VECTOR_SCALAR_TRANSFORMATION(operator==, equal)
207REMORA_VECTOR_SCALAR_TRANSFORMATION(operator!=, not_equal)
208REMORA_VECTOR_SCALAR_TRANSFORMATION(min, min)
209REMORA_VECTOR_SCALAR_TRANSFORMATION(max, max)
210REMORA_VECTOR_SCALAR_TRANSFORMATION(pow, pow)
211#undef REMORA_VECTOR_SCALAR_TRANSFORMATION
212
213// operations of the form op(t,v)[i] = op(t,v[i])
214#define REMORA_VECTOR_SCALAR_TRANSFORMATION_2(name, F)\
215template<class T, class VecV, class Device> \
216typename std::enable_if< \
217 std::is_convertible<T, typename VecV::value_type >::value,\
218 vector_binary<scalar_vector<typename VecV::value_type, Device>, VecV, typename device_traits<Device>:: template F<typename VecV::value_type> > \
219>::type \
220name (T t, vector_expression<VecV, Device> const& v){ \
221 typedef typename VecV::value_type type;\
222 typedef typename device_traits<Device>:: template F<type> functor_type;\
223 return vector_binary<scalar_vector<T, Device>, VecV, functor_type >(scalar_vector<T, Device>(v().size(),t), v() ,functor_type()); \
224}
225REMORA_VECTOR_SCALAR_TRANSFORMATION_2(min, min)
226REMORA_VECTOR_SCALAR_TRANSFORMATION_2(max, max)
227#undef REMORA_VECTOR_SCALAR_TRANSFORMATION_2
228
229template<class VecV1, class VecV2, class Device>
230vector_binary<VecV1, VecV2,
231 typename device_traits<Device>:: template safe_divide<typename common_value_type<VecV1,VecV2>::type >
232>
233safe_div(
234 vector_expression<VecV1, Device> const& v1,
235 vector_expression<VecV2, Device> const& v2,
236 typename common_value_type<VecV1,VecV2>::type defaultValue
237){
238 REMORA_SIZE_CHECK(v1().size() == v2().size());
239 typedef typename common_value_type<VecV1,VecV2>::type result_type;
240
241 typedef typename device_traits<Device>:: template safe_divide<result_type> functor_type;
242 return vector_binary<VecV1, VecV2, functor_type>(v1(),v2(), functor_type(defaultValue));
243}
244
245/////VECTOR REDUCTIONS
246
247/// \brief sum v = sum_i v_i
248template<class VecV, class Device>
249typename VecV::value_type
250sum(vector_expression<VecV, Device> const& v) {
251 typedef typename VecV::value_type value_type;
252 typedef typename device_traits<Device>:: template add<value_type> functor_type;
253 auto const& elem_result = eval_block(v);
254 value_type result = 0;
255 kernels::vector_fold<functor_type>(elem_result,result);
256 return result;
257}
258
259/// \brief max v = max_i v_i
260template<class VecV, class Device>
261typename VecV::value_type
262max(vector_expression<VecV, Device> const& v) {
263 typedef typename VecV::value_type value_type;
264 typedef typename device_traits<Device>:: template max<value_type> functor_type;
265 auto const& elem_result = eval_block(v);
266 value_type result = std::numeric_limits<value_type>::lowest();
267 kernels::vector_fold<functor_type>(elem_result,result);
268 return result;
269}
270
271/// \brief min v = min_i v_i
272template<class VecV, class Device>
273typename VecV::value_type
274min(vector_expression<VecV, Device> const& v) {
275 typedef typename VecV::value_type value_type;
276 typedef typename device_traits<Device>:: template min<value_type> functor_type;
277 auto const& elem_result = eval_block(v);
278 value_type result = std::numeric_limits<value_type>::max();
279 kernels::vector_fold<functor_type>(elem_result,result);
280 return result;
281}
282
283/// \brief arg_max v = arg max_i v_i
284template<class VecV, class Device>
285std::size_t arg_max(vector_expression<VecV, Device> const& v) {
286 REMORA_SIZE_CHECK(v().size() > 0);
287 auto const& elem_result = eval_block(v);
288 return kernels::vector_max(elem_result);
289}
290
291/// \brief arg_min v = arg min_i v_i
292template<class VecV, class Device>
293std::size_t arg_min(vector_expression<VecV, Device> const& v) {
294 REMORA_SIZE_CHECK(v().size() > 0);
295 return arg_max(-v);
296}
297
298/// \brief soft_max v = ln(sum(exp(v)))
299///
300/// Be aware that this is NOT the same function as used in machine learning: exp(v)/sum(exp(v))
301///
302/// The function is computed in an numerically stable way to prevent that too high values of v_i produce inf or nan.
303/// The name of the function comes from the fact that it behaves like a continuous version of max in the respect that soft_max v <= v.size()*max(v)
304/// max is reached in the limit as the gap between the biggest value and the rest grows to infinity.
305template<class VecV, class Device>
306typename VecV::value_type
307soft_max(vector_expression<VecV, Device> const& v) {
308 typename VecV::value_type maximum = max(v);
309 using std::log;
310 return log(sum(exp(v - maximum))) + maximum;
311}
312
313
314////implement all the norms based on sum!
315
316/// \brief norm_1 v = sum_i |v_i|
317template<class VecV, class Device>
318typename real_traits<typename VecV::value_type >::type
319norm_1(vector_expression<VecV, Device> const& v) {
320 return sum(abs(eval_block(v)));
321}
322
323/// \brief norm_2 v = sum_i |v_i|^2
324template<class VecV, class Device>
325typename real_traits<typename VecV::value_type >::type
326norm_sqr(vector_expression<VecV, Device> const& v) {
327 return sum(sqr(eval_block(v)));
328}
329
330/// \brief norm_2 v = sqrt (sum_i |v_i|^2 )
331template<class VecV, class Device>
332typename real_traits<typename VecV::value_type >::type
333norm_2(vector_expression<VecV, Device> const& v) {
334 using std::sqrt;
335 return sqrt(norm_sqr(v));
336}
337
338/// \brief norm_inf v = max_i |v_i|
339template<class VecV, class Device>
340typename real_traits<typename VecV::value_type >::type
341norm_inf(vector_expression<VecV, Device> const& v){
342 return max(abs(eval_block(v)));
343}
344
345/// \brief index_norm_inf v = arg max_i |v_i|
346template<class VecV, class Device>
347std::size_t index_norm_inf(vector_expression<VecV, Device> const& v){
348 return arg_max(abs(eval_block(v)));
349}
350
351// inner_prod (v1, v2) = sum_i v1_i * v2_i
352template<class VecV1, class VecV2, class Device>
353decltype(
354 typename VecV1::value_type() * typename VecV2::value_type()
355)
356inner_prod(
357 vector_expression<VecV1, Device> const& v1,
358 vector_expression<VecV2, Device> const& v2
359) {
360 REMORA_SIZE_CHECK(v1().size() == v2().size());
361 typedef decltype(
362 typename VecV1::value_type() * typename VecV2::value_type()
363 ) value_type;
364 value_type result = value_type();
365 kernels::dot(eval_block(v1),eval_block(v2),result);
366 return result;
367}
368
369
370// Vector Concatenation
371
372
373///\brief Concatenates two vectors
374///
375/// Given two vectors v and w, forms the vector (v,w).
376template<class VecV1, class VecV2, class Device>
377vector_concat<VecV1,VecV2> operator|(
378 vector_expression<VecV1, Device> const& v1,
379 vector_expression<VecV2, Device> const& v2
380){
381 return vector_concat<VecV1,VecV2>(v1(),v2());
382}
383
384///\brief Concatenates a vector with a scalar
385///
386/// Given a vector v and a scalar t, forms the vector (v,t)
387template<class VecV, class T, class Device>
388typename std::enable_if<
389 std::is_convertible<T, typename VecV::value_type>::value,
390 vector_concat<VecV, scalar_vector<T, Device> >
391>::type operator| (
392 vector_expression<VecV, Device> const& v,
393 T t
394){
395 return v | scalar_vector<T, Device>(1,t);
396}
397
398///\brief Concatenates a vector with a scalar
399///
400/// Given a vector v and a scalar t, forms the vector (v,t)
401template<class T, class VecV, class Device>
402typename std::enable_if<
403 std::is_convertible<T, typename VecV::value_type>::value,
404 vector_concat<scalar_vector<T, Device>,VecV >
405>::type operator| (
406 T t,
407 vector_expression<VecV, Device> const& v
408){
409 return scalar_vector<T, Device>(1,t) | v;
410}
411
412
413}
414
415#endif