MetaTools.h 15.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/** \file MetaTools.h */

#ifndef META_TOOLS_H
#define META_TOOLS_H

// #include "AMDiS.h"
#include <boost/mpl/if.hpp>
#include <boost/mpl/arithmetic.hpp>
#include <boost/mpl/comparison.hpp>
#include <boost/mpl/logical.hpp>
#include <boost/utility/enable_if.hpp>

#define GET_LOOP_INDEX(VAR, SEQ, IDX) static const long VAR = boost::mpl::at_c<SEQ,IDX>::type::value;

// namespace boost {
//   namespace fusion {
//     template<int I> inline double at_c(WorldVector<double> &x) { return x[I]; }
//   }
// }

namespace tools {
  
Praetorius, Simon's avatar
Praetorius, Simon committed
23
24
25
26
  /**
    * \defgroup binomial_coefficients binomial coefficients
    * @{
    */
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
  namespace details {

    template <int N, int K, bool N_pos = true, bool K_pos = true>
    struct Binomial
    {
      BOOST_STATIC_CONSTANT(int, value = (details::Binomial<N-1, K-1, (N-1>0), (K-1>0) >::value
					+ details::Binomial<N-1, K,   (N-1>0), (K>0)   >::value) );
    };
    
    template <int N, bool N_pos>
    struct Binomial<N, N, N_pos, true> // K > 0, N > 0, K==N
    {
      BOOST_STATIC_CONSTANT(int, value = 1);
    };
    
    template <int N, int K, bool N_pos>
    struct Binomial<N, K, N_pos, false> // K <= 0
    {
      BOOST_STATIC_CONSTANT(int, value = 1);
    };
    
    template <int N, int K>
    struct Binomial<N, K, false, true> // K > 0, N <= 0
    {
      BOOST_STATIC_CONSTANT(int, value = 0);
    };
  }
  
Praetorius, Simon's avatar
Praetorius, Simon committed
55
56
57
  /// \brief
  /// interface to calc binomial coefficients (N, K) at compile-time for N,K >= 0
  ///
58
59
60
61
62
  template <int N, int K>
  struct Binomial
  {
    BOOST_STATIC_CONSTANT(int, value = (details::Binomial<N, K, (N>0), (K>0)>::value) );
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
63
  /**@}*/
64

Praetorius, Simon's avatar
Praetorius, Simon committed
65
66
67
68
69
70
// ______________________________________________________________________________________________________
  
  /**
    * \defgroup math_basics some mathematical basics
    * @{
    */
71
  
Praetorius, Simon's avatar
Praetorius, Simon committed
72
73
74
  /// \brief
  /// check divisibility
  ///
75
76
77
78
79
80
  template <int u, int v>
  struct Divisible
  {
    BOOST_STATIC_CONSTANT(bool, value = u % v == 0 ? true : false );
  };
    
Praetorius, Simon's avatar
Praetorius, Simon committed
81
82
83
  /// \brief
  /// check divide by zero condition
  ///
84
85
86
87
88
89
  template <int u>
  struct Divisible<u, 0>
  {
    BOOST_STATIC_CONSTANT(bool, value = false ); // better: static assert
  };
  
Praetorius, Simon's avatar
Praetorius, Simon committed
90
91
92
  /// \brief
  /// check number is even or not
  ///
93
94
95
96
97
98
  template <int u>
  struct IsEven
  {
    BOOST_STATIC_CONSTANT(bool, value = (Divisible<u, 2>::value) );
  };
    
Praetorius, Simon's avatar
Praetorius, Simon committed
99
100
101
  /// \brief
  /// check number is odd or not
  ///
102
103
104
105
106
107
  template <int u>
  struct IsOdd
  {
    BOOST_STATIC_CONSTANT(bool, value = (-Divisible<u, 2>::value) );
  };
  
Praetorius, Simon's avatar
Praetorius, Simon committed
108
109
110
  /// \brief
  /// absolute value calculation
  ///
111
112
113
114
115
116
  template <int N>
  struct Abs
  {
    BOOST_STATIC_CONSTANT(bool, value = (N >= 0) ? N : -N );
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
117
118
119
120
121
122
// ______________________________________________________________________________________________________
//  greatest common devisor and least common multiple

  /// \brief
  /// calculate gcd of template parameters u and v
  ///
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
  template <int u, int v>
  struct Gcd
  {
    BOOST_STATIC_CONSTANT(int, value = (Gcd<v, u % v>::value) );
  };
    
  template <int u>
  struct Gcd<u, 0>
  {
    BOOST_STATIC_CONSTANT(int, value = u );
  };
    
