Tomographer  v5.4
Tomographer C++ Framework Documentation
ezmatio.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_EZMATIO_H
29 #define _TOMOGRAPHER_TOOLS_EZMATIO_H
30 
31 #include <cstdint>
32 #include <cerrno>
33 
34 #include <complex>
35 #include <string>
36 #include <vector>
37 #include <algorithm>
38 #include <iostream>
39 #include <stdexcept>
40 #include <sstream>
41 #include <initializer_list>
42 #include <utility>
43 
44 extern "C" {
45 # include <matio.h>
46 }
47 
48 #include <Eigen/Core>
49 
52 #include <tomographer/tools/fmt.h>
54 
55 
63 namespace Tomographer {
64 // namespace doc in doc/doxdocs/namespaces.cxx
65 namespace MAT {
66 
67 
70 class TOMOGRAPHER_EXPORT Exception : public std::exception
71 {
72  std::string p_heading;
73  std::string p_message;
74  std::string p_final;
75 public:
76  Exception(std::string msg) : p_message(msg) { update_final(); }
77  Exception(std::string heading, std::string msg) : p_heading(heading), p_message(msg) { update_final(); }
78  virtual ~Exception() noexcept { }
79 
80  virtual const char * what() const noexcept { return p_final.c_str(); }
81 
82 protected:
83  inline void setHeading(const std::string heading) {
84  p_heading = heading;
85  update_final();
86  }
87  inline void setMessage(const std::string message) {
88  p_message = message;
89  update_final();
90  }
91 private:
92  inline void update_final() {
93  p_final = p_heading + p_message;
94  }
95 };
96 
99 class TOMOGRAPHER_EXPORT VarError : public Exception
100 {
101 public:
102  VarError(std::string msg) : Exception("", msg) { }
103  VarError(std::string varname, std::string msg) : Exception(heading(varname), msg) { }
104  virtual ~VarError() noexcept { }
105 
106  void setVarName(std::string varname) { setHeading(heading(varname)); }
107 
108 private:
109  static std::string heading(std::string varname) { return "Variable " + varname + ": "; }
110 };
111 
113 class TOMOGRAPHER_EXPORT VarReadError : public VarError
114 {
115 public:
116  VarReadError(const std::string varname)
117  : VarError(varname, "Can't read variable")
118  {
119  }
120  virtual ~VarReadError() noexcept { }
121 };
122 
124 class TOMOGRAPHER_EXPORT VarTypeError : public VarError
125 {
126 public:
127  VarTypeError(const std::string varname, const std::string msg)
128  : VarError(varname, msg)
129  {
130  }
131  virtual ~VarTypeError() noexcept { }
132 };
133 
135 class TOMOGRAPHER_EXPORT VarMatTypeError : public VarError {
136 public:
137  VarMatTypeError(const std::string msg)
138  : VarError(msg)
139  {
140  }
141  VarMatTypeError(const std::string varname, const std::string msg)
142  : VarError(varname, msg)
143  {
144  }
145  virtual ~VarMatTypeError() noexcept { }
146 };
147 
149 class TOMOGRAPHER_EXPORT FileOpenError : public Exception
150 {
151 public:
152  FileOpenError(const std::string fname, const std::string errmsg = std::string())
153  : Exception(heading(fname), "Error opening file" + (errmsg.size() ? ": "+errmsg : ""))
154  {
155  }
156  virtual ~FileOpenError() noexcept { }
157 
158  void setFileName(const std::string fname) {
159  setHeading(heading(fname));
160  }
161 
162 private:
163  static std::string heading(std::string fname) { return "File `"+fname+"`: "; }
164 };
165 
166 
168 class TOMOGRAPHER_EXPORT InvalidIndexError : public Exception {
169 public:
170  InvalidIndexError(const std::string msg) : Exception("Invalid index: ", msg) { }
171  virtual ~InvalidIndexError() noexcept { }
172 };
173 
174 
175 
176 
177 
178 class Var;
179 
180 
192 class TOMOGRAPHER_EXPORT File
193 {
194 public:
200  File(const std::string fname)
201  {
202  errno = 0;
203  p_matfp = Mat_Open(fname.c_str(), MAT_ACC_RDONLY);
204  if ( p_matfp == NULL ) {
205  throw FileOpenError(fname, strerror(errno));
206  }
207  }
208  File(const File& other) = delete;
209 
215  File(File&& other)
216  {
217  // move constructor. Steal the other object's MAT file handle
218  p_matfp = other.p_matfp;
219  other.p_matfp = NULL;
220  }
221 
222  ~File()
223  {
224  if ( p_matfp != NULL ) {
225  Mat_Close(p_matfp);
226  }
227  }
228 
238  inline Var var(const std::string varname, bool load_data = true);
239 
240 
246  inline std::vector<Var> getVarInfoList();
247 
253  mat_t * getMatPtr()
254  {
255  return p_matfp;
256  }
257 
263  File& operator=(File&& other)
264  {
265  // move assignment operator. Steal the other object's MAT file handle
266  p_matfp = other.p_matfp;
267  other.p_matfp = NULL;
268  return *this;
269  }
270  File& operator=(const File& other) = delete;
271 
272 private:
273  mat_t *p_matfp;
274 };
275 
276 
282 template<typename It, typename ValueType = typename std::iterator_traits<It>::value_type>
283 inline ValueType getNumEl(It begin, It end)
284 {
285  ValueType n = 1;
286  for (It it = begin; it != end; ++it) {
287  n *= (*it);
288  }
289  return n;
290 }
291 
292 
306 class TOMOGRAPHER_EXPORT DimList : public std::vector<int>
307 {
308 public:
310  DimList() : std::vector<int>() { }
311 
317  template<typename VectorType>
318  DimList(VectorType&& dims)
319  : std::vector<int>(std::forward<VectorType>(dims))
320  {
321  }
335  : std::vector<int>(init)
336  {
337  }
343  template<class It>
344  DimList(It b, It e) : std::vector<int>(b, e)
345  {
346  }
347 
352  inline int numel() const { return getNumEl(begin(), end()); }
353 
359  inline int ndims() const { return (int)size(); }
360 
364  inline bool matchesWanted(const DimList& wanted) const
365  {
366  if (size() != wanted.size()) {
367  return false;
368  }
369  for (std::size_t k = 0; k < size(); ++k) {
370  if (wanted[k] >= 0 && this->operator[](k) != wanted[k]) {
371  return false;
372  }
373  }
374  return true;
375  }
376 
382  DimList& operator<<(int dim) {
383  push_back(dim);
384  return *this;
385  }
391  DimList& operator<<(const std::vector<int>& moredims) {
392  insert(end(), moredims.begin(), moredims.end());
393  return *this;
394  }
395 };
396 
398 inline std::ostream& operator<<(std::ostream& out, const DimList& dlist)
399 {
400  out << "[";
401  for (DimList::const_iterator it = dlist.begin(); it != dlist.end(); ++it) {
402  if (it != dlist.begin()) {
403  out << " ";
404  }
405  if (*it == -1) {
406  // usually meant for arbitrary size requirement. In any case, it is not a physical
407  // array size.
408  out << '*';
409  } else {
410  out << *it;
411  }
412  }
413  out << "]";
414  return out;
415 }
416 
439 template<bool IsRowMajor_ = false>
440 class TOMOGRAPHER_EXPORT IndexList : public std::vector<int>
441 {
442 public:
444  static constexpr bool IsRowMajor = IsRowMajor_;
447 
453  IndexList(const std::vector<int>& dims = std::vector<int>(), int linearindex = -1)
454  : std::vector<int>(dims.size()), p_dims(dims)
455  {
456  if (getNumEl(dims.begin(), dims.end()) <= 0) {
457  throw InvalidIndexError("zero-sized array given by dimension list");
458  }
459  if (linearindex >= 0) {
460  setLinearIndex(linearindex);
461  }
462  }
469  // the reason we dont have template<typename... VectorArgs> for std::vector<>
470  // constructor is that g++ 4.6 if I remember correctly can't unpack arguments into a
471  // fixed number of arguments (TODO: CHECK THIS!)
472  template<typename VectorIntInitializer>
473  IndexList(const std::vector<int>& dims, VectorIntInitializer&& index)
474  : std::vector<int>(std::forward<VectorIntInitializer>(index)), p_dims(dims)
475  {
476  if (getNumEl(dims.begin(), dims.end()) <= 0) {
477  throw InvalidIndexError("zero-sized array given by dimension list");
478  }
479  }
480 
488  TOMOGRAPHER_ENABLED_IF(IsRowMajor)
489  void setLinearIndex(int linearindex)
490  {
491  const int ndims = (int)p_dims.size();
492  for (int k = ndims-1; k >= 0; --k) {
493  this->at((std::size_t)k) = linearindex % p_dims[(std::size_t)k];
494  linearindex /= p_dims[(std::size_t)k]; // integer division
495  }
496  }
504  TOMOGRAPHER_ENABLED_IF(!IsRowMajor)
505  void setLinearIndex(int linearindex)
506  {
507  const int ndims = (int)p_dims.size();
508  for (int k = 0; k < ndims; ++k) {
509  this->at((std::size_t)k) = linearindex % p_dims[(std::size_t)k];
510  linearindex /= p_dims[(std::size_t)k]; // integer division
511  // std::cout << "k = " << k << "; p_dims = " << p_dims << "; at(k) = " << this->at(k)
512  // << "; linearindex=" << linearindex << "\n";
513  }
514  }
515 
520  TOMOGRAPHER_ENABLED_IF(IsRowMajor)
521  int linearIndex() const
522  {
523  int linindex = 0;
524  for (int k = 0; k < (int)p_dims.size(); ++k) {
525  linindex *= p_dims[(std::size_t)k];
526  linindex += this->at((std::size_t)k);
527  }
528  return linindex;
529  }
534  TOMOGRAPHER_ENABLED_IF(!IsRowMajor)
535  int linearIndex() const
536  {
537  int linindex = 0;
538  for (int k = (int)p_dims.size()-1; k >= 0; --k) {
539  linindex *= p_dims[(std::size_t)k];
540  linindex += this->at((std::size_t)k);
541  }
542  return linindex;
543  }
544 
547  inline const std::vector<int> & index() const
548  {
549  return *this;
550  }
551 
556  inline const std::vector<int> & dims() const
557  {
558  return p_dims;
559  }
560 
568  IndexList& operator<<(int ind) {
569  push_back(ind);
570  return *this;
571  }
572 
581  IndexList& operator<<(const std::vector<int>& moredims) {
582  insert(end(), moredims.begin(), moredims.end());
583  return *this;
584  }
585 
586 
587 
588  // for internal use:
589 
590  template<bool IsRowMajor2 = false>
591  static const std::vector<int> & forward_index(const IndexList<IsRowMajor2> & index)
592  {
593  return index;
594  }
595  template<bool IsRowMajor2 = false>
596  static std::vector<int> && forward_index(IndexList<IsRowMajor2> && index)
597  {
598  return std::move(index);
599  }
600 
601 private:
602  const std::vector<int> p_dims;
603 };
604 
605 
607 template<bool IsRowMajor>
608 inline std::ostream& operator<<(std::ostream& str, const IndexList<IsRowMajor> & indexlist)
609 {
610  str << "[";
611  for (std::size_t j = 0; j < indexlist.size(); ++j) {
612  str << indexlist[j];
613  if (j < indexlist.size()) {
614  str << ", ";
615  } else {
616  str << ";";
617  }
618  }
619  return str << "==" << indexlist.linearIndex() << "]";
620 }
621 
622 
626 template<bool IsRowMajor_ = false, typename IntType_ = int>
627 class TOMOGRAPHER_EXPORT IndexListIterator
628 {
629 public:
630  static constexpr bool IsRowMajor = IsRowMajor_;
631  typedef IntType_ IntType;
633 
634 private:
635  const std::vector<IntType> p_dims;
636  const IntType p_numel;
637 
638  std::vector<IntType> p_index;
639  IntType p_linearIndex;
640 
641 public:
643  : p_dims(dims),
644  p_numel(getNumEl(dims.begin(), dims.end())),
645  p_index(dims.size(), 0),
646  p_linearIndex(0)
647  {
648  if (p_numel <= 0) {
650  ss << "Invalid dimensions: " << p_dims;
651  throw InvalidIndexError(ss.str());
652  }
653  }
654 
655  inline const std::vector<IntType>& index() const
656  {
657  return p_index;
658  }
659 
660  inline IntType linearIndex() const
661  {
662  return p_linearIndex;
663  }
664 
665  inline IntType numel() const
666  {
667  return p_numel;
668  }
669 
670 
671  TOMOGRAPHER_ENABLED_IF(IsRowMajor)
672  IntType increment()
673  {
674  for (int k = (int)p_dims.size() - 1; k >= 0; --k) {
675  p_index[(std::size_t)k]++;
676  if (p_index[(std::size_t)k] < p_dims[(std::size_t)k]) {
677  // if this increment succeeded and stays in range, ok and stop.
678  break;
679  } else {
680  // otherwise continue the loop and increment the next value, while resetting
681  // this one to zero.
682  p_index[(std::size_t)k] = 0;
683  }
684  }
685  return ++p_linearIndex;
686  }
687 
688  TOMOGRAPHER_ENABLED_IF(!IsRowMajor)
689  IntType increment()
690  {
691  for (std::size_t k = 0; k < p_dims.size(); ++k) {
692  p_index[k]++;
693  if (p_index[k] < p_dims[k]) {
694  // if this increment succeeded and stays in range, ok and stop.
695  break;
696  } else {
697  // otherwise continue the loop and increment the next value, while resetting
698  // this one to zero.
699  p_index[k] = 0;
700  }
701  }
702  return ++p_linearIndex;
703  }
704 
705  IntType operator++() {
706  return increment();
707  }
708 
709  bool valid() const
710  {
711  return Tools::isPositive<IntType>(p_linearIndex) && p_linearIndex < p_numel;
712  }
713 
714 
715  // for internal use:
716 
717  template<bool IsRowMajor2 = false, typename IntType2 = int>
718  static const std::vector<IntType2> &
719  forward_index(const IndexListIterator<IsRowMajor2,IntType2> & index)
720  {
721  return index.p_index;
722  }
723  template<bool IsRowMajor2 = false, typename IntType2 = int>
724  static std::vector<IntType2> &&
725  forward_index(IndexListIterator<IsRowMajor2,IntType2> && index)
726  {
727  return std::move(index.p_index);
728  }
729 
730 };
731 
732 
734 template<bool IsRowMajor, typename IntType>
735 inline std::ostream& operator<<(std::ostream& str, const IndexListIterator<IsRowMajor, IntType> & indexlistit)
736 {
737  std::vector<int> indexlist{indexlistit.index()};
738 
739  str << "[";
740  for (std::size_t j = 0; j < indexlist.size(); ++j) {
741  str << indexlist[j];
742  if (j < indexlist.size()) {
743  str << ", ";
744  } else {
745  str << ";";
746  }
747  }
748  return str << "==" << indexlistit.linearIndex() << "]";
749 }
750 
751 
752 
753 
754 
801 template<typename T, typename Enabled = void>
802 class TOMOGRAPHER_EXPORT VarValueDecoder
803 {
804 public:
810  typedef T RetType;
811 
818  static inline void checkShape(const Var & var)
819  {
820  (void)var; // silence "unused parameter" warnings
821  }
822 
837  static inline RetType decodeValue(const Var & var);
838 };
839 
840 
841 namespace tomo_internal {
842 
843  template<typename T, typename Res = void>
844  struct enable_if_is_type { typedef T check_type; typedef Res type; };
845 
846  // has_params_member<T>::value is true if T has a type member named Params,
847  // i.e. typename T::Params, else it is false.
848  template<typename T, typename Enabled = void>
849  struct has_params_member {
850  enum { value = 0 };
851  };
852  // if T has no Params member, following will fail by SFINAE (the default arguemnt of
853  // parameter #2 always stays int)
854  template<typename T>
855  struct has_params_member<T, typename enable_if_is_type<typename T::Params>::type >
856  {
857  enum { value = 1 };
858  };
859 
860 template<typename Scalar>
861 struct isComplex { // use Eigen's existing implementation
863 };
864 
865 
866 }; // namespace tomo_internal
867 
868 
869 
870 class Var;
871 
877 template<typename T>
878 inline typename VarValueDecoder<T>::RetType value(const Var& var)
879 {
882 }
883 
892 template<typename T, TOMOGRAPHER_ENABLED_IF_TMPL( tomo_internal::has_params_member<VarValueDecoder<T> >::value )>
893 inline typename VarValueDecoder<T>::RetType value(const Var& var,
894  const typename VarValueDecoder<T>::Params & params)
895 {
896  VarValueDecoder<T>::checkShape(var, params);
897  return VarValueDecoder<T>::decodeValue(var, params);
898 }
899 
900 
901 
918 class TOMOGRAPHER_EXPORT Var
919 {
920 private:
921  struct VarData {
922  VarData() : refcount(0), p_matvar(NULL) { }
923  int refcount;
924  matvar_t * p_matvar;
925  std::string p_varname;
926  };
927  VarData *p_vardata;
928 
929  Var(matvar_t * varinfo) // take over the struct ownership
930  {
931  if (varinfo == NULL) {
932  throw VarReadError("<unknown>");
933  }
934  p_vardata = new VarData;
935  p_vardata->refcount++;
936  p_vardata->p_matvar = varinfo;
937  p_vardata->p_varname = std::string(varinfo->name);
938  }
939 
940 public:
951  Var(File & matf, const std::string& varname, bool load_data = true)
952  {
953  p_vardata = new VarData;
954  p_vardata->refcount++;
955  if (load_data) {
956  p_vardata->p_matvar = Mat_VarRead(matf.getMatPtr(), varname.c_str());
957  } else {
958  p_vardata->p_matvar = Mat_VarReadInfo(matf.getMatPtr(), varname.c_str());
959  }
960  if (p_vardata->p_matvar == NULL) {
961  delete p_vardata;
962  p_vardata = NULL;
963  throw VarReadError(varname);
964  }
965  p_vardata->p_varname = varname;
966  }
967 
969  Var(const Var& copy)
970  {
971  p_vardata = copy.p_vardata;
972  p_vardata->refcount++;
973  }
974 
976  Var(Var&& other)
977  {
978  // steal the data
979  p_vardata = other.p_vardata;
980  other.p_vardata = NULL;
981  }
982 
983  ~Var() {
984  if (p_vardata != NULL) {
985  p_vardata->refcount--;
986  if (!p_vardata->refcount) {
987  Mat_VarFree(p_vardata->p_matvar);
988  delete p_vardata;
989  }
990  }
991  }
992 
999  static Var takeOver(matvar_t * varinfo)
1000  {
1001  return Var(varinfo);
1002  }
1003 
1008  inline const std::string & varName() const
1009  {
1010  tomographer_assert(p_vardata != NULL);
1011  return p_vardata->p_varname;
1012  }
1013 
1018  inline int ndims() const
1019  {
1020  tomographer_assert(p_vardata != NULL);
1021  return p_vardata->p_matvar->rank;
1022  }
1027  inline DimList dims() const
1028  {
1029  tomographer_assert(p_vardata != NULL);
1030  return DimList(&p_vardata->p_matvar->dims[0],
1031  &p_vardata->p_matvar->dims[p_vardata->p_matvar->rank]);
1032  }
1037  inline int numel() const
1038  {
1039  tomographer_assert(p_vardata != NULL);
1040  return dims().numel();
1041  }
1050  inline bool isComplex() const
1051  {
1052  tomographer_assert(p_vardata != NULL);
1053  return p_vardata->p_matvar->isComplex;
1054  }
1059  inline bool isSquareMatrix() const
1060  {
1061  tomographer_assert(p_vardata != NULL);
1062  return p_vardata->p_matvar->rank == 2 && p_vardata->p_matvar->dims[0] == p_vardata->p_matvar->dims[1];
1063  }
1064 
1070  inline bool hasData() const
1071  {
1072  tomographer_assert(p_vardata != NULL);
1073  return (p_vardata->p_matvar->data != NULL);
1074  }
1075 
1115  template<typename T>
1116  inline typename VarValueDecoder<T>::RetType value() const
1117  {
1118  return Tomographer::MAT::value<T>(*this);
1119  }
1130  template<typename T, TOMOGRAPHER_ENABLED_IF_TMPL( tomo_internal::has_params_member<VarValueDecoder<T> >::value )>
1131  inline typename VarValueDecoder<T>::RetType value(const typename VarValueDecoder<T>::Params & params)
1132  {
1133  return Tomographer::MAT::value<T>(*this, params);
1134  }
1135 
1141  const matvar_t * getMatvarPtr() const
1142  {
1143  tomographer_assert(p_vardata != NULL);
1144  return p_vardata->p_matvar;
1145  }
1146 
1150  Var& operator=(Var&& other)
1151  {
1152  // steal the data
1153  p_vardata = other.p_vardata;
1154  other.p_vardata = NULL;
1155  return *this;
1156  }
1158  Var& operator=(const Var& copy)
1159  {
1160  p_vardata = copy.p_vardata;
1161  p_vardata->refcount++;
1162  return *this;
1163  }
1164 };
1165 
1166 
1167 // ---------------
1168 
1169 
1170 inline Var File::var(const std::string varname, bool load_data)
1171 {
1172  return Var(*this, varname, load_data);
1173 }
1174 
1176 {
1177  std::vector<Var> varlist;
1178 
1179  Mat_Rewind(p_matfp);
1180  matvar_t * p = NULL;
1181  while ( (p = Mat_VarReadNextInfo(p_matfp)) != NULL ) {
1182  varlist.push_back(Var::takeOver(p));
1183  }
1184 
1185  return varlist;
1186 }
1187 
1188 
1189 template<typename T, typename Enabled>
1190 // static
1192 {
1193  throw std::runtime_error(std::string("Not Implemented: Please specialize MAT::VarValueDecoder<> for ")
1194  + typeid(T).name() + " to decode variable " + var.varName());
1195 }
1196 
1197 
1198 
1199 
1200 
1201 // =============================================================================
1202 // Some basic utilities for dealing with MATLAB variables
1203 // =============================================================================
1204 
1205 
1206 
1218 template<int MatTypeId = -1> struct MatType { };
1220 template<> struct MatType<MAT_T_DOUBLE> { typedef double Type; };
1222 template<> struct MatType<MAT_T_SINGLE> { typedef float Type; };
1224 template<> struct MatType<MAT_T_INT64> { typedef int64_t Type; };
1226 template<> struct MatType<MAT_T_INT32> { typedef int32_t Type; };
1228 template<> struct MatType<MAT_T_INT16> { typedef int16_t Type; };
1230 template<> struct MatType<MAT_T_INT8> { typedef int8_t Type; };
1232 template<> struct MatType<MAT_T_UINT64> { typedef uint64_t Type; };
1234 template<> struct MatType<MAT_T_UINT32> { typedef uint32_t Type; };
1236 template<> struct MatType<MAT_T_UINT16> { typedef uint16_t Type; };
1238 template<> struct MatType<MAT_T_UINT8> { typedef uint8_t Type; };
1239 
1240 
1247 #define MAT_SWITCH_REAL_TYPE(typ, ...) \
1248  do { switch (typ) { \
1249  case MAT_T_DOUBLE: { typedef typename MatType<MAT_T_DOUBLE>::Type Type; { __VA_ARGS__; } break; } \
1250  case MAT_T_SINGLE: { typedef typename MatType<MAT_T_SINGLE>::Type Type; { __VA_ARGS__; } break; } \
1251  case MAT_T_INT64: { typedef typename MatType<MAT_T_INT64>::Type Type; { __VA_ARGS__; } break; } \
1252  case MAT_T_INT32: { typedef typename MatType<MAT_T_INT32>::Type Type; { __VA_ARGS__; } break; } \
1253  case MAT_T_INT16: { typedef typename MatType<MAT_T_INT16>::Type Type; { __VA_ARGS__; } break; } \
1254  case MAT_T_INT8: { typedef typename MatType<MAT_T_INT8>::Type Type; { __VA_ARGS__; } break; } \
1255  case MAT_T_UINT64: { typedef typename MatType<MAT_T_UINT64>::Type Type; { __VA_ARGS__; } break; } \
1256  case MAT_T_UINT32: { typedef typename MatType<MAT_T_UINT32>::Type Type; { __VA_ARGS__; } break; } \
1257  case MAT_T_UINT16: { typedef typename MatType<MAT_T_UINT16>::Type Type; { __VA_ARGS__; } break; } \
1258  case MAT_T_UINT8: { typedef typename MatType<MAT_T_UINT8>::Type Type; { __VA_ARGS__; } break; } \
1259  default: \
1260  throw VarMatTypeError( streamstr("Uknown/unsupported encoded type from matio: " \
1261  <<(typ)) ); \
1262  } \
1263  } while (false)
1264 
1271 #define MAT_SWITCH_COMPLEX_TYPE(typ, ...) \
1272  do { switch (typ) { \
1273  case MAT_T_DOUBLE: { typedef std::complex<typename MatType<MAT_T_DOUBLE>::Type> Type; { __VA_ARGS__; } break; } \
1274  case MAT_T_SINGLE: { typedef std::complex<typename MatType<MAT_T_SINGLE>::Type> Type; { __VA_ARGS__; } break; } \
1275  default: \
1276  throw VarMatTypeError( streamstr("Uknown/unsupported encoded type from matio: " \
1277  <<(typ)) ); \
1278  } \
1279  } while (false)
1280 
1281 
1289 #define MAT_SWITCH_TYPE(matvar_ptr, ...) \
1290  do { \
1291  if (!(matvar_ptr)->isComplex) { \
1292  switch ((matvar_ptr)->data_type) { \
1293  case MAT_T_DOUBLE: { typedef typename MatType<MAT_T_DOUBLE>::Type Type; { __VA_ARGS__; } break; } \
1294  case MAT_T_SINGLE: { typedef typename MatType<MAT_T_SINGLE>::Type Type; { __VA_ARGS__; } break; } \
1295  case MAT_T_INT64: { typedef typename MatType<MAT_T_INT64>::Type Type; { __VA_ARGS__; } break; } \
1296  case MAT_T_INT32: { typedef typename MatType<MAT_T_INT32>::Type Type; { __VA_ARGS__; } break; } \
1297  case MAT_T_INT16: { typedef typename MatType<MAT_T_INT16>::Type Type; { __VA_ARGS__; } break; } \
1298  case MAT_T_INT8: { typedef typename MatType<MAT_T_INT8>::Type Type; { __VA_ARGS__; } break; } \
1299  case MAT_T_UINT64: { typedef typename MatType<MAT_T_UINT64>::Type Type; { __VA_ARGS__; } break; } \
1300  case MAT_T_UINT32: { typedef typename MatType<MAT_T_UINT32>::Type Type; { __VA_ARGS__; } break; } \
1301  case MAT_T_UINT16: { typedef typename MatType<MAT_T_UINT16>::Type Type; { __VA_ARGS__; } break; } \
1302  case MAT_T_UINT8: { typedef typename MatType<MAT_T_UINT8>::Type Type; { __VA_ARGS__; } break; } \
1303  default: \
1304  throw VarMatTypeError( streamstr("Uknown/unsupported encoded type from matio: " \
1305  << (matvar_ptr)->data_type) ); \
1306  } \
1307  } else { \
1308  switch ((matvar_ptr)->data_type) { \
1309  case MAT_T_DOUBLE: { typedef std::complex<typename MatType<MAT_T_DOUBLE>::Type> Type; { __VA_ARGS__; } break; } \
1310  case MAT_T_SINGLE: { typedef std::complex<typename MatType<MAT_T_SINGLE>::Type> Type; { __VA_ARGS__; } break; } \
1311  default: \
1312  throw VarMatTypeError( streamstr("Uknown/unsupported encoded type from matio: " \
1313  << (matvar_ptr)->data_type) ); \
1314  } \
1315  } \
1316  } while (false)
1317 
1318 
1319 
1320 
1321 // =============================================================================
1322 // now, let's define some ways to get variable values.
1323 // =============================================================================
1324 
1325 
1326 
1330 template<typename T>
1331 class TOMOGRAPHER_EXPORT VarValueDecoder<T,
1332 #ifdef TOMOGRAPHER_PARSED_BY_DOXYGEN
1333  _IS_NUMERIC_TYPE
1334 #else
1335  typename std::enable_if<(std::numeric_limits<T>::is_specialized ||
1336  tomo_internal::isComplex<T>::value)>::type
1337 #endif
1338  >
1339 {
1340 public:
1342  typedef T RetType;
1343 
1345  static inline void checkShape(const Var & var)
1346  {
1348  throw VarTypeError(var.varName(),
1349  streamstr("Can't store complex matrix in type " << typeid(T).name()));
1350  }
1351  if (var.numel() != 1) {
1352  throw VarTypeError(var.varName(),
1353  streamstr("Expected scalar but got matrix"));
1354  }
1355  }
1356 
1358  static inline RetType decodeValue(const Var & var)
1359  {
1360  const matvar_t * matvar_ptr = var.getMatvarPtr();
1361 
1362  MAT_SWITCH_TYPE(matvar_ptr,
1363  return get_value<Type>(matvar_ptr, var.varName());
1364  );
1365  // should never reach here
1366  tomographer_assert(false);
1367  }
1368 
1369 private:
1370  template<typename MATType,
1371  TOMOGRAPHER_ENABLED_IF_TMPL(!tomo_internal::isComplex<RetType>::value &&
1373  static inline RetType get_value(const matvar_t * matvar_ptr, const std::string & )
1374  {
1375  return RetType( ((const MATType *) matvar_ptr->data)[0] );
1376  }
1377 
1378  template<typename MATType,
1379  TOMOGRAPHER_ENABLED_IF_TMPL(tomo_internal::isComplex<RetType>::value &&
1381  static inline RetType get_value(const matvar_t * matvar_ptr, const std::string & )
1382  {
1383  return RetType( typename Tools::ComplexRealScalar<RetType>::type(((const MATType *) matvar_ptr->data)[0]),
1384  0 );
1385  }
1386 
1387  template<typename MATType,
1388  TOMOGRAPHER_ENABLED_IF_TMPL(!tomo_internal::isComplex<RetType>::value &&
1390  static inline RetType get_value(const matvar_t * /*matvar_ptr*/, const std::string & varname)
1391  {
1392  throw VarTypeError(varname, "Expected real scalar, got complex type");
1393  }
1394 
1395  template<typename MATType,
1396  TOMOGRAPHER_ENABLED_IF_TMPL(tomo_internal::isComplex<RetType>::value &&
1398  static inline RetType get_value(const matvar_t * matvar_ptr, const std::string & )
1399  {
1400  typedef typename Tools::ComplexRealScalar<MATType>::type MATRealType;
1401  const mat_complex_split_t * cdata = (mat_complex_split_t*) matvar_ptr->data;
1402 
1403  return RetType( typename Tools::ComplexRealScalar<RetType>::type( ((const MATRealType *) cdata->Re)[0] ),
1404  typename Tools::ComplexRealScalar<RetType>::type( ((const MATRealType *) cdata->Im)[0] ) );
1405  }
1406 
1407 };
1408 
1409 
1410 
1411 // -----------------------------------------------
1412 
1413 
1414 namespace tomo_internal {
1415 
1417 template<typename OutType, typename MatInnerT>
1418 class VarMatDataAccessor
1419 {
1420  const Var & p_var;
1421 
1422  typedef typename Tools::ComplexRealScalar<MatInnerT>::type MatRealInnerT;
1423 
1424  const MatInnerT * p_r_ptr;
1425  const MatRealInnerT * p_cre_ptr;
1426  const MatRealInnerT * p_cim_ptr;
1427 
1428 public:
1429  VarMatDataAccessor(const Var & var)
1430  : p_var(var)
1431  {
1432  const matvar_t * matvar_ptr = var.getMatvarPtr();
1433  tomographer_assert(matvar_ptr->data != NULL);
1434  if (!matvar_ptr->isComplex) {
1435  // real type
1436  p_r_ptr = (const MatInnerT*) matvar_ptr->data;
1437  p_cre_ptr = NULL;
1438  p_cim_ptr = NULL;
1439  } else {
1440  const mat_complex_split_t * cdata = (mat_complex_split_t*) matvar_ptr->data;
1441  p_r_ptr = NULL;
1442  p_cre_ptr = (const MatRealInnerT*) cdata->Re;
1443  p_cim_ptr = (const MatRealInnerT*) cdata->Im;
1444  tomographer_assert(p_cre_ptr != NULL);
1445  tomographer_assert(p_cim_ptr != NULL);
1446  }
1447  }
1448 
1449  virtual ~VarMatDataAccessor() { }
1450 
1451  template<typename IndexListType,
1452  typename OutType__ = OutType, typename MatInnerT__ = MatInnerT,
1453  TOMOGRAPHER_ENABLED_IF_TMPL(!tomo_internal::isComplex<OutType__>::value &&
1455  inline OutType value(IndexListType&& index) const
1456  {
1457  tomographer_assert(p_r_ptr != NULL);
1458 
1459  // real value.
1460  std::size_t lin = linear_index(std::forward<IndexListType>(index));
1461  return (OutType)p_r_ptr[lin];
1462  }
1463 
1464  template<typename IndexListType,
1465  typename OutType__ = OutType, typename MatInnerT__ = MatInnerT,
1466  TOMOGRAPHER_ENABLED_IF_TMPL(!tomo_internal::isComplex<OutType__>::value &&
1468  inline OutType value(IndexListType&& ) const
1469  {
1470  throw VarTypeError(p_var.varName(), "Expected real type, got complex");
1471  }
1472 
1473  template<typename IndexListType,
1474  typename OutType__ = OutType, typename MatInnerT__ = MatInnerT,
1475  TOMOGRAPHER_ENABLED_IF_TMPL(tomo_internal::isComplex<OutType__>::value &&
1477  inline OutType value(IndexListType&& index) const
1478  {
1479  tomographer_assert(p_r_ptr != NULL);
1480 
1481  // real value.
1482  std::size_t lin = linear_index(std::forward<IndexListType>(index));
1483  return OutType( typename Eigen::NumTraits<OutType>::Real(p_r_ptr[lin]) , 0 );
1484  }
1485 
1486  template<typename IndexListType,
1487  typename OutType__ = OutType, typename MatInnerT__ = MatInnerT,
1488  TOMOGRAPHER_ENABLED_IF_TMPL(tomo_internal::isComplex<OutType__>::value &&
1490  inline OutType value(IndexListType&& index) const
1491  {
1492  tomographer_assert(p_cre_ptr != NULL);
1493  tomographer_assert(p_cim_ptr != NULL);
1494 
1495  // complex value
1496  std::size_t lin = linear_index(std::forward<IndexListType>(index));
1497  return OutType( typename Eigen::NumTraits<OutType>::Real(p_cre_ptr[lin]),
1498  typename Eigen::NumTraits<OutType>::Real(p_cim_ptr[lin]) );
1499  }
1500 
1501 private:
1502  template<typename IndexListType>
1503  std::size_t linear_index(IndexListType && index) const
1504  {
1505  IndexList<false> ind_cmaj{
1506  p_var.dims(),
1508  };
1509  // std::cout << "linear_index: ind_cmaj=" << ind_cmaj << ", -> linear index = " << ind_cmaj.linearIndex() << "\n";
1510  return (std::size_t)ind_cmaj.linearIndex();
1511  }
1512 };
1513 
1514 } // namespace tomo_internal
1515 
1516 
1517 // -----------------------------------------------
1518 
1559 struct TOMOGRAPHER_EXPORT VarShape
1560 {
1568  const bool is_complex;
1569 
1584  const DimList dims;
1585 
1600  const bool is_square;
1601 
1606  template<typename DimListType>
1607  VarShape(bool is_complex_, DimListType&& dims_, bool is_square_)
1608  : is_complex(is_complex_), dims(std::forward<DimListType>(dims_)), is_square(is_square_)
1609  {
1610  _check_consistency();
1611  }
1612 
1618  VarShape(const Var & var)
1619  : is_complex(var.isComplex()), dims(var.dims()), is_square(var.isSquareMatrix())
1620  {
1621  }
1622 
1623 
1624 private:
1625  void _check_consistency()
1626  {
1627  if (is_square) {
1628  tomographer_assert(dims.size() == 0 || (dims.size() == 2 && (dims[0] == -1 || dims[1] == -1 || dims[0] == dims[1])));
1629  }
1630  }
1631 
1632 public:
1633 
1641  inline void checkShape(const VarShape & shape);
1642 
1650  inline void checkShape(const Var & var)
1651  {
1652  try {
1653  checkShape(VarShape(var));
1654  } catch (VarTypeError & err) {
1655  err.setVarName(var.varName());
1656  throw;
1657  }
1658  }
1659 };
1660 
1662 inline std::ostream& operator<<(std::ostream& str, const VarShape & varshape)
1663 {
1664  str << ((varshape.is_complex) ? std::string("complex ") : std::string("real "));
1665  if (varshape.dims.size() == 1) {
1666  str << "vector";
1667  } else if (varshape.dims.size() == 2) {
1668  if (varshape.is_square) {
1669  str << "square ";
1670  }
1671  str << "matrix";
1672  } else if (varshape.dims.size() > 0) {
1673  str << varshape.dims.size() << "-D array";
1674  if (varshape.is_square) {
1675  str << " with first two dimensions square";
1676  }
1677  } else {
1678  // empty dims means no requirement
1679  str << "array";
1680  }
1681  if (varshape.dims.size()) {
1682  str << " of shape " << varshape.dims;
1683  }
1684  return str;
1685 }
1686 
1687 inline void VarShape::checkShape(const VarShape& other)
1688 {
1689  const DimList mvardims = other.dims;
1690 
1691  if ((other.is_complex && !is_complex) ||
1692  (mvardims.size() != dims.size() && dims.size() > 0) ||
1693  (is_square && (mvardims.size() != 2 || mvardims[0] != mvardims[1])) ||
1694  (dims.size() > 0 && !mvardims.matchesWanted(dims))) {
1695 
1696  std::stringstream errstr;
1697 
1698  errstr << "Expected "
1699  << *this
1700  << ", got "
1701  << other;
1702 
1703  //fprintf(stderr, "Bad var type for variable %s\n", var.varName().c_str());
1704  throw VarTypeError(std::string(), errstr.str());
1705 
1706  }
1707 }
1708 
1709 
1710 
1711 // -----------------------------------------------------------------------------
1712 
1713 // get an std::vector<T> for a variable.
1714 
1724 template<typename T_, bool IsRowMajor_ = false>
1725 struct TOMOGRAPHER_EXPORT GetStdVector {
1726  typedef T_ type;
1727  static constexpr bool IsRowMajor = IsRowMajor_;
1728 };
1729 
1731 template<typename T, bool IsRowMajor>
1732 class TOMOGRAPHER_EXPORT VarValueDecoder<GetStdVector<T, IsRowMajor> >
1733 {
1734 public:
1735  typedef std::vector<T> RetType;
1736 
1737  static inline void checkShape(const Var & var)
1738  {
1740  throw VarTypeError(var.varName(),
1741  std::string("can't store complex matrix in type ")
1742  + typeid(T).name());
1743  }
1744  }
1745 
1746  static inline std::vector<T> decodeValue(const Var & var)
1747  {
1748  const matvar_t * matvar_ptr = var.getMatvarPtr();
1749  std::size_t numel = (std::size_t)var.numel();
1750  std::vector<T> val(numel);
1751 
1752  MAT_SWITCH_TYPE(matvar_ptr,
1753  tomo_internal::VarMatDataAccessor<T, Type> acc(var);
1754 
1755  for(IndexListIterator<IsRowMajor> il(var.dims()); il.valid(); ++il) {
1756  // std::cout << "index: " << il << "\n";
1757  val[(std::size_t)il.linearIndex()] = acc.value(il);
1758  }
1759 
1760  return val;
1761  );
1762  // should never reach here
1763  tomographer_assert(false);
1764  }
1765 
1766 };
1767 
1768 
1769 
1770 
1771 
1772 
1773 
1774 // =============================================================================
1775 
1776 // get Eigen types.
1777 
1778 namespace tomo_internal {
1779 
1780 
1781 inline DimList dims_stackedcols(DimList vdims)
1782 {
1783  tomographer_assert(vdims.size() >= 1);
1784  DimList vdimsreshaped;
1785  if (vdims.size() == 1) {
1786  vdimsreshaped = vdims;
1787  vdimsreshaped << 1;
1788  } else if (vdims.size() == 2) {
1789  vdimsreshaped = vdims;
1790  } else if (vdims.size() > 2) {
1791  vdimsreshaped << getNumEl(vdims.data(), vdims.data()+vdims.size()-1) << vdims[vdims.size()-1];
1792  }
1793  tomographer_assert(vdimsreshaped[0] != -1 && vdimsreshaped[1] != -1);
1794  return vdimsreshaped;
1795 }
1796 
1797 template<typename MatrixType, typename MatType,
1798  TOMOGRAPHER_ENABLED_IF_TMPL(Eigen::NumTraits<typename MatrixType::Scalar>::IsComplex &&
1800 inline void init_eigen_matrix(MatrixType & matrix, const DimList & vdims,
1801  const Var & var, const std::ptrdiff_t data_offset = 0)
1802 {
1803  typedef typename MatrixType::Scalar Scalar;
1804  typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
1805 
1806  typedef typename Eigen::NumTraits<MatType>::Real MatRealType;
1807 
1808  const matvar_t * matvar_ptr = var.getMatvarPtr();
1809 
1810  const mat_complex_split_t * cdata = (mat_complex_split_t*) matvar_ptr->data;
1811 
1812  DimList vdimsreshaped = dims_stackedcols(vdims);
1813 
1814  matrix = (
1816  (const MatRealType *) cdata->Re + data_offset, vdimsreshaped[0], vdimsreshaped[1]
1817  ).template cast<std::complex<RealScalar> >()
1818  + std::complex<RealScalar>(0,1) *
1820  (const MatRealType *) cdata->Im + data_offset, vdimsreshaped[0], vdimsreshaped[1]
1821  ).template cast<std::complex<RealScalar> >()
1822  );
1823 }
1824 
1825 template<typename MatrixType, typename MatRealType,
1826  TOMOGRAPHER_ENABLED_IF_TMPL(!Eigen::NumTraits<typename MatrixType::Scalar>::IsComplex &&
1828 inline void init_eigen_matrix(MatrixType & /*matrix*/, const DimList & /*vdims*/,
1829  const Var & var, const std::ptrdiff_t /*data_offset*/ = 0)
1830 {
1831  throw VarTypeError(var.varName(), "Expected real type, but got complex.");
1832 }
1833 
1834 template<typename MatrixType, typename MatRealType,
1835  TOMOGRAPHER_ENABLED_IF_TMPL(!Eigen::NumTraits<MatRealType>::IsComplex)>
1836 inline void init_eigen_matrix(MatrixType & matrix, const DimList & vdims,
1837  const Var & var, const std::ptrdiff_t data_offset = 0)
1838 {
1839  typedef typename MatrixType::Scalar Scalar;
1840 
1841  const matvar_t * matvar_ptr = var.getMatvarPtr();
1842 
1843  DimList vdimsreshaped = dims_stackedcols(vdims);
1844 
1845  matrix = (
1847  (const MatRealType *) matvar_ptr->data + data_offset, vdimsreshaped[0], vdimsreshaped[1]
1848  )
1849  ).template cast<Scalar>();
1850 }
1851 
1852 } // namespace tomo_internal
1853 
1854 
1855 
1864 template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
1865 class TOMOGRAPHER_EXPORT VarValueDecoder<Eigen::Matrix<Scalar,Rows,Cols,Options,MaxRows,MaxCols> >
1866 {
1867 public:
1868 
1870  typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
1871  typedef MatrixType RetType;
1872 
1873  /* not needed in fact.
1874  struct Params {
1875  Params() : dims{Rows, Cols}
1876  {
1877  }
1878  Params(int rows) : dims{rows, 1}
1879  {
1880  tomographer_assert(Cols == 1);
1881  tomographer_assert(Rows == Eigen::Dynamic || Rows == rows);
1882  }
1883  Params(int rows, int cols) : dims{rows, cols}
1884  {
1885  tomographer_assert(Rows == Eigen::Dynamic || Rows == rows);
1886  tomographer_assert(Cols == Eigen::Dynamic || Cols == cols);
1887  }
1888  Params(const DimList& dims_) : dims{dims_}
1889  {
1890  tomographer_assert(dims.size() == 2);
1891  tomographer_assert(Rows == Eigen::Dynamic || Rows == dims[0]);
1892  tomographer_assert(Cols == Eigen::Dynamic || Cols == dims[1]);
1893  }
1894 
1895  const DimList dims;
1896  };
1897  */
1898 
1899  static inline void checkShape(const Var & var) //, const Params & p = Params())
1900  {
1901  DimList matdims;
1902  matdims << (Rows!=Eigen::Dynamic ? Rows : -1)
1903  << (Cols!=Eigen::Dynamic ? Cols : -1);
1904 
1906  matdims,
1907  matdims[0] != -1 && matdims[0] == matdims[1]);
1908 
1909  DimList vdims = var.dims();
1910  if (vdims.size() > 2) {
1911  vdims = (DimList() << getNumEl(vdims.data(), vdims.data()+vdims.size()-1) << vdims[vdims.size()-1]);
1912  }
1913  try {
1914  shape.checkShape(VarShape(var.isComplex(), vdims, (vdims[0] != -1 && vdims[0] == vdims[1])));
1915  } catch (VarTypeError & err) {
1916  err.setVarName(var.varName());
1917  throw; // re-throw
1918  }
1919  }
1920 
1921  static inline RetType decodeValue(const Var & var) //, const Params & p = Params())
1922  {
1923  const matvar_t * matvar_ptr = var.getMatvarPtr();
1924  DimList vdims{var.dims()};
1925 
1926  if (vdims.size() < 2) {
1927  throw VarTypeError(var.varName(), streamstr("Expected matrix, but variable shape is " << vdims));
1928  }
1929  // if we are stacking several dimensions in the column, force column-major ordering.
1930  if (vdims.size() > 2 &&
1931  ((MatrixType::Options & Eigen::RowMajorBit) == Eigen::RowMajor)) {
1932  throw VarTypeError(var.varName(),
1933  "When collapsing several dimensions into Eigen columns, you must use "
1934  "column-major ordering (sorry).");
1935  }
1936 
1937  const int cols = vdims[vdims.size()-1];
1938  // rest of dimensions.
1939  const int rows = getNumEl(vdims.data(), vdims.data()+vdims.size()-1);
1940 
1941  MatrixType matrix(rows, cols);
1942 
1944  matvar_ptr,
1945  tomo_internal::init_eigen_matrix<MatrixType, Type>(
1946  matrix, vdims, var
1947  );
1948  );
1949 
1950  return matrix;
1951  }
1952 
1953 };
1954 
1955 
1956 
1957 
1958 
1968 template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols, typename Alloc>
1969 class TOMOGRAPHER_EXPORT VarValueDecoder<std::vector<Eigen::Matrix<Scalar,Rows,Cols,Options,MaxRows,MaxCols>, Alloc> >
1970 {
1971 public:
1972 
1974  typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
1975 
1977 
1978  static inline void checkShape(const Var & var)
1979  {
1980  DimList matdims;
1981  matdims << (Rows!=Eigen::Dynamic ? Rows : -1)
1982  << (Cols!=Eigen::Dynamic ? Cols : -1)
1983  << -1;
1984 
1986  matdims,
1987  false);
1988 
1989  DimList vdims = var.dims();
1990  if (vdims.size() < 1) {
1991  throw VarTypeError(var.varName(), "Invalid (empty) variable dimensions");
1992  }
1993  if (vdims.size() == 1) {
1994  vdims << 1 << 1;
1995  } else if (vdims.size() == 2) {
1996  vdims << 1;
1997  } else { //if (vdims.size() > 2) {
1998  vdims = (DimList() << vdims[0] << vdims[1] << getNumEl(vdims.data()+2, vdims.data()+vdims.size()));
1999  }
2000 
2001  // check shape now:
2002  try {
2003  shape.checkShape(VarShape(var.isComplex(), vdims, false));
2004  } catch (VarTypeError & err) {
2005 
2006  err.setVarName(var.varName());
2007  throw; // re-throw
2008  }
2009 
2010  }
2011 
2012  static inline RetType decodeValue(const Var & var)
2013  {
2014  DimList vardims{var.dims()};
2015 
2016  tomographer_assert(vardims.size() >= 1);
2017 
2018  DimList innerdims;
2019  std::size_t outerdim = 1;
2020  if (vardims.size() == 1) {
2021  innerdims << vardims[0] << 1;
2022  } else if (vardims.size() == 2) {
2023  innerdims = vardims;
2024  } else {
2025  // more dimensions
2026  innerdims << vardims[0] << vardims[1];
2027  outerdim = (std::size_t)getNumEl(vardims.data()+2, vardims.data()+vardims.size());
2028  }
2029 
2030  RetType value(outerdim);
2031 
2032  std::ptrdiff_t innernumel = innerdims[0]*innerdims[1];
2033 
2034  std::size_t j;
2035  const matvar_t * matvar_ptr = var.getMatvarPtr();
2036 
2038  matvar_ptr,
2039  for (j = 0; j < outerdim; ++j) {
2040  tomo_internal::init_eigen_matrix<MatrixType,Type>(
2041  value[j], // matrix reference
2042  innerdims, // dimensions of matrix reference
2043  var, // data
2044  (std::ptrdiff_t)j*innernumel // offset
2045  );
2046  }
2047  );
2048 
2049  return value;
2050  }
2051 
2052 };
2053 
2054 
2055 // get a (guaranteed) positive semidefinite matrix from the variable, along with its
2056 // matrix square root
2057 
2058 template<typename EigenType> class EigenPosSemidefMatrixWithSqrt;
2059 template<typename EigenType> class VarValueDecoder<EigenPosSemidefMatrixWithSqrt<EigenType> >;
2060 
2061 
2077 template<typename EigenType_>
2078 class TOMOGRAPHER_EXPORT EigenPosSemidefMatrixWithSqrt
2079 {
2080 public:
2081  typedef EigenType_ EigenType;
2082  typedef typename Eigen::NumTraits<typename EigenType::Scalar>::Real RealScalarType;
2083 
2096  struct Params {
2099  : tolerance(tolerance_) { }
2101  const RealScalarType tolerance;
2102  };
2103 
2106 
2108  EigenType mat;
2110  EigenType sqrt;
2111 
2112 private:
2113  friend class VarValueDecoder<EigenPosSemidefMatrixWithSqrt<EigenType> >;
2114 
2115  EigenPosSemidefMatrixWithSqrt(EigenType && m, const Params & p)
2116  : eig( m )
2117  {
2118  typedef typename Eigen::SelfAdjointEigenSolver<EigenType>::RealVectorType RealVectorType;
2119 
2120  auto U = eig.eigenvectors();
2121  auto d = eig.eigenvalues();
2122 
2123  Tomographer::MathTools::forcePosVecKeepSum<RealVectorType>(d, p.tolerance);
2124 
2125  mat = U * d.asDiagonal() * U.adjoint();
2126  sqrt = U * d.cwiseSqrt().asDiagonal() * U.adjoint();
2127  }
2128 };
2130 template<typename EigenType_>
2131 class TOMOGRAPHER_EXPORT VarValueDecoder<EigenPosSemidefMatrixWithSqrt<EigenType_> >
2132 {
2133 public:
2134  typedef EigenType_ EigenType;
2136 
2138 
2140 
2141  static inline void checkShape(const Var & var, const Params & = Params())
2142  {
2143  // make sure we can decode the given variable
2144  EigenVarDecoder::checkShape(var);
2145  // make sure that var is a square matrix
2146  const DimList d = var.dims();
2147  if (d.size() != 2 || d[0] != d[1]) {
2148  throw VarTypeError(var.varName(), std::string("variable to decode is not a square matrix"));
2149  }
2150  }
2151 
2152  static inline RetType decodeValue(const Var & var, const Params & p = Params())
2153  {
2154  return EigenPosSemidefMatrixWithSqrt<EigenType>(EigenVarDecoder::decodeValue(var), p);
2155  }
2156 
2157 };
2158 
2159 
2160 
2161 
2162 
2163 } // namespace MAT
2164 } // namespace Tomographer
2165 
2166 
2167 
2168 
2169 #endif
Utilities for formatting strings.
Utility to iterate over a multidim array by increasing linear index.
Definition: ezmatio.h:627
File(File &&other)
Move constructor.
Definition: ezmatio.h:215
DimList()
Initializes an empty dimension list.
Definition: ezmatio.h:310
DimList dims() const
Specific dimensions of this numeric array or tensor.
Definition: ezmatio.h:1027
const std::string & varName() const
The variable name.
Definition: ezmatio.h:1008
Base namespace for the Tomographer project.
Definition: densellh.h:45
IndexList(const std::vector< int > &dims=std::vector< int >(), int linearindex=-1)
Constructor with linear index.
Definition: ezmatio.h:453
Var var(const std::string varname, bool load_data=true)
Find a variable by name.
Definition: ezmatio.h:1170
Invalid index or index list provided to a routine.
Definition: ezmatio.h:168
const RealVectorType & eigenvalues() const
Unknown type of a variable present in the data file.
Definition: ezmatio.h:135
void checkShape(const Var &var)
Shorthand for checkShape(const VarShape& other) for a shape given by a variable.
Definition: ezmatio.h:1650
Error while reading a variable from the MATLAB data file.
Definition: ezmatio.h:113
const std::vector< int > & index() const
Return a reference to *this.
Definition: ezmatio.h:547
DimList(VectorType &&dims)
Initialize a dimension list with another list of dimensions.
Definition: ezmatio.h:318
int numel() const
Get the total number of elements in an array of these dimensions.
Definition: ezmatio.h:352
A MATLAB file open for reading data.
Definition: ezmatio.h:192
EigenType mat
The matrix itself.
Definition: ezmatio.h:2108
int ndims() const
Get the number of dimensions in this array.
Definition: ezmatio.h:359
STL namespace.
STL class.
Var & operator=(Var &&other)
Move assignment operator as this object implements C++11 move semantics.
Definition: ezmatio.h:1150
T end(T... args)
File(const std::string fname)
Open a data file for reading.
Definition: ezmatio.h:200
VarValueDecoder< T >::RetType value() const
Read this variable data as native C++ object.
Definition: ezmatio.h:1116
A list of indices with an API for linear or subindices access.
Definition: ezmatio.h:440
Base Exception class for errors within our MAT routines.
Definition: ezmatio.h:70
Specializable template which takes care of decoding values.
Definition: ezmatio.h:802
IndexList & operator<<(int ind)
Append index to list.
Definition: ezmatio.h:568
IndexList(const std::vector< int > &dims, VectorIntInitializer &&index)
Constructor with multidimensional index specification.
Definition: ezmatio.h:473
std::vector< Var > getVarInfoList()
Get a list of all variables in this data file.
Definition: ezmatio.h:1175
STL class.
Map matio&#39;s constants to C/C++ types.
Definition: ezmatio.h:1218
Tools for dealing with positive semidefinite matrices.
Error while opening a MATLAB file.
Definition: ezmatio.h:149
T push_back(T... args)
static Var takeOver(matvar_t *varinfo)
Take in charge the given C matvar_t pointer.
Definition: ezmatio.h:999
T data(T... args)
Var(Var &&other)
Var implements C++11 move semantics.
Definition: ezmatio.h:976
Basic utilities for dealing with Eigen matrices and other types.
void checkShape(const VarShape &shape)
Verify that our requirements match the given other shape.
Definition: ezmatio.h:1687
const bool is_complex
Whether the variable is or should be complex.
Definition: ezmatio.h:1568
A MATLAB variable in the MAT file.
Definition: ezmatio.h:918
std::vector< int > VectorType
Base vector type (superclass)
Definition: ezmatio.h:446
EigenType sqrt
The matrix square root of the read matrix.
Definition: ezmatio.h:2110
Params(RealScalarType tolerance_=Eigen::NumTraits< RealScalarType >::dummy_precision())
Constructor, specify the tolerance.
Definition: ezmatio.h:2098
const RealScalarType tolerance
The tolerance.
Definition: ezmatio.h:2101
VarValueDecoder< T >::RetType value(const Var &var, const typename VarValueDecoder< T >::Params &params)
Access the value of the given variable, as a C++ type.
Definition: ezmatio.h:893
T str(T... args)
Ask for this type in Var::value<typename T>() to get an std::vector of the data.
Definition: ezmatio.h:1725
Parameters to read the positive semidefinite matrix.
Definition: ezmatio.h:2096
int ndims() const
Number of dimensions of this object.
Definition: ezmatio.h:1018
Some C++ utilities, with a tad of C++11 tricks.
STL class.
T move(T... args)
Read a positive semidefinite matrix along with its matrix square root.
Definition: ezmatio.h:2058
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
static RetType decodeValue(const Var &var)
Decode the variable var into the C++ type T.
Definition: ezmatio.h:1191
T size(T... args)
static void checkShape(const Var &var)
Check that the type and shape of the var are compatible with T.
Definition: ezmatio.h:818
std::ostream & operator<<(std::ostream &out, const DimList &dlist)
C++ output stream operators for a DimList.
Definition: ezmatio.h:398
#define MAT_SWITCH_TYPE(matvar_ptr,...)
Useful hack to get C++ type from dynamical MAT type ID.
Definition: ezmatio.h:1289
STL class.
An array of ints which specifies a list of dimensions.
Definition: ezmatio.h:306
Var & operator=(const Var &copy)
Var objects are copyable. Beware though that the data is shared.
Definition: ezmatio.h:1158
const std::vector< int > & dims() const
Get the underlying dimensions given to the constructor.
Definition: ezmatio.h:556
VarShape(bool is_complex_, DimListType &&dims_, bool is_square_)
Construct a VarShape object from given characteristics.
Definition: ezmatio.h:1607
bool isComplex() const
Whether this variable is complex or real.
Definition: ezmatio.h:1050
File & operator=(File &&other)
Move assignment operator.
Definition: ezmatio.h:263
T begin(T... args)
VarShape(const Var &var)
Construct a VarShape object reflecting the actual shape of a variable in the data file...
Definition: ezmatio.h:1618
bool isSquareMatrix() const
Whether this is a square matrix.
Definition: ezmatio.h:1059
const MatrixType & eigenvectors() const
::value DimList(std::initializer_list< T > init)
Initialize the dimension list with an explicit initializer.
Definition: ezmatio.h:334
VarValueDecoder< T >::RetType value(const Var &var)
Access the value of the given variable, as a C++ type.
Definition: ezmatio.h:878
T c_str(T... args)
const matvar_t * getMatvarPtr() const
Access the underlying C pointer to the MatIO structure. Use with care.
Definition: ezmatio.h:1141
static void checkShape(const Var &var)
See VarValueDecoder::checkShape()
Definition: ezmatio.h:1345
Type mismatch (wrong type requested) in a variable read from the MATLAB data file.
Definition: ezmatio.h:124
Exception relating to a MATLAB variable in the data file.
Definition: ezmatio.h:99
Describe shape of variable and whether it is complex.
Definition: ezmatio.h:1559
mat_t * getMatPtr()
Direct access to the underlying C pointer used for the MatIO library.
Definition: ezmatio.h:253
DimList & operator<<(int dim)
Add (append) a dimension to the dimension list.
Definition: ezmatio.h:382
static RetType decodeValue(const Var &var)
See VarValueDecoder::decodeValue()
Definition: ezmatio.h:1358
const unsigned int RowMajorBit
T RetType
Type returned by decodeValue()
Definition: ezmatio.h:810
#define streamstr(tokens)
Utility macro to format stream tokens to a std::string.
Definition: fmt.h:149
Var(File &matf, const std::string &varname, bool load_data=true)
Read variable from MATLAB data file.
Definition: ezmatio.h:951
bool matchesWanted(const DimList &wanted) const
Test whether our dimensions match the given required list of dimensions.
Definition: ezmatio.h:364
Var(const Var &copy)
Var objects are copyable. Beware though that the data is shared.
Definition: ezmatio.h:969
DimList(It b, It e)
Initialize the dimension list with an iterator range.
Definition: ezmatio.h:344
STL class.
const DimList dims
The dimensions of the variable, or requested dimensions.
Definition: ezmatio.h:1584
#define tomographer_assert(...)
Assertion test macro.
Definition: cxxdefs.h:84
RowMajor
ValueType getNumEl(It begin, It end)
Calculate the product of all dimensions.
Definition: ezmatio.h:283
const bool is_square
Whether the variable&#39;s two first dimensions are (or should be) the same.
Definition: ezmatio.h:1600
VarValueDecoder< T >::RetType value(const typename VarValueDecoder< T >::Params &params)
Read this variable data as native C++ object.
Definition: ezmatio.h:1131
Eigen::SelfAdjointEigenSolver< EigenType > eig
The readily-initialized eigenvalue decomposition of the matrix read from the MAT file.
Definition: ezmatio.h:2105
bool hasData() const
Whether data for this Var object has been loaded.
Definition: ezmatio.h:1070
int numel() const
The total number of elements in this array or tensor.
Definition: ezmatio.h:1037