Tomographer  v5.4
Tomographer C++ Framework Documentation
eigenutil.h
Go to the documentation of this file.
1 /* This file is part of the Tomographer project, which is distributed under the
2  * terms of the MIT license.
3  *
4  * The MIT License (MIT)
5  *
6  * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
7  * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a copy
10  * of this software and associated documentation files (the "Software"), to deal
11  * in the Software without restriction, including without limitation the rights
12  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13  * copies of the Software, and to permit persons to whom the Software is
14  * furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included in
17  * all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  * SOFTWARE.
26  */
27 
28 #ifndef TOMOGRAPHER_TOOLS_EIGENUTIL_H
29 #define TOMOGRAPHER_TOOLS_EIGENUTIL_H
30 
38 #include <complex>
39 #include <vector>
40 
41 #include <Eigen/Core>
42 #include <Eigen/Dense>
43 #include <Eigen/StdVector>
44 
45 #include <boost/serialization/serialization.hpp>
46 #include <boost/serialization/split_free.hpp>
47 #include <boost/serialization/complex.hpp>
48 
49 
50 // -----------------------------------------------------------------------------
51 // Eigen and std::vector
52 // -----------------------------------------------------------------------------
53 
54 namespace Tomographer {
55 namespace Tools {
56 
71 template<typename EigenType>
73 {
75 };
76 
77 } // namespace Tools
78 } // namespace Tomographer
79 
80 // -----------------------------------------------------------------------------
81 // Random matrices in Eigen
82 // -----------------------------------------------------------------------------
83 
84 namespace Tomographer {
85 namespace Tools {
86 
87 namespace tomo_internal {
88 
89  template<typename Rng, typename RndDist, typename Scalar>
90  struct random_generator
91  {
92  typedef Scalar result_type;
93 
94  Rng & rng;
95  RndDist & rnddist;
96 
97  random_generator(Rng & rng_, RndDist & rnddist_)
98  : rng(rng_), rnddist(rnddist_)
99  {
100  }
101 
102  template<typename Index>
103  inline const result_type operator() (Index, Index = 0) const {
104  return rnddist(rng);
105  }
106  };
107 
108  template<typename Rng, typename RndDist, typename RealScalar>
109  struct random_generator<Rng, RndDist, std::complex<RealScalar> >
110  {
111  typedef std::complex<RealScalar> result_type;
112 
113  Rng & rng;
114  RndDist & rnddist;
115 
116  random_generator(Rng & rng_, RndDist & rnddist_)
117  : rng(rng_), rnddist(rnddist_)
118  {
119  }
120 
121  template<typename Index>
122  inline const result_type operator() (Index, Index = 0) const {
123  return result_type(rnddist(rng), rnddist(rng));
124  }
125  };
126 } // end namespace tomo_internal
127 
128 } // namespace Tools
129 } // namespace Tomographer
130 
131 namespace Eigen {
132  namespace internal {
134  template<typename Rng, typename RndDist, typename Scalar>
135  struct functor_traits<Tomographer::Tools::tomo_internal::random_generator<Rng, RndDist, Scalar> >
136  { enum { Cost = 50 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
137  }
138 } // end namespace Eigen
139 
140 
141 namespace Tomographer {
142 namespace Tools {
143 
144 
153 template<typename Der, typename Rng, typename RndDist, typename... IndexTypes>
154 inline auto denseRandom(Rng & rng, RndDist &rnddist, IndexTypes... sizes)
155  -> const Eigen::CwiseNullaryOp<
156  tomo_internal::random_generator<Rng, RndDist, typename Eigen::internal::traits<Der>::Scalar>,
157  Der
158  >
159 {
160  typedef typename Der::Scalar Scalar;
161 
163  sizes..., tomo_internal::random_generator<Rng, RndDist, Scalar>(rng, rnddist)
164  );
165 }
166 
167 
168 
169 // ---------------------------
170 
171 
172 
173 namespace tomo_internal {
175  template<typename Scalar, typename IndexType>
176  struct can_basis_vec_generator
177  {
178  typedef Scalar result_type;
179 
180  const IndexType k;
181  const IndexType j;
182 
183  can_basis_vec_generator(IndexType k_, IndexType j_ = 0)
184  : k(k_), j(j_)
185  {
186  }
187 
188  inline const result_type operator() (IndexType a, IndexType b = 0) const {
189  return (a == k) && (b == j) ? result_type(1) : result_type(0);
190  }
191  };
192 } // namespace tomo_internal
193 } // namespace Tools
194 } // namespace Tomographer
195 namespace Eigen {
196  namespace internal {
198  template<typename Scalar, typename IndexType>
199  struct functor_traits<Tomographer::Tools::tomo_internal::can_basis_vec_generator<Scalar, IndexType> >
200  { enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = true }; };
201  template<typename Scalar, typename IndexType>
202  struct functor_has_linear_access<Tomographer::Tools::tomo_internal::can_basis_vec_generator<Scalar, IndexType> >
203  { enum { ret = 0 }; };
204  }
205 } // end namespace Eigen
206 
207 namespace Tomographer {
208 namespace Tools {
209 
213 template<typename Der, typename IndexType>
214 inline auto canonicalBasisVec(IndexType k, IndexType size)
215  -> const Eigen::CwiseNullaryOp<
216  tomo_internal::can_basis_vec_generator<typename Eigen::internal::traits<Der>::Scalar, IndexType>,
217  Der
218  >
219 {
220  typedef typename Der::Scalar Scalar;
221 
223  size, tomo_internal::can_basis_vec_generator<Scalar, IndexType>(k)
224  );
225 }
226 
231 template<typename Der, typename IndexType>
232 inline auto canonicalBasisVec(IndexType k, IndexType j, IndexType rows, IndexType cols)
233  -> const Eigen::CwiseNullaryOp<
234  tomo_internal::can_basis_vec_generator<typename Eigen::internal::traits<Der>::Scalar, IndexType>,
235  Der
236  >
237 {
238  typedef typename Der::Scalar Scalar;
239 
241  rows, cols, tomo_internal::can_basis_vec_generator<Scalar, IndexType>(k, j)
242  );
243 }
244 
245 
246 
247 
248 // ---------------------------
249 
250 
251 
252 namespace tomo_internal {
254  template<typename Scalar>
255  struct powers_of_two_generator
256  {
257  typedef typename Eigen::NumTraits<Scalar>::Real result_type;
258 
259  powers_of_two_generator() { }
260 
261  template<typename IndexType>
262  inline const result_type operator() (IndexType a) const {
263  return std::ldexp(result_type(1), (int)a);
264  }
265 
266  // Don't expose an operator()(i,j) method, because otherwise Eigen might think that we
267  // don't have linear access
268  // https://eigen.tuxfamily.org/dox/NullaryFunctors_8h_source.html#l00147
269  //
270  // template<typename IndexType>
271  // inline const result_type operator() (IndexType a, IndexType b) const {
272  // eigen_assert(b == 0 && "powers_of_two_generator may only be used with 1-D objects or with linear access!");
273  // (void)b; // silence unused variable warning if eigen_assert is optimized out
274  // return std::ldexp(result_type(1), a);
275  // }
276 
277  };
278 } // namespace tomo_internal
279 } // namespace Tools
280 } // namespace Tomographer
281 namespace Eigen {
282  namespace internal {
284  template<typename Scalar>
285  struct functor_traits<Tomographer::Tools::tomo_internal::powers_of_two_generator<Scalar> >
286  { enum { Cost = 8 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = true }; };
287  }
288 } // end namespace Eigen
289 
290 namespace Tomographer {
291 namespace Tools {
292 
293 
306 template<typename Der, typename... IndexTypes>
307 inline auto powersOfTwo(IndexTypes... sizes)
308  -> const Eigen::CwiseNullaryOp<
309  tomo_internal::powers_of_two_generator<typename Eigen::internal::traits<Der>::Scalar>,
310  Der
311  >
312 {
313  typedef typename Der::Scalar Scalar;
314 
316  sizes..., tomo_internal::powers_of_two_generator<Scalar>()
317  );
318 }
319 
320 
321 
322 
323 
324 // -----------------------------------------------------------------------------
325 
341 template<int RowFactorCTime, int ColFactorCTime, typename Derived,
342  typename std::enable_if<(RowFactorCTime == Eigen::Dynamic || ColFactorCTime == Eigen::Dynamic),
343  bool>::type dummy = true>
344 inline auto replicated(const Eigen::DenseBase<Derived> & x, int row_factor, int col_factor)
346 {
347  eigen_assert(RowFactorCTime == Eigen::Dynamic || row_factor == RowFactorCTime);
348  eigen_assert(ColFactorCTime == Eigen::Dynamic || col_factor == ColFactorCTime);
349  return x.replicate(row_factor, col_factor);
350 }
355 template<int RowFactorCTime, int ColFactorCTime, typename Derived,
356  typename std::enable_if<(RowFactorCTime != Eigen::Dynamic && ColFactorCTime != Eigen::Dynamic),
357  bool>::type dummy2 = true>
358 inline auto replicated(const Eigen::DenseBase<Derived> & x, int row_factor, int col_factor)
360 {
361  eigen_assert(row_factor == RowFactorCTime); (void)row_factor; // "unused argument" warning
362  eigen_assert(col_factor == ColFactorCTime); (void)col_factor; // "unused argument" warning
363  return x.template replicate<RowFactorCTime, ColFactorCTime>();
364 }
365 
366 
367 
368 
369 } // namespace Tools
370 } // namespace Tomographer
371 
372 
373 
374 
375 
376 // -----------------------------------------------------------------------------
377 // Eigen Matrix/Array and boost::serialization
378 // -----------------------------------------------------------------------------
379 namespace boost {
380 namespace serialization {
381 
382 template<typename Archive, typename Derived>
383 void eigen_helper_save(Archive & a, const Derived & m, unsigned int /* version */)
384 {
385  Eigen::Index rows = m.rows();
386  Eigen::Index cols = m.cols();
387  a << rows;
388  a << cols;
389  for (Eigen::Index i = 0; i < m.rows(); ++i) {
390  for (Eigen::Index j = 0; j < m.cols(); ++j) {
391  a << m(i,j);
392  }
393  }
394 }
395 
396 template<typename Archive, typename Derived>
397 void eigen_helper_load(Archive & a, Derived & m, unsigned int /* version */)
398 {
399  Eigen::Index rows, cols;
400  a >> rows;
401  a >> cols;
402  m.resize(rows, cols); // legal also on fixed-size matrices
403  for (Eigen::Index i = 0; i < m.rows(); ++i) {
404  for (Eigen::Index j = 0; j < m.cols(); ++j) {
405  a >> m(i,j);
406  }
407  }
408 
409 }
410 
411 
412 template<typename Archive,
413  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
415  const unsigned int version)
416 {
417  eigen_helper_save(a, m, version);
418 }
419 template<typename Archive,
420  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
422  const unsigned int version)
423 {
424  eigen_helper_save(a, m, version);
425 }
426 template<typename Archive,
427  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
429  const unsigned int version)
430 {
431  eigen_helper_load(a, m, version);
432 }
433 template<typename Archive,
434  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
436  const unsigned int version)
437 {
438  eigen_helper_load(a, m, version);
439 }
440 
441 template<typename Archive,
442  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
443 inline void serialize(Archive & a, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
444  const unsigned int version)
445 {
446  split_free(a, m, version);
447 }
448 template<typename Archive,
449  typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
450 inline void serialize(Archive & a, Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & m,
451  const unsigned int version)
452 {
453  split_free(a, m, version);
454 }
455 
456 
457 } // serialization
458 } // boost
459 
460 
461 
462 
463 
464 
465 #endif
auto canonicalBasisVec(IndexType k, IndexType size) -> const Eigen::CwiseNullaryOp< tomo_internal::can_basis_vec_generator< typename Eigen::internal::traits< Der >::Scalar, IndexType >, Der >
Expression for the k-th canonical basis vector of given dimension.
Definition: eigenutil.h:214
static const CwiseNullaryOp< CustomNullaryOp, Derived > NullaryExpr(Index rows, Index cols, const CustomNullaryOp &func)
Base namespace for the Tomographer project.
Definition: densellh.h:45
STL class.
Use this if you need an std::vector of any Eigen type.
Definition: eigenutil.h:72
T ldexp(T... args)
auto powersOfTwo(IndexTypes... sizes) -> const Eigen::CwiseNullaryOp< tomo_internal::powers_of_two_generator< typename Eigen::internal::traits< Der >::Scalar >, Der >
Expression for a 1-D expression of powers of two.
Definition: eigenutil.h:307
auto replicated(const Eigen::DenseBase< Derived > &x, int row_factor, int col_factor) -> const Eigen::Replicate< Derived, Eigen::Dynamic, Eigen::Dynamic >
Replicate a Eigen Dense object; same call for compile-time & run-time dimensions. ...
Definition: eigenutil.h:344
auto denseRandom(Rng &rng, RndDist &rnddist, IndexTypes... sizes) -> const Eigen::CwiseNullaryOp< tomo_internal::random_generator< Rng, RndDist, typename Eigen::internal::traits< Der >::Scalar >, Der >
a matrix populated with random entries using C++&#39;s random framework
Definition: eigenutil.h:154
STL class.
const Replicate< Derived, RowFactor, ColFactor > replicate() const