  template <>
  struct Gcd<0, 0>
  {
    BOOST_STATIC_CONSTANT(int, value = -1 );
  };
  
Praetorius, Simon's avatar
Praetorius, Simon committed
141
142
143
  /// \brief
  /// calculate the least common multiple of template parameters u and v
  ///
144
145
146
147
148
  template <int u, int v>
  struct Lcm
  {
    BOOST_STATIC_CONSTANT(int, value = (Abs<u * v>::value / Gcd<u, v>::value) );
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
149
150
151
152
153
154
155

// ______________________________________________________________________________________________________
// powers of values
 
  /// \brief
  /// calculate the power of template parameters a<sup>p</sup>
  ///
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  template <int a, int p>
  struct Pow
  {
    BOOST_STATIC_CONSTANT(int, value = (a * Pow<a, p-1>::value) );
  };
  
  template <int a>
  struct Pow<a, 1>
  {
    BOOST_STATIC_CONSTANT(int, value = a );
  };
  
  template <int a>
  struct Pow<a, 0>
  {
    BOOST_STATIC_CONSTANT(int, value = 1 );
  };
  
  template <int n>
  struct Sqr
  {
    BOOST_STATIC_CONSTANT(int, value = (Pow<n, 2>::value) );
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
179
180
181
182
183
184
185

// ______________________________________________________________________________________________________
// roots of values
 
  /// \brief
  /// template to compute sqrt(a) via iteration
  ///
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
  template <int a, int I=1>
  struct Sqrt {
      // instantiate next step or result type as branch
      typedef typename boost::mpl::if_c< (I*I < a),
					Sqrt<a, I+1>,
					boost::mpl::int_<I>
				      >::type impl;

      // use the result of branch type
      BOOST_STATIC_CONSTANT(int, value = impl::value );
  };
  
  // partial specialization to end the iteration
  template<int a>
  struct Sqrt<a, a> {
    BOOST_STATIC_CONSTANT(int, value = a );
  };
  template<int I>
  struct Sqrt<0, I> {
    BOOST_STATIC_CONSTANT(int, value = 0 );
  };
  
Praetorius, Simon's avatar
Praetorius, Simon committed
208
209
210
  /// \brief
  /// template to compute the Nth root of a via iteration
  ///
211
#ifndef _WIN32
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
  template <int a, int N, int I=1>
  struct Root {
    // instantiate next step or result type as branch
    typedef typename boost::mpl::if_c<Pow<I, N>::value < a,
				      Root<a, N, I+1>,
				      boost::mpl::int_<I>
				    >::type impl;

    // use the result of branch type
    BOOST_STATIC_CONSTANT(int, value = impl::value );
  };
  
  // partial specialization to end the iteration
  template<int a, int N>
  struct Root<a, N, a> {
    BOOST_STATIC_CONSTANT(int, value = a );
  };
  template<int N, int I>
  struct Root<0, N, I> {
    BOOST_STATIC_CONSTANT(int, value = 0 );
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
233
  /**@}*/
234
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
235
236
// ______________________________________________________________________________________________________
// monomials and polynomials
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287

  template<int N>
  struct Monomial {
    typedef double result_type;
    result_type operator()(double x) { return x*Monomial<N-1>()(x); }
  };
  
  template<>
  struct Monomial<1> {
    typedef double result_type;
    result_type operator()(double x) { return x; }
  };
  
  template<>
  struct Monomial<0> {
    typedef double result_type;
    result_type operator()(double x) { return 1.0; }
  };

  template<typename Coefficients, int I, int N>
  struct Horner {
    typedef double result_type;
    Horner(Coefficients& coefficients_) : coefficients(coefficients_) {}
    result_type operator()(double x) {
      return boost::fusion::at_c<I>(coefficients) + x*(Horner<Coefficients, I+1, N>(coefficients)(x));
    }
  private:
    Coefficients& coefficients;
  };

  template<typename Coefficients, int N>
  struct Horner<Coefficients, N, N> {
    typedef double result_type;
    Horner(Coefficients& coefficients_) : coefficients(coefficients_) {}
    result_type operator()(double x) {
      return boost::fusion::at_c<N>(coefficients);
    }
  private:
    Coefficients& coefficients;
  };

  template<typename Coefficients>
  struct Polynomial {
    typedef double result_type;
    Polynomial(Coefficients& coefficients_) : coefficients(coefficients_) {}
    result_type operator()(double x) {
      return Horner<Coefficients, 0, boost::fusion::result_of::size<Coefficients>::value - 1>(coefficients)(x);
    }
  private:
    Coefficients& coefficients;
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
288
289
290
291

// ______________________________________________________________________________________________________
// Bernstein polynomials and Bezier-Surfaces

292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
  template<int I, int order>
  struct BernsteinPolynomial 
  {
    static inline double value(double x) {
      return Binomial<order, I>::value
	      * std::pow(x, I) * std::pow(1.0 - x, order - I);
    }
  };
  
