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 SHARK_LINALG_BLAS_VECTOR_EXPRESSION_HPP
29 #define SHARK_LINALG_BLAS_VECTOR_EXPRESSION_HPP
30 
31 #include <boost/type_traits/is_convertible.hpp>
32 #include <boost/utility/enable_if.hpp>
33 #include "vector_proxy.hpp"
34 #include "kernels/dot.hpp"
35 
36 namespace shark {
37 namespace blas {
38 
39 ///\brief Implements multiplications of a vector by a scalar
40 template<class E>
42  public vector_expression<vector_scalar_multiply <E> > {
44 public:
45  typedef typename E::const_closure_type expression_closure_type;
46  typedef typename E::size_type size_type;
47  typedef typename E::difference_type difference_type;
48  typedef typename E::value_type value_type;
49  typedef typename E::scalar_type scalar_type;
50  typedef value_type const_reference;
51  typedef value_type reference;
52  typedef value_type const * const_pointer;
53  typedef value_type const* pointer;
54 
55  typedef typename E::index_type index_type;
56  typedef typename E::const_index_pointer const_index_pointer;
58 
59  typedef self_type const const_closure_type;
60  typedef self_type closure_type;
61  typedef unknown_storage_tag storage_category;
62  typedef typename E::evaluation_category evaluation_category;
63 
64  // Construction and destruction
65  // May be used as mutable expression.
66  vector_scalar_multiply(vector_expression<E> const &e, scalar_type scalar):
67  m_expression(e()), m_scalar(scalar) {}
68 
69  // Accessors
70  size_type size() const {
71  return m_expression.size();
72  }
73 
74  // Expression accessors
75  expression_closure_type const &expression() const {
76  return m_expression;
77  }
78 
79 public:
80  // Element access
81  const_reference operator()(index_type i) const {
82  return m_scalar * m_expression(i);
83  }
84 
85  const_reference operator[](index_type i) const {
86  return m_scalar * m_expression(i);
87  }
88 
89 
90  //computation kernels
91  template<class VecX>
92  void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
93  m_expression.assign_to(x,alpha*m_scalar);
94  }
95  template<class VecX>
96  void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
97  m_expression.plus_assign_to(x,alpha*m_scalar);
98  }
99 
100  template<class VecX>
101  void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
102  m_expression.minus_assign_to(x,alpha*m_scalar);
103  }
104 
105 
106  //iterators
107  typedef transform_iterator<typename E::const_iterator,scalar_multiply1<value_type, scalar_type> > const_iterator;
108  typedef const_iterator iterator;
109 
110  const_iterator begin() const {
111  return const_iterator(m_expression.begin(),scalar_multiply1<value_type, scalar_type>(m_scalar));
112  }
113  const_iterator end() const {
114  return const_iterator(m_expression.end(),scalar_multiply1<value_type, scalar_type>(m_scalar));
115  }
116 private:
117  expression_closure_type m_expression;
118  scalar_type m_scalar;
119 };
120 
121 
122 template<class T, class E>
123 typename boost::enable_if<
124  boost::is_convertible<T, typename E::scalar_type >,
126 >::type
127 operator* (vector_expression<E> const& e, T scalar){
128  typedef typename E::scalar_type scalar_type;
129  return vector_scalar_multiply<E>(e, scalar_type(scalar));
130 }
131 template<class T, class E>
132 typename boost::enable_if<
133  boost::is_convertible<T, typename E::scalar_type >,
135 >::type
136 operator* (T scalar, vector_expression<E> const& e){
137  typedef typename E::scalar_type scalar_type;
138  return vector_scalar_multiply<E>(e, scalar_type(scalar));//explicit cast prevents warning, alternative would be to template vector_scalar_multiply on T as well
139 }
140 
141 template<class E>
143  typedef typename E::scalar_type scalar_type;
144  return vector_scalar_multiply<E>(e, scalar_type(-1));//explicit cast prevents warning, alternative would be to template vector_scalar_multiply on T as well
145 }
146 
147 /// \brief Vector expression representing a constant valued vector.
148 template<class T>
149 class scalar_vector:public vector_expression<scalar_vector<T> > {
150 
151  typedef scalar_vector<T> self_type;
152 public:
153  typedef std::size_t size_type;
154  typedef std::ptrdiff_t difference_type;
155  typedef T value_type;
156  typedef T scalar_type;
157  typedef const T& const_reference;
158  typedef const_reference reference;
159  typedef T const* const_pointer;
160  typedef const_pointer pointer;
161 
162 
163  typedef std::size_t index_type;
164  typedef index_type const* const_index_pointer;
165  typedef index_type index_pointer;
166 
167  typedef self_type const_closure_type;
168  typedef self_type closure_type;
169  typedef unknown_storage_tag storage_category;
170  typedef elementwise_tag evaluation_category;
171 
172  // Construction and destruction
174  :m_size(0), m_value() {}
175  explicit scalar_vector(size_type size, value_type value)
176  :m_size(size), m_value(value) {}
178  :m_size(v.m_size), m_value(v.m_value) {}
179 
180  // Accessors
181  size_type size() const {
182  return m_size;
183  }
184 
185  // Element access
186  const_reference operator()(index_type /*i*/) const {
187  return m_value;
188  }
189 
190  const_reference operator [](index_type /*i*/) const {
191  return m_value;
192  }
193 
194 public:
195  typedef constant_iterator<T> iterator;
196  typedef constant_iterator<T> const_iterator;
197 
198  const_iterator begin() const {
199  return const_iterator(0,m_value);
200  }
201  const_iterator end() const {
202  return const_iterator(m_size,m_value);
203  }
204 
205 private:
206  size_type m_size;
207  value_type m_value;
208 };
209 
210 ///\brief Creates a vector having a constant value.
211 ///
212 ///@param scalar the value which is repeated
213 ///@param elements the size of the resulting vector
214 template<class T>
215 typename boost::enable_if<boost::is_arithmetic<T>, scalar_vector<T> >::type
216 repeat(T scalar, std::size_t elements){
217  return scalar_vector<T>(elements,scalar);
218 }
219 
220 
221 ///\brief Class implementing vector transformation expressions.
222 ///
223 ///transforms a vector Expression e of type E using a Function f of type F as an elementwise transformation f(e(i))
224 ///This transformation needs f to be constant, meaning that applying f(x), f(y), f(z) yields the same results independent of the
225 ///order of application. Also F must provide a type F::result_type indicating the result type of the functor.
226 template<class E, class F>
228  public vector_expression<vector_unary<E, F> > {
229  typedef vector_unary<E, F> self_type;
230 public:
231  typedef F functor_type;
232  typedef typename E::const_closure_type expression_closure_type;
233  typedef typename E::size_type size_type;
234  typedef typename E::difference_type difference_type;
235  typedef typename F::result_type value_type;
236  typedef value_type scalar_type;
237  typedef value_type const_reference;
238  typedef value_type reference;
239  typedef value_type const * const_pointer;
240  typedef value_type const* pointer;
241 
242  typedef typename E::index_type index_type;
243  typedef typename E::const_index_pointer const_index_pointer;
245 
246  typedef self_type const const_closure_type;
247  typedef self_type closure_type;
248  typedef unknown_storage_tag storage_category;
249  typedef typename E::evaluation_category evaluation_category;
250 
251  // Construction and destruction
252  // May be used as mutable expression.
253  vector_unary(vector_expression<E> const &e, F const &functor):
254  m_expression(e()), m_functor(functor) {}
255 
256  // Accessors
257  size_type size() const {
258  return m_expression.size();
259  }
260 
261  // Expression accessors
262  expression_closure_type const &expression() const {
263  return m_expression;
264  }
265 
266  //computation kernels
267  template<class VecX>
268  void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
269  //compute this by first assigning the result of the argument and then applying
270  //the function to every element
271  assign(x,m_expression);
272  typename VecX::iterator end=x().end();
273  for(typename VecX::iterator pos =x().begin(); pos != end; ++pos){
274  *pos= alpha * m_functor(*pos);
275  }
276  }
277  template<class VecX>
278  void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
279  //First assign result of this expression to a temporary and then perform plus_assignment to x
280  typename vector_temporary<self_type>::type temporary(size());
281  assign_to(temporary,alpha);
282  plus_assign_to(x,temporary);
283  }
284 
285  template<class VecX>
286  void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
287  //First assign result of this expression to a temporary and then perform minus_assignment to x
288  typename vector_temporary<self_type>::type temporary(size());
289  assign_to(temporary,alpha);
290  minus_assign_to(x,temporary);
291  }
292 
293 public:
294  // Element access
295  const_reference operator()(index_type i) const {
296  return m_functor(m_expression(i));
297  }
298 
299  const_reference operator[](index_type i) const {
300  return m_functor(m_expression[i]);
301  }
302 
303  typedef transform_iterator<typename E::const_iterator,functor_type> const_iterator;
304  typedef const_iterator iterator;
305 
306  // Element lookup
307  const_iterator begin() const {
308  return const_iterator(m_expression.begin(),m_functor);
309  }
310  const_iterator end() const {
311  return const_iterator(m_expression.end(),m_functor);
312  }
313 private:
314  expression_closure_type m_expression;
315  F m_functor;
316 };
317 
318 
319 
320 #define SHARK_UNARY_VECTOR_TRANSFORMATION(name, F)\
321 template<class E>\
322 vector_unary<E,F<typename E::value_type> >\
323 name(vector_expression<E> const& e){\
324  typedef F<typename E::value_type> functor_type;\
325  return vector_unary<E, functor_type>(e, functor_type());\
326 }
327 SHARK_UNARY_VECTOR_TRANSFORMATION(abs, scalar_abs)
328 SHARK_UNARY_VECTOR_TRANSFORMATION(log, scalar_log)
329 SHARK_UNARY_VECTOR_TRANSFORMATION(exp, scalar_exp)
330 SHARK_UNARY_VECTOR_TRANSFORMATION(cos, scalar_cos)
331 SHARK_UNARY_VECTOR_TRANSFORMATION(sin, scalar_sin)
332 SHARK_UNARY_VECTOR_TRANSFORMATION(tanh,scalar_tanh)
333 SHARK_UNARY_VECTOR_TRANSFORMATION(atanh,scalar_atanh)
335 SHARK_UNARY_VECTOR_TRANSFORMATION(abs_sqr, scalar_abs_sqr)
336 SHARK_UNARY_VECTOR_TRANSFORMATION(sqrt, scalar_sqrt)
339 SHARK_UNARY_VECTOR_TRANSFORMATION(elem_inv, scalar_inverse)
340 #undef SHARK_UNARY_VECTOR_TRANSFORMATION
341 
342 
343 //operations of the form op(v,t)[i] = op(v[i],t)
344 #define SHARK_VECTOR_SCALAR_TRANSFORMATION(name, F)\
345 template<class T, class E> \
346 typename boost::enable_if< \
347  boost::is_convertible<T, typename E::value_type >,\
348  vector_unary<E,F<typename E::value_type,T> > \
349 >::type \
350 name (vector_expression<E> const& e, T scalar){ \
351  typedef F<typename E::value_type,T> functor_type; \
352  return vector_unary<E, functor_type>(e, functor_type(scalar)); \
353 }
354 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator/, scalar_divide)
355 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator<, scalar_less_than)
356 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator<=, scalar_less_equal_than)
357 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator>, scalar_bigger_than)
358 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator>=, scalar_bigger_equal_than)
359 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator==, scalar_equal)
360 SHARK_VECTOR_SCALAR_TRANSFORMATION(operator!=, scalar_not_equal)
363 SHARK_VECTOR_SCALAR_TRANSFORMATION(pow, scalar_pow)
364 #undef SHARK_VECTOR_SCALAR_TRANSFORMATION
365 
366 // operations of the form op(t,v)[i] = op(t,v[i])
367 #define SHARK_VECTOR_SCALAR_TRANSFORMATION_2(name, F)\
368 template<class T, class E> \
369 typename boost::enable_if< \
370  boost::is_convertible<T, typename E::value_type >,\
371  vector_unary<E,F<typename E::value_type,T> > \
372 >::type \
373 name (T scalar, vector_expression<E> const& e){ \
374  typedef F<typename E::value_type,T> functor_type; \
375  return vector_unary<E, functor_type>(e, functor_type(scalar)); \
376 }
379 #undef SHARK_VECTOR_SCALAR_TRANSFORMATION_2
380 
381 template<class E1, class E2>
382 class vector_addition: public vector_expression<vector_addition<E1,E2> > {
383 private:
384  typedef scalar_binary_plus<
385  typename E1::value_type,
386  typename E2::value_type
387  > functor_type;
388 public:
389  typedef std::size_t size_type;
390  typedef std::ptrdiff_t difference_type;
391  typedef typename functor_type::result_type value_type;
392  typedef value_type scalar_type;
393  typedef value_type const_reference;
394  typedef value_type reference;
395  typedef value_type const * const_pointer;
396  typedef value_type const* pointer;
397 
398  typedef typename E1::index_type index_type;
399  typedef typename E1::const_index_pointer const_index_pointer;
401 
402  typedef typename E1::const_closure_type expression_closure1_type;
403  typedef typename E2::const_closure_type expression_closure2_type;
404 
405  typedef vector_addition<E1,E2> const_closure_type;
406  typedef const_closure_type closure_type;
407  typedef unknown_storage_tag storage_category;
408  typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
409 
410  // Construction and destruction
411  explicit vector_addition (
412  expression_closure1_type e1,
413  expression_closure2_type e2
414  ):m_expression1(e1),m_expression2(e2){
415  SIZE_CHECK(e1.size() == e2.size());
416  }
417 
418  // Accessors
419  size_type size() const {
420  return m_expression1.size();
421  }
422 
423  // Expression accessors
424  expression_closure1_type const& expression1() const {
425  return m_expression1;
426  }
427  expression_closure2_type const& expression2() const {
428  return m_expression2;
429  }
430 
431  // Element access
432  const_reference operator() (index_type i) const {
433  SIZE_CHECK(i < size());
434  return m_expression1(i) + m_expression2(i);
435  }
436 
437  const_reference operator[] (index_type i) const {
438  SIZE_CHECK(i < size());
439  return m_expression1(i) + m_expression2(i);
440  }
441 
442  //computation kernels
443  template<class VecX>
444  void assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
445  assign(x,alpha*m_expression1);
446  plus_assign(x,alpha*m_expression2);
447  }
448  template<class VecX>
449  void plus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
450  plus_assign(x,alpha*m_expression1);
451  plus_assign(x,alpha*m_expression2);
452  }
453 
454  template<class VecX>
455  void minus_assign_to(vector_expression<VecX>& x, scalar_type alpha = scalar_type(1) )const{
456  minus_assign(x,alpha*m_expression1);
457  minus_assign(x,alpha*m_expression2);
458  }
459 
460  // Iterator types
461  typedef binary_transform_iterator<
462  typename E1::const_iterator,
463  typename E2::const_iterator,
464  functor_type
467 
469  return const_iterator(functor_type(),
470  m_expression1.begin(),m_expression1.end(),
471  m_expression2.begin(),m_expression2.end()
472  );
473  }
474  const_iterator end() const {
475  return const_iterator(functor_type(),
476  m_expression1.end(),m_expression1.end(),
477  m_expression2.end(),m_expression2.end()
478  );
479  }
480 
481 private:
482  expression_closure1_type m_expression1;
483  expression_closure2_type m_expression2;
484 };
485 
486 ///\brief Adds two vectors
487 template<class E1, class E2>
489  vector_expression<E1> const& e1,
490  vector_expression<E2> const& e2
491 ){
492  return vector_addition<E1, E2>(e1(),e2());
493 }
494 ///\brief Subtracts two vectors
495 template<class E1, class E2>
497  vector_expression<E1> const& e1,
498  vector_expression<E2> const& e2
499 ){
501 }
502 
503 ///\brief Adds a vector plus a scalr which is interpreted as a constant vector
504 template<class E, class T>
505 typename boost::enable_if<
506  boost::is_convertible<T, typename E::value_type>,
508 >::type operator+ (
509  vector_expression<E> const& e,
510  T t
511 ){
512  return e + scalar_vector<T>(e().size(),t);
513 }
514 
515 ///\brief Adds a vector plus a scalar which is interpreted as a constant vector
516 template<class T, class E>
517 typename boost::enable_if<
518  boost::is_convertible<T, typename E::value_type>,
520 >::type operator+ (
521  T t,
522  vector_expression<E> const& e
523 ){
524  return e + scalar_vector<T>(e().size(),t);
525 }
526 
527 ///\brief Subtracts a scalar which is interpreted as a constant vector from a vector.
528 template<class E, class T>
529 typename boost::enable_if<
530  boost::is_convertible<T, typename E::value_type> ,
532 >::type operator- (
533  vector_expression<E> const& e,
534  T t
535 ){
536  return e - scalar_vector<T>(e().size(),t);
537 }
538 
539 ///\brief Subtracts a vector from a scalar which is interpreted as a constant vector
540 template<class E, class T>
541 typename boost::enable_if<
542  boost::is_convertible<T, typename E::value_type>,
544 >::type operator- (
545  T t,
546  vector_expression<E> const& e
547 ){
548  return scalar_vector<T>(e().size(),t) - e;
549 }
550 
551 
552 template<class E1, class E2, class F>
554  public vector_expression<vector_binary<E1,E2, F> > {
555  typedef vector_binary<E1,E2, F> self_type;
556  typedef E1 const expression1_type;
557  typedef E2 const expression2_type;
558 public:
559  typedef std::size_t size_type;
560  typedef std::ptrdiff_t difference_type;
561  typedef typename F::result_type value_type;
562  typedef value_type scalar_type;
563  typedef value_type const_reference;
564  typedef value_type reference;
565  typedef value_type const * const_pointer;
566  typedef value_type const* pointer;
567 
568  typedef typename E1::index_type index_type;
569  typedef typename E1::const_index_pointer const_index_pointer;
571 
572  typedef F functor_type;
573  typedef typename E1::const_closure_type expression_closure1_type;
574  typedef typename E2::const_closure_type expression_closure2_type;
575 
576  typedef self_type const const_closure_type;
577  typedef self_type closure_type;
578  typedef unknown_storage_tag storage_category;
579  typedef typename evaluation_restrict_traits<E1,E2>::type evaluation_category;
580 
581  // Construction and destruction
582  explicit vector_binary (
583  expression_closure1_type e1,
584  expression_closure2_type e2,
585  F functor
586  ):m_expression1(e1),m_expression2(e2), m_functor(functor) {
587  SIZE_CHECK(e1.size() == e2.size());
588  }
589 
590  // Accessors
591  size_type size() const {
592  return m_expression1.size ();
593  }
594 
595  // Expression accessors
596  expression_closure1_type const& expression1() const {
597  return m_expression1;
598  }
599  expression_closure2_type const& expression2() const {
600  return m_expression2;
601  }
602 
603  // Element access
604  const_reference operator() (index_type i) const {
605  SIZE_CHECK(i < size());
606  return m_functor(m_expression1(i),m_expression2(i));
607  }
608 
609  const_reference operator[] (index_type i) const {
610  SIZE_CHECK(i < size());
611  return m_functor(m_expression1(i),m_expression2(i));
612  }
613 
614  // Iterator types
615 
616  // Iterator enhances the iterator of the referenced expressions
617  // with the unary functor.
618  typedef binary_transform_iterator<
619  typename E1::const_iterator,
620  typename E2::const_iterator,
621  functor_type
624 
626  return const_iterator (m_functor,
627  m_expression1.begin(),m_expression1.end(),
628  m_expression2.begin(),m_expression2.end()
629  );
630  }
631  const_iterator end() const {
632  return const_iterator (m_functor,
633  m_expression1.end(),m_expression1.end(),
634  m_expression2.end(),m_expression2.end()
635  );
636  }
637 
638 private:
639  expression_closure1_type m_expression1;
640  expression_closure2_type m_expression2;
641  F m_functor;
642 };
643 
644 #define SHARK_BINARY_VECTOR_EXPRESSION(name, F)\
645 template<class E1, class E2>\
646 vector_binary<E1, E2, F<typename E1::value_type, typename E2::value_type> >\
647 name(vector_expression<E1> const& e1, vector_expression<E2> const& e2){\
648  typedef F<typename E1::value_type, typename E2::value_type> functor_type;\
649  return vector_binary<E1, E2, functor_type>(e1(),e2(), functor_type());\
650 }
651 SHARK_BINARY_VECTOR_EXPRESSION(operator*, scalar_binary_multiply)
652 SHARK_BINARY_VECTOR_EXPRESSION(element_prod, scalar_binary_multiply)
653 SHARK_BINARY_VECTOR_EXPRESSION(operator/, scalar_binary_divide)
654 SHARK_BINARY_VECTOR_EXPRESSION(element_div, scalar_binary_divide)
655 SHARK_BINARY_VECTOR_EXPRESSION(min, scalar_binary_min)
656 SHARK_BINARY_VECTOR_EXPRESSION(max, scalar_binary_max)
657 #undef SHARK_BINARY_VECTOR_EXPRESSION
658 
659 template<class E1, class E2>
660 vector_binary<E1, E2,
661  scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type>
662 >
664  vector_expression<E1> const& e1,
665  vector_expression<E2> const& e2,
666  typename promote_traits<
667  typename E1::value_type,
668  typename E2::value_type
669  >::promote_type defaultValue
670 ){
671  typedef scalar_binary_safe_divide<typename E1::value_type, typename E2::value_type> functor_type;
672  return vector_binary<E1, E2, functor_type>(e1(),e2(), functor_type(defaultValue));
673 }
674 
675 /////VECTOR REDUCTIONS
676 
677 /// \brief sum v = sum_i v_i
678 template<class E>
679 typename E::value_type
681  typedef typename E::value_type value_type;
682  vector_fold<scalar_binary_plus<value_type, value_type> > kernel;
683  return kernel(e,value_type());
684 }
685 
686 /// \brief max v = max_i v_i
687 template<class E>
688 typename E::value_type
690  typedef typename E::value_type value_type;
691  vector_fold<scalar_binary_max<value_type, value_type> > kernel;
692  return kernel(e,e()(0));
693 }
694 
695 /// \brief min v = min_i v_i
696 template<class E>
697 typename E::value_type
699  typedef typename E::value_type value_type;
700  vector_fold<scalar_binary_min<value_type, value_type> > kernel;
701  return kernel(e,e()(0));
702 }
703 
704 /// \brief arg_max v = arg max_i v_i
705 template<class E>
706 std::size_t arg_max(const vector_expression<E> &e) {
707  SIZE_CHECK(e().size() > 0);
708  return std::max_element(e().begin(),e().end()).index();
709 }
710 
711 /// \brief arg_min v = arg min_i v_i
712 template<class E>
713 std::size_t arg_min(const vector_expression<E> &e) {
714  SIZE_CHECK(e().size() > 0);
715  return arg_max(-e);
716 }
717 
718 /// \brief soft_max v = ln(sum(exp(v)))
719 ///
720 /// Be aware that this is NOT the same function as used in machine learning: exp(v)/sum(exp(v))
721 ///
722 /// The function is computed in an numerically stable way to prevent that too high values of v_i produce inf or nan.
723 /// 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)
724 /// max is reached in the limit as the gap between the biggest value and the rest grows to infinity.
725 template<class E>
726 typename E::value_type
728  typename E::value_type maximum = max(e);
729  return std::log(sum(exp(e - maximum))) + maximum;
730 }
731 
732 
733 ////implement all the norms based on sum!
734 
735 /// \brief norm_1 v = sum_i |v_i|
736 template<class E>
737 typename real_traits<typename E::value_type >::type
739  return sum(abs(eval_block(e)));
740 }
741 
742 /// \brief norm_2 v = sum_i |v_i|^2
743 template<class E>
744 typename real_traits<typename E::value_type >::type
746  return sum(abs_sqr(eval_block(e)));
747 }
748 
749 /// \brief norm_2 v = sqrt (sum_i |v_i|^2 )
750 template<class E>
751 typename real_traits<typename E::value_type >::type
753  using std::sqrt;
754  return sqrt(norm_sqr(e));
755 }
756 
757 /// \brief norm_inf v = max_i |v_i|
758 template<class E>
759 typename real_traits<typename E::value_type >::type
761  return max(abs(eval_block(e)));
762 }
763 
764 /// \brief index_norm_inf v = arg max_i |v_i|
765 template<class E>
766 std::size_t index_norm_inf(vector_expression<E> const &e){
767  return arg_max(abs(eval_block(e)));
768 }
769 
770 // inner_prod (v1, v2) = sum_i v1_i * v2_i
771 template<class E1, class E2>
772 typename promote_traits<
773  typename E1::value_type,
774  typename E2::value_type
775 >::promote_type
777  vector_expression<E1> const& e1,
778  vector_expression<E2> const& e2
779 ) {
780  typedef typename promote_traits<
781  typename E1::value_type,
782  typename E2::value_type
783  >::promote_type value_type;
784  value_type result = value_type();
785  kernels::dot(eval_block(e1),eval_block(e2),result);
786  return result;
787 }
788 
789 }
790 
791 }
792 
793 #endif