1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // Antioch - A Gas Dynamics Thermochemistry Library
5 //
6 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
7 // Sylvain Plessis, Roy H. Stonger
8 //
9 // Copyright (C) 2013 The PECOS Development Team
10 //
11 // This library is free software; you can redistribute it and/or
12 // modify it under the terms of the Version 2.1 GNU Lesser General
13 // Public License as published by the Free Software Foundation.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
23 // Boston, MA 02110-1301 USA
24 //
25 //-----------------------------------------------------------------------el-
26
27 #ifndef ANTIOCH_EIGEN_UTILS_H
28 #define ANTIOCH_EIGEN_UTILS_H
29
30 #ifdef ANTIOCH_METAPROGRAMMING_H
31 # ifndef ANTIOCH_EIGEN_UTILS_DECL_H
32 # error eigen_utils_decl.h must be included before metaprogramming.h
33 # endif
34 #endif
35
36 #include "antioch_config.h"
37
38 #ifdef ANTIOCH_HAVE_EIGEN
39 // Though the following implementations are all valid without <Eigen/Dense>,
40 // successfully using them with Eigen types requires Eigen be included first.
41 // Configure-time Eigen support enforces this constraint but header-only Eigen
42 // may be mixed with header-only Antioch without configure-time flags.
43 #include <Eigen/Dense>
44 #endif
45
46 #include "antioch/metaprogramming.h"
47
48 // Notice _Matrix template templates might be Eigen::Matrix or Eigen::Array.
49 // Therefore, always use .array()- or .matrix()-like operations for robustness.
50 // Otherwise Eigen will complain with YOU_CANNOT_MIX_ARRAYS_AND_MATRICES.
51
52 namespace std
53 {
54
55 template <
56 template <typename, int, int, int, int, int> class _Matrix,
57 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
58 >
59 inline
60 ANTIOCH_AUTO(_Matrix<_Scalar ANTIOCH_COMMA _Rows ANTIOCH_COMMA _Cols ANTIOCH_COMMA _Options ANTIOCH_COMMA _MaxRows ANTIOCH_COMMA _MaxCols>)
61 max(const _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & a,
62 const _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & b)
63 ANTIOCH_AUTOFUNC(_Matrix<_Scalar ANTIOCH_COMMA _Rows ANTIOCH_COMMA _Cols ANTIOCH_COMMA _Options ANTIOCH_COMMA _MaxRows ANTIOCH_COMMA _MaxCols>,
64 a.array().max(b.array()))
65
66 template <
67 template <typename, int, int, int, int, int> class _Matrix,
68 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
69 >
70 inline
71 ANTIOCH_AUTO(_Matrix<_Scalar ANTIOCH_COMMA _Rows ANTIOCH_COMMA _Cols ANTIOCH_COMMA _Options ANTIOCH_COMMA _MaxRows ANTIOCH_COMMA _MaxCols>)
72 min(const _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & a,
73 const _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> & b)
74 ANTIOCH_AUTOFUNC(_Matrix<_Scalar ANTIOCH_COMMA _Rows ANTIOCH_COMMA _Cols ANTIOCH_COMMA _Options ANTIOCH_COMMA _MaxRows ANTIOCH_COMMA _MaxCols>,
75 a.array().min(b.array()))
76
77 } // end namespace std
78
79
80 namespace Antioch
81 {
82
83 template <typename T>
84 inline
85 typename Antioch::enable_if_c<is_eigen<T>::value,
86 typename value_type<T>::type
87 >::type
max(const T & in)88 max(const T& in)
89 {
90 return in.maxCoeff();
91 }
92
93 template <typename T>
94 inline
95 typename Antioch::enable_if_c<is_eigen<T>::value,
96 typename value_type<T>::type
97 >::type
min(const T & in)98 min(const T& in)
99 {
100 return in.minCoeff();
101 }
102
103 template <typename T>
104 inline
105 typename Antioch::enable_if_c<is_eigen<T>::value,
106 bool>::type
has_nan(const T & in)107 has_nan(const T& in)
108 {
109 using std::isnan;
110
111 const size_t n_rows = in.rows();
112 const size_t n_cols = in.cols();
113 for (size_t i=0; i != n_rows; ++i)
114 for (size_t j=0; j != n_cols; ++j)
115 if(isnan(in(i,j)))
116 return true;
117 return false;
118 }
119
120 template <typename T>
121 struct has_size<T, typename Antioch::enable_if_c<is_eigen<T>::value,void>::type>
122 {
123 const static bool value = true;
124 };
125
126 template <typename T>
127 struct return_auto<T, typename Antioch::enable_if_c<is_eigen<T>::value,void>::type>
128 {
129 const static bool value = true;
130 };
131
132 template <typename T>
133 struct size_type<T, typename Antioch::enable_if_c<is_eigen<T>::value,void>::type>
134 {
135 typedef typename T::Index type;
136 };
137
138 template <
139 template <typename, int, int, int, int, int> class _Matrix,
140 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
141 >
142 struct value_type<_Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
143 {
144 typedef _Scalar type;
145 };
146
147 template <
148 template <typename, int, int, int, int, int> class _Matrix,
149 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
150 >
151 struct raw_value_type<_Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
152 {
153 typedef typename raw_value_type<_Scalar>::type type;
154 };
155
156 template <
157 template <typename, int, int, int, int, int> class _Matrix,
158 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
159 >
160 inline
161 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
162 zero_clone(const _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& ex)
163 {
164 // We can't just use setZero here with arbitrary _Scalar types
165 if (ex.size())
166 return
167 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
168 (ex.rows(), ex.cols()).setConstant(zero_clone(ex[0]));
169
170 return
171 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
172 (ex.rows(), ex.cols());
173 }
174
175 template <
176 template <typename, int, int, int, int, int> class _Matrix,
177 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols,
178 typename Scalar2
179 >
180 inline
181 void
182 zero_clone(_Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& output,
183 const _Matrix<Scalar2, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& ex)
184 {
185 // We can't just use setZero here with arbitrary _Scalar types
186 output.resize(ex.rows(), ex.cols());
187 output.setConstant(zero_clone(ex[0]));
188 }
189
190
191 template <
192 template <typename, int, int, int, int, int> class _Matrix,
193 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols,
194 typename Scalar
195 >
196 inline
197 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
198 constant_clone(const _Matrix<_Scalar, _Rows, _Cols, _Options,
199 _MaxRows, _MaxCols>& ex,
200 const Scalar& value)
201 {
202 // We can't just use setZero here with arbitrary _Scalar types
203 if (ex.size())
204 return
205 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
206 (ex.rows(), ex.cols()).setConstant(value);
207
208 return
209 _Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
210 (ex.rows(), ex.cols());
211 }
212
213
214 template <
215 template <typename, int, int, int, int, int> class _Matrix,
216 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols,
217 typename Scalar
218 >
219 inline
220 void
221 constant_fill(_Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows,
222 _MaxCols>& output,
223 const Scalar& value)
224 {
225 output.fill(value);
226 }
227
228 template <
229 template <typename, int, int, int, int, int> class _Matrix,
230 typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
231 >
232 inline
233 void
234 set_zero(_Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& a)
235 {
236 // We can't just use setZero here with arbitrary value_types
237 if (a.size())
238 a.setConstant (zero_clone(a[0]));
239 }
240
241 /*
242 template <template <typename, int, int, int, int, int> class _Matrix,
243 typename T, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols,
244 typename VectorScalar
245 >
246 inline
247 _Matrix<T,_Rows,_Cols,_Options,_MaxRows,_MaxCols>
248 custom_clone(const _Matrix<T,_Rows,_Cols,_Options,_MaxRows,_MaxCols> & /example/, const VectorScalar& values, const _Matrix<unsigned int,_Rows,_Cols,_Options,_MaxRows,_MaxCols> & indexes)
249 {
250 _Matrix<T,_Rows,_Cols,_Options,_MaxRows,_MaxCols> returnval;
251 returnval.resize(indexes.size());
252 for(std::size_t i = 0; i < indexes.size(); i++)
253 {
254 returnval[i] = values[indexes[i]];
255 }
256 return returnval;
257 }
258 */
259
260 template <typename Condition, typename T1, typename T2>
261 inline
262 typename enable_if_c<
263 is_eigen<T1>::value &&
264 is_eigen<T2>::value,
265 typename state_type<T1>::type
266 >::type
267 if_else(
268 const Condition& condition,
269 const T1& if_true,
270 const T2& if_false)
271 {
272 return condition.select(if_true, if_false);
273 }
274
275
276 template <typename VectorT,
277 template <typename , int, int, int, int, int> class _Matrix,
278 typename _UIntT, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols
279 >
280 inline
281 typename enable_if_c<
282 is_eigen<typename value_type<VectorT>::type>::value,
283 typename value_type<VectorT>::type
284 >::type
285 eval_index(const VectorT& vec, const _Matrix<_UIntT, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& index)
286 {
287 typename value_type<VectorT>::type returnval = vec[0];
288 for (unsigned int i=0; i != index.size(); ++i)
289 returnval[i] = vec[index[i]][i];
290 return returnval;
291 }
292
293 template <typename T>
294 inline
295 bool conjunction_root(const T & vec, eigen_library_tag)
296 {
297 return vec.all();
298 }
299
300 template <typename T>
301 inline
302 bool disjunction_root(const T & vec, eigen_library_tag)
303 {
304 return vec.any();
305 }
306
307
308 } // end namespace Antioch
309
310 #endif //ANTIOCH_EIGEN_UTILS_H
311