  template<int I, int order, bool valid = true>
  struct BernsteinPolynomial_
  {
    static inline double value(double x) {
      return (1.0 - x) * BernsteinPolynomial_<I, order - 1, I >= 0 && I <= (order-1)>::value(x)
	+ x * BernsteinPolynomial_<I - 1, order - 1, (I-1) >= 0 && (I-1) <= (order-1)>::value(x);
    }
  };
  
  template<>
  struct BernsteinPolynomial_<0, 0, true>
  {
    static inline double value(double x) {
      return 1.0;
    }
  };
  
  template<int I, int order>
  struct BernsteinPolynomial_<I, order, false>
  {
    static inline double value(double x) {
      return 0.0;
    }
  };
  
  template<int I, int J, int order>
  struct BernsteinPolynomial2
  {
    static inline double value(double x, double y) {
      return BernsteinPolynomial<I, order>::value(x) * BernsteinPolynomial<J, order>::value(x);
    }
  };
  
  template<int I, int J, int K, int order>
  struct BernsteinPolynomial3
  {
    static inline double value(double x, double y, double z) {
      return BernsteinPolynomial<I, order>::value(x) * BernsteinPolynomial<J, order>::value(x) * BernsteinPolynomial<K, order>::value(z);
    }
  };
Praetorius, Simon's avatar
Praetorius, Simon committed
341
342
343

// ______________________________________________________________________________________________________
// Bezier-Polynomial
344
#ifndef _WIN32
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
  namespace details {
    
    // sum_{I=0}^N B_I^N(x)*c_I
    template<int I, int N>
    struct BezierPolynomial
    {
      template<typename Coefficients>
      static double value(Coefficients& coefficients, double x) {
	return boost::fusion::at_c<I>(coefficients)*BernsteinPolynomial<I, N>::value(x) + details::BezierPolynomial<I+1, N>::value(coefficients, x);
      }
    };
    
    template<int N>
    struct BezierPolynomial<N, N>
    {
      template<typename Coefficients>
      static double value(Coefficients& coefficients, double x) {
	return boost::fusion::at_c<N>(coefficients)*BernsteinPolynomial<N,N>::value(x);
      }
    };
  }
  
  template<typename Coefficients, int DOW=1>
  struct BezierPolynomial
  {
    BOOST_STATIC_CONSTANT(int, N = (boost::mpl::minus< boost::mpl::int_< Root< boost::fusion::result_of::size<Coefficients>::value, 
									       DOW >::value >, 
						      boost::mpl::int_<1> >::type::value) ); // N = rootN(coeff.size, DOW)-1
    
    BezierPolynomial(Coefficients& coefficients_) : coefficients(coefficients_) {}
    double value(double x) {
      return details::BezierPolynomial<0, N>::value(coefficients, x);
    }
  private:
    Coefficients& coefficients;
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
382
383
// ______________________________________________________________________________________________________

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
  
  namespace details {
    
    template<int IJ, int N, bool finished=false> // IJ = I + J*N, max(IJ) = N*N - 2
    struct BezierPolynomial2
    {
      template<typename Coefficients>
      static double value(Coefficients& coefficients, double x, double y) {
	return boost::fusion::at_c<N>(coefficients) * BernsteinPolynomial2<IJ%N, IJ/N, N>::value(x, y)
	  + details::BezierPolynomial2<IJ + 1, N, IJ == (N+1)*(N+1)>::value(coefficients, x, y);
      }
    };
    
    template<int IJ, int N>
    struct BezierPolynomial2<IJ, N, true>
    {
      template<typename Coefficients>
      static double value(Coefficients& coefficients, double x, double y) {
	return 0.0;
      }
    };
  }
    
Praetorius, Simon's avatar
Praetorius, Simon committed
407
408
409
  /// \brief
  /// interface for bezier-surfaces of order (N,N)
  ///
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
  template<typename Coefficients>
  struct BezierPolynomial2
  {
    BOOST_STATIC_CONSTANT(int, N = (boost::mpl::minus< 
				      boost::mpl::int_<
					Root<
					  boost::fusion::result_of::size<Coefficients>::value, 
					  2
					>::value
				      >, 
				      boost::mpl::int_<1>
				    >::type::value) ); // N = rootN(coeff.size, 2)-1
    
    BezierPolynomial2(Coefficients& coefficients_) : coefficients(coefficients_) {} 
    double value(double x, double y) {
      return details::BezierPolynomial2<0, N, N == 0>::value(coefficients, x, y);
    }
  private:
    Coefficients& coefficients; // p in points is element of result_type
  };
430
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
431
432
// ______________________________________________________________________________________________________
// for loops: for<i, n [,{incr,decr}]>, for<mpl::range_c>, for<mpl::vector_c>
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465

  template<int I> struct incr { BOOST_STATIC_CONSTANT(int, value = I+1); };
  template<int I> struct decr { BOOST_STATIC_CONSTANT(int, value = I-1); };

