1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // MetaPhysicL - A metaprogramming library for physics calculations
5 //
6 // Copyright (C) 2013 The PECOS Development Team
7 //
8 // This library is free software; you can redistribute it and/or
9 // modify it under the terms of the Version 2.1 GNU Lesser General
10 // Public License as published by the Free Software Foundation.
11 //
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // Lesser General Public License for more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public
18 // License along with this library; if not, write to the Free Software
19 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
20 // Boston, MA 02110-1301 USA
21 //
22 //-----------------------------------------------------------------------el-
23 //
24 // $Id: core.h 37197 2013-02-21 05:49:09Z roystgnr $
25 //
26 //--------------------------------------------------------------------------
27
28
29 #ifndef METAPHYSICL_NAMEDINDEXARRAY_H
30 #define METAPHYSICL_NAMEDINDEXARRAY_H
31
32 #include <algorithm>
33 #include <functional>
34 #include <stdexcept>
35 #include <ostream>
36 #include <type_traits>
37
38 #include "metaphysicl/compare_types.h"
39 #include "metaphysicl/ct_set.h"
40 #include "metaphysicl/metaprogramming.h"
41 #include "metaphysicl/metaphysicl_asserts.h"
42 #include "metaphysicl/raw_type.h"
43
44 namespace MetaPhysicL {
45
46 template <typename DataType, typename SizeType, typename PermType>
47 const typename boostcopy::enable_if<
48 BuiltinTraits<DataType>, DataType&>::type
reshape(DataType & data,const SizeType & sizes,const PermType &)49 reshape(DataType& data,
50 const SizeType& sizes,
51 const PermType& /* perm */)
52 {
53 #ifndef NDEBUG
54 for (unsigned int i = 0; i != sizes.size(); ++i)
55 metaphysicl_assert_equal_to(sizes[i], 1);
56 #endif
57 return data;
58 }
59
60 template <typename DataType, typename SizeType, typename PermType>
61 const typename boostcopy::enable_if<
62 BuiltinTraits<DataType>, const DataType&>::type
reshape(const DataType & data,const SizeType & sizes,const PermType &)63 reshape(const DataType& data,
64 const SizeType& sizes,
65 const PermType& /* perm */)
66 {
67 #ifndef NDEBUG
68 for (unsigned int i = 0; i != sizes.size(); ++i)
69 metaphysicl_assert_equal_to(sizes[i], 1);
70 #endif
71 return data;
72 }
73
74 template <typename DataVector, typename SparseSizeVector>
75 class NamedIndexArray
76 {
77 public:
78 typedef typename SparseSizeVector::index_set index_set;
79
80 typedef typename index_set::head_type index_type;
81
82 static const size_t index_size = index_set::size;
83
size()84 std::size_t size() const
85 { return _data_vector.size(); }
86
NamedIndexArray()87 NamedIndexArray() :
88 _data_vector(), _size_vector() {}
89
NamedIndexArray(DataVector vec_in,SparseSizeVector size_in)90 NamedIndexArray(DataVector vec_in, SparseSizeVector size_in) :
91 _data_vector(vec_in), _size_vector(size_in) {}
92
93 template <typename DataVector2, typename SparseSizeVector2>
NamedIndexArray(const NamedIndexArray<DataVector2,SparseSizeVector2> & arr_in)94 NamedIndexArray(const NamedIndexArray<DataVector2,SparseSizeVector2>& arr_in) :
95 _data_vector(arr_in.raw_data()), _size_vector(arr_in.raw_sizes()) {}
96
97 template <typename DataVector2, typename SparseSizeVector2>
98 NamedIndexArray<DataVector, SparseSizeVector>&
99 operator= (const NamedIndexArray<DataVector2,SparseSizeVector2>& arr_in) {
100 _data_vector = arr_in.raw_data();
101 _size_vector = arr_in.raw_sizes();
102 return *this;
103 }
104
raw_data()105 DataVector& raw_data()
106 { return _data_vector; }
107
raw_data()108 const DataVector& raw_data() const
109 { return _data_vector; }
110
raw_sizes()111 SparseSizeVector& raw_sizes()
112 { return _size_vector; }
113
raw_sizes()114 const SparseSizeVector& raw_sizes() const
115 { return _size_vector; }
116
117 auto
118 operator! () const
119 -> NamedIndexArray <decltype(!this->raw_data()), SparseSizeVector>
120 {
121 return NamedIndexArray
122 <decltype(!_data_vector), SparseSizeVector>
123 (!_data_vector, _size_vector);
124 }
125
126 auto
127 operator- () const
128 -> NamedIndexArray <decltype(-this->raw_data()), SparseSizeVector>
129 {
130 return NamedIndexArray
131 <decltype(-_data_vector), SparseSizeVector>
132 (-_data_vector, _size_vector);
133 }
134
135
136 #define NamedIndexArray_opequals(opname) \
137 template <typename DataVector2, typename SparseSizeVector2> \
138 NamedIndexArray<DataVector,SparseSizeVector>& \
139 operator opname (const NamedIndexArray<DataVector2,SparseSizeVector2>& a) { \
140 typedef typename SparseSizeVector::index_set IndexSet; \
141 typedef typename SparseSizeVector2::index_set IndexSet2; \
142 using MetaPhysicL::PermutationArray; \
143 ctassert<IndexSet2::template Difference<IndexSet>::type::size == 0>::apply(); \
144 _data_vector opname reshape(a.raw_data(), \
145 _size_vector.raw_data_array(), \
146 PermutationArray<IndexSet2,IndexSet>::value()); \
147 return *this; \
148 } \
149 \
150 template <typename T2> \
151 NamedIndexArray<DataVector,SparseSizeVector>& \
152 operator opname (const T2& a) \
153 { _data_vector opname a; return *this; }
154
155 NamedIndexArray_opequals(+=)
156 NamedIndexArray_opequals(-=)
157 NamedIndexArray_opequals(*=)
158 NamedIndexArray_opequals(/=)
159
160 private:
161 DataVector _data_vector;
162 SparseSizeVector _size_vector;
163 };
164
165
166 //
167 // Non-member functions
168 //
169
170
171 #define NamedIndexArray_op_ab(opname) \
172 template <typename DataVector, typename DataVector2, \
173 typename SparseSizeVector, typename SparseSizeVector2> \
174 inline \
175 auto \
176 operator opname (const NamedIndexArray<DataVector, SparseSizeVector>& a, \
177 const NamedIndexArray<DataVector2, SparseSizeVector2>& b) \
178 { \
179 typedef typename SparseSizeVector::index_set IndexSet; \
180 typedef typename SparseSizeVector2::index_set IndexSet2; \
181 typedef typename IndexSet2::template Union<IndexSet>::type UnionSet; \
182 using std::max; \
183 using MetaPhysicL::PermutationArray; \
184 const auto final_sizes = max(a.raw_sizes(), b.raw_sizes()); \
185 return NamedIndexArray<decltype( \
186 reshape(a.raw_data(), \
187 final_sizes.raw_data_array(), \
188 PermutationArray<IndexSet,UnionSet>::value()) \
189 opname \
190 reshape(b.raw_data(), \
191 final_sizes.raw_data_array(), \
192 PermutationArray<IndexSet2,UnionSet>::value())), \
193 typename std::remove_const<decltype(final_sizes)>::type> ( \
194 reshape(a.raw_data(), \
195 final_sizes.raw_data_array(), \
196 PermutationArray<IndexSet,UnionSet>::value()) \
197 opname \
198 reshape(b.raw_data(), \
199 final_sizes.raw_data_array(), \
200 PermutationArray<IndexSet2,UnionSet>::value()), \
201 final_sizes); \
202 } \
203 \
204 template <typename DataVector, typename SparseSizeVector, typename T2> \
205 inline \
206 auto \
207 operator opname (const NamedIndexArray<DataVector, SparseSizeVector>& a, \
208 const T2& b) \
209 -> typename enable_if_c<DefinesSupertype<CompareTypes< \
210 NamedIndexArray<DataVector, SparseSizeVector>, \
211 T2 \
212 > \
213 >::value, \
214 NamedIndexArray<decltype(a.raw_data() opname b), SparseSizeVector> \
215 >::type \
216 { \
217 return NamedIndexArray<decltype(a.raw_data() opname b), SparseSizeVector> \
218 (a.raw_data() opname b, a.raw_sizes()); \
219 } \
220 \
221 template <typename DataVector, typename SparseSizeVector, typename T> \
222 inline \
223 auto \
224 operator opname (const T& a, \
225 const NamedIndexArray<DataVector, SparseSizeVector>& b) \
226 -> typename enable_if_c<DefinesSupertype<CompareTypes< \
227 NamedIndexArray<DataVector, SparseSizeVector>, \
228 T \
229 > \
230 >::value, \
231 NamedIndexArray<decltype(a opname b.raw_data()), SparseSizeVector> \
232 >::type \
233 { \
234 return NamedIndexArray<decltype(a opname b.raw_data()), SparseSizeVector> \
235 (a opname b.raw_data(), b.raw_sizes()); \
236 }
237
238 NamedIndexArray_op_ab(+)
239 NamedIndexArray_op_ab(-)
240 NamedIndexArray_op_ab(*)
241 NamedIndexArray_op_ab(/)
242 NamedIndexArray_op_ab(<)
243 NamedIndexArray_op_ab(<=)
244 NamedIndexArray_op_ab(>)
245 NamedIndexArray_op_ab(>=)
246 NamedIndexArray_op_ab(==)
247 NamedIndexArray_op_ab(!=)
248
249
250 template <typename DataVector, typename SparseSizeVector>
251 inline
252 std::ostream&
253 operator<< (std::ostream& output,
254 const NamedIndexArray<DataVector, SparseSizeVector>& a)
255 {
256 output << "Sizes: " << a.raw_sizes() << '\n';
257 output << "Data: " << a.raw_data() << '\n';
258 return output;
259 }
260
261 // ScalarTraits, RawType, CompareTypes specializations
262
263 template <typename DataVector, typename SparseSizeVector>
264 struct ScalarTraits<NamedIndexArray<DataVector, SparseSizeVector> >
265 {
266 static const bool value = ScalarTraits<DataVector>::value;
267 };
268
269 template <typename DataVector, typename SparseSizeVector>
270 struct RawType<NamedIndexArray<DataVector, SparseSizeVector> >
271 {
272 typedef typename RawType<DataVector>::value_type value_type;
273
274 static value_type value(const NamedIndexArray<DataVector, SparseSizeVector>& a)
275 { return raw_value(a.raw_data()); }
276 };
277
278 // NamedIndexArray CompareTypes values aren't accurate, but we can
279 // still test for their existence to determine what types are
280 // comparable
281 template <typename DataVector, typename SparseSizeVector, typename T2, bool reverseorder>
282 struct CompareTypes<NamedIndexArray<DataVector, SparseSizeVector>, T2, reverseorder,
283 typename boostcopy::enable_if<BuiltinTraits<T2> >::type> {
284 typedef bool supertype;
285 };
286
287 template <typename DataVector, typename SparseSizeVector,
288 typename DataVector2, typename SparseSizeVector2, bool reverseorder>
289 struct CompareTypes<NamedIndexArray<DataVector, SparseSizeVector>,
290 NamedIndexArray<DataVector2, SparseSizeVector2>,
291 reverseorder> {
292 typedef bool supertype;
293 };
294
295 template <typename DataVector, typename SparseSizeVector, bool reverseorder>
296 struct CompareTypes<NamedIndexArray<DataVector, SparseSizeVector>,
297 NamedIndexArray<DataVector, SparseSizeVector>,
298 reverseorder> {
299 typedef bool supertype;
300 };
301
302
303
304 template <typename DataVector, typename SparseSizeVector>
305 struct
306 copy_or_reference<NamedIndexArray<DataVector, SparseSizeVector> &>
307 {
308 typedef typename
309 IfElse<copy_or_reference<DataVector&>::copy,
310 NamedIndexArray<DataVector, SparseSizeVector>,
311 NamedIndexArray<DataVector, SparseSizeVector>&>::type type;
312
313 static const bool copy = copy_or_reference<DataVector&>::copy;
314 };
315
316
317 template <typename DataVector, typename SparseSizeVector>
318 struct
319 copy_or_reference<const NamedIndexArray<DataVector, SparseSizeVector> &>
320 {
321 typedef typename
322 IfElse<copy_or_reference<DataVector&>::copy,
323 NamedIndexArray<DataVector, SparseSizeVector>,
324 const NamedIndexArray<DataVector, SparseSizeVector>&>::type type;
325
326 static const bool copy = copy_or_reference<DataVector&>::copy;
327 };
328
329
330
331 } // namespace MetaPhysicL
332
333
334 namespace std {
335
336 using MetaPhysicL::NamedIndexArray;
337 using MetaPhysicL::CompareTypes;
338 using MetaPhysicL::enable_if_c;
339 using MetaPhysicL::DefinesSupertype;
340
341 #define NamedIndexArray_std_unary(funcname) \
342 template <typename DataVector, typename SparseSizeVector> \
343 inline \
344 auto \
345 funcname (const NamedIndexArray<DataVector, SparseSizeVector> &a) \
346 { \
347 return NamedIndexArray<decltype(std::funcname(a.raw_data())), \
348 SparseSizeVector> \
349 (std::funcname(a.raw_data()), a.raw_sizes()); \
350 }
351
352
353 #define NamedIndexArray_fl_unary(funcname) \
354 NamedIndexArray_std_unary(funcname##f) \
355 NamedIndexArray_std_unary(funcname##l)
356
357
358 #define NamedIndexArray_stdfl_unary(funcname) \
359 NamedIndexArray_std_unary(funcname) \
360 NamedIndexArray_fl_unary(funcname)
361
362
363 #define NamedIndexArray_std_binary(funcname) \
364 template <typename DataVector, typename DataVector2, \
365 typename SparseSizeVector, typename SparseSizeVector2> \
366 inline \
367 auto \
368 funcname (const NamedIndexArray<DataVector , SparseSizeVector >& a, \
369 const NamedIndexArray<DataVector2, SparseSizeVector2>& b) \
370 { \
371 typedef typename SparseSizeVector::index_set IndexSet; \
372 typedef typename SparseSizeVector2::index_set IndexSet2; \
373 typedef typename IndexSet2::template Union<IndexSet>::type UnionSet; \
374 using std::max; \
375 const auto final_sizes = max(a.raw_sizes(), b.raw_sizes()); \
376 \
377 using std::funcname; \
378 using MetaPhysicL::PermutationArray; \
379 return NamedIndexArray<decltype( \
380 funcname( \
381 reshape(a.raw_data(), \
382 final_sizes, \
383 PermutationArray<IndexSet,UnionSet>::value()), \
384 reshape(b.raw_data(), \
385 final_sizes, \
386 PermutationArray<IndexSet2,UnionSet>::value()))), \
387 decltype(final_sizes)> (\
388 funcname( \
389 reshape(a.raw_data(), \
390 final_sizes, \
391 PermutationArray<IndexSet,UnionSet>::value()), \
392 reshape(b.raw_data(), \
393 final_sizes, \
394 PermutationArray<IndexSet2,UnionSet>::value())) \
395 , final_sizes); \
396 } \
397 \
398 template <typename DataVector, typename SparseSizeVector, typename T2> \
399 inline \
400 auto \
401 funcname (const NamedIndexArray<DataVector, SparseSizeVector>& a, \
402 const T2& b) \
403 -> typename enable_if_c<DefinesSupertype<CompareTypes< \
404 NamedIndexArray<DataVector, SparseSizeVector>, \
405 T2 \
406 > \
407 >::value, \
408 NamedIndexArray<decltype(funcname(a.raw_data(), b)), SparseSizeVector> \
409 >::type \
410 { \
411 using std::funcname; \
412 return NamedIndexArray<decltype(funcname(a.raw_data(), b)), SparseSizeVector> \
413 (funcname(a.raw_data(), b), a.raw_sizes()); \
414 } \
415 \
416 template <typename DataVector, typename SparseSizeVector, typename T> \
417 inline \
418 auto \
419 funcname (const T& a, \
420 const NamedIndexArray<DataVector, SparseSizeVector>& b) \
421 -> typename enable_if_c<DefinesSupertype<CompareTypes< \
422 NamedIndexArray<DataVector, SparseSizeVector>, \
423 T \
424 > \
425 >::value, \
426 NamedIndexArray<decltype(funcname(a, b.raw_data())), SparseSizeVector> \
427 >::type \
428 { \
429 using std::funcname; \
430 return NamedIndexArray<decltype(funcname(a, b.raw_data())), SparseSizeVector> \
431 (funcname(a, b.raw_data()), b.raw_sizes()); \
432 }
433
434
435 #define NamedIndexArray_fl_binary(funcname) \
436 NamedIndexArray_std_binary(funcname##f) \
437 NamedIndexArray_std_binary(funcname##l)
438
439
440 #define NamedIndexArray_stdfl_binary(funcname) \
441 NamedIndexArray_std_binary(funcname) \
442 NamedIndexArray_fl_binary(funcname)
443
444
445 NamedIndexArray_std_binary(pow)
446 NamedIndexArray_std_unary(exp)
447 NamedIndexArray_std_unary(log)
448 NamedIndexArray_std_unary(log10)
449 NamedIndexArray_std_unary(sin)
450 NamedIndexArray_std_unary(cos)
451 NamedIndexArray_std_unary(tan)
452 NamedIndexArray_std_unary(asin)
453 NamedIndexArray_std_unary(acos)
454 NamedIndexArray_std_unary(atan)
455 NamedIndexArray_std_binary(atan2)
456 NamedIndexArray_std_unary(sinh)
457 NamedIndexArray_std_unary(cosh)
458 NamedIndexArray_std_unary(tanh)
459 NamedIndexArray_std_unary(sqrt)
460 NamedIndexArray_std_unary(abs)
461 NamedIndexArray_std_unary(fabs)
462 NamedIndexArray_std_binary(max)
463 NamedIndexArray_std_binary(min)
464 NamedIndexArray_std_unary(ceil)
465 NamedIndexArray_std_unary(floor)
466 NamedIndexArray_std_binary(fmod)
467
468 #if __cplusplus >= 201103L
469 NamedIndexArray_std_unary(llabs)
470 NamedIndexArray_std_unary(imaxabs)
471 NamedIndexArray_fl_unary(fabs)
472 NamedIndexArray_fl_unary(exp)
473 NamedIndexArray_stdfl_unary(exp2)
474 NamedIndexArray_stdfl_unary(expm1)
475 NamedIndexArray_fl_unary(log)
476 NamedIndexArray_fl_unary(log10)
477 NamedIndexArray_stdfl_unary(log2)
478 NamedIndexArray_stdfl_unary(log1p)
479 NamedIndexArray_fl_unary(sqrt)
480 NamedIndexArray_stdfl_unary(cbrt)
481 NamedIndexArray_fl_unary(sin)
482 NamedIndexArray_fl_unary(cos)
483 NamedIndexArray_fl_unary(tan)
484 NamedIndexArray_fl_unary(asin)
485 NamedIndexArray_fl_unary(acos)
486 NamedIndexArray_fl_unary(atan)
487 NamedIndexArray_fl_unary(sinh)
488 NamedIndexArray_fl_unary(cosh)
489 NamedIndexArray_fl_unary(tanh)
490 NamedIndexArray_stdfl_unary(asinh)
491 NamedIndexArray_stdfl_unary(acosh)
492 NamedIndexArray_stdfl_unary(atanh)
493 NamedIndexArray_stdfl_unary(erf)
494 NamedIndexArray_stdfl_unary(erfc)
495 NamedIndexArray_stdfl_unary(tgamma)
496 NamedIndexArray_stdfl_unary(lgamma)
497 NamedIndexArray_fl_unary(ceil)
498 NamedIndexArray_fl_unary(floor)
499 NamedIndexArray_stdfl_unary(trunc)
500 NamedIndexArray_stdfl_unary(round)
501 NamedIndexArray_stdfl_unary(nearbyint)
502 NamedIndexArray_stdfl_unary(rint)
503
504 NamedIndexArray_fl_binary(pow)
505 NamedIndexArray_fl_binary(fmod)
506 NamedIndexArray_stdfl_binary(remainder)
507 NamedIndexArray_stdfl_binary(fmax)
508 NamedIndexArray_stdfl_binary(fmin)
509 NamedIndexArray_stdfl_binary(fdim)
510 NamedIndexArray_stdfl_binary(hypot)
511 NamedIndexArray_fl_binary(atan2)
512 #endif // __cplusplus >= 201103L
513
514 template <typename DataVector, typename SparseSizeVector>
515 class numeric_limits<NamedIndexArray<DataVector, SparseSizeVector> > :
516 public MetaPhysicL::raw_numeric_limits<NamedIndexArray<DataVector, SparseSizeVector>, DataVector> {};
517
518 } // namespace std
519
520
521 #endif // METAPHYSICL_NAMEDINDEXARRAY_H
522