  template<int I, int N, template<int> class op, bool finished, template<int> class F>
  struct FOR_LOOP {
    static void loop() {
      F<I>::eval(); // your implementation here
      FOR_LOOP<op<I>::value, N, op, op<I>::value == N, F>::loop();
    }
  };

  // initial/break condition for 'for-loop': e.g. I==N
  template<int I, int N, template<int> class op, template<int> class F> 
  struct FOR_LOOP<I, N, op, true, F> 
  {
    static void loop() {}
  };

  // 'for-loop', for I<N
  template<int I, int N, bool plus, bool finished, template<int> class F> 
  struct FOR_LOOP_CHECK 
  {
    static void loop() { FOR_LOOP<I, N, incr, finished, F>::loop(); }
  };

  // 'for-loop', for I>N
  template<int I, int N, bool finished, template<int> class F> 
  struct FOR_LOOP_CHECK<I,N,false,finished, F> 
  {
    static void loop() { FOR_LOOP<I, N, decr, finished, F>::loop(); }
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
466
467
468
469
  /// \brief
  /// interface for classical for loop: for j=I..N
  /// the index will be incremented, if I&lt;N and decremented otherwise
  ///
470
471
472
473
474
475
  template<int I, int N, template<int> class F> 
  struct FOR 
  { 
    static void loop() { FOR_LOOP_CHECK<I, N, I < N, I == N, F>::loop(); } 
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
476
  // for-loop over sequence of values
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
  template<typename Seq, long I, long N, template<int> class F>
  struct FOR_SUBSEQ
  {
    static void loop() {
      GET_LOOP_INDEX (i, Seq, I);
      F<i>::eval(); // your implementation here
      FOR_SUBSEQ<Seq, I+1, N, F>::loop();
    }
  };

  // initial/break condition for 'for-loop' over sequences: I==N
  template<typename Seq, long N, template<int> class F> 
  struct FOR_SUBSEQ<Seq,N,N,F> 
  { 
    static void loop() {} 
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
494
495
496
497
  /// \brief
  /// interface for loop over sequence of values
  /// for i in {seq[0]..seq[end]}
  ///
498
499
500
501
502
503
  template<typename Seq, template<int> class F> 
  struct FOR_SEQ 
  { 
    static void loop() { FOR_SUBSEQ<Seq, 0, boost::mpl::size<Seq>::value, F>::loop(); } 
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
504
505
//____________________________________________________________________________
// for-each loops
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529

  template<int I> struct for_each
  {
    template<typename Container1, typename Container2, typename Binary_Functor>
    static void eval(Container1& vec1, Container2& vec2, Binary_Functor& f);
    
    template<typename Container1, typename Container2, typename T, typename Tertiary_Functor>
    static void accumulate(Container1& vec1, Container2& vec2, T& value, Tertiary_Functor& f);
  };

  template<> struct for_each<0>
  {
    template<typename Container1, typename Container2, typename Binary_Functor>
    static void eval(Container1& vec1, Container2& vec2, Binary_Functor& f);

    template<typename Container1, typename Container2, typename T, typename Tertiary_Functor>
    static void accumulate(Container1& vec1, Container2& vec2, T& value, Tertiary_Functor& f);
  };

  template<int I>
  template<typename Container1, typename Container2, typename Binary_Functor>
  void for_each<I>::eval(Container1& vec1, Container2& vec2, Binary_Functor& f) {
    f(boost::fusion::at_c<I>(vec1), boost::fusion::at_c<I>(vec2));
    for_each<I-1>::eval(vec1,vec2,f);
530
  }
531
532
533
534
535
536

  template<int I>
  template<typename Container1, typename Container2, typename T, typename Tertiary_Functor>
  void for_each<I>::accumulate(Container1& vec1, Container2& vec2, T& value, Tertiary_Functor& f) {
    f(boost::fusion::at_c<I>(vec1), boost::fusion::at_c<I>(vec2), value);
    for_each<I-1>::accumulate(vec1,vec2,value,f);
537
  }
538
539
540
541

  template<typename Container1, typename Container2, typename Binary_Functor>
  void for_each<0>::eval(Container1& vec1, Container2& vec2, Binary_Functor& f) {
    f(boost::fusion::at_c<0>(vec1), boost::fusion::at_c<0>(vec2));
542
  }
543
544
545
546

  template<typename Container1, typename Container2, typename T, typename Tertiary_Functor>
  void for_each<0>::accumulate(Container1& vec1, Container2& vec2, T& value, Tertiary_Functor& f) {
    f(boost::fusion::at_c<0>(vec1), boost::fusion::at_c<0>(vec2),value);
547
  }
548
549
550
551
  
} // namespace tools

#endif // META_TOOLS_H