Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "Moose.h"
13 : #include "MooseError.h"
14 : #include "DataIO.h"
15 : #include <algorithm>
16 :
17 : /**
18 : * Implements a container class for multi-indexed objects
19 : * with an arbitrary number of indices.
20 : */
21 : template <class T>
22 : class MultiIndex
23 : {
24 : public:
25 : /// MultiIndex container iterator
26 : template <bool is_const = false>
27 : class const_noconst_iterator;
28 :
29 : ///@{ container related types and categories
30 : typedef T value_type;
31 : typedef std::vector<unsigned int> size_type;
32 : using iterator = const_noconst_iterator<false>;
33 : using const_iterator = const_noconst_iterator<true>;
34 : ///@}
35 :
36 : /// construct zero initialized container of a given shape
37 : MultiIndex(const size_type & shape);
38 :
39 : /// construct container of a given shape initialized from a flat data blob
40 : MultiIndex(const size_type & shape, const std::vector<T> & data);
41 :
42 : ///@{ element access operators
43 : T & operator()(const size_type & indices);
44 : const T & operator()(const size_type & indices) const;
45 : ///@}
46 :
47 : ///@{ direct data access via bracket operator
48 8 : T & operator[](unsigned int j) { return _data[j]; }
49 21305 : const T & operator[](unsigned int j) const { return _data[j]; }
50 1 : T & at(unsigned int j) { return _data.at(j); }
51 : ///@}
52 :
53 : ///@{ accesses a slice of the multi index object
54 : MultiIndex<T> slice(unsigned int dimension, unsigned int index) const;
55 : MultiIndex<T> slice(size_type dimension, size_type index) const;
56 : ///@}
57 :
58 : /// container size as dim dimensional vector
59 5361 : const size_type & size() const { return _shape; }
60 :
61 : /// dimension of the container
62 32301 : unsigned int dim() const { return _shape.size(); }
63 :
64 : /// total number of values stored in the container
65 1 : unsigned int nEntries() const { return _nentries; }
66 :
67 : /// get the raw data vector
68 : std::vector<T> getRawData() const { return _data; }
69 :
70 : /// Resize container. Must keep dimensionality constant.
71 : void resize(const size_type & shape);
72 :
73 : /// Reshape container arbitrarily and initialize with value
74 : void assign(const size_type & shape, T value);
75 :
76 : ///@{ access the stride
77 : unsigned int stride(unsigned int j) const;
78 2664 : const size_type & stride() const { return _stride; }
79 : ///@}
80 :
81 : ///@{ Implement loadHelper and storeHelper for easier data (de)serialization
82 : void dataStore(std::ostream & stream, void * context);
83 : void dataLoad(std::istream & stream, void * context);
84 : ///@}
85 :
86 : ///@{ iterators for begin and end of this container
87 2678 : iterator begin() { return const_noconst_iterator<false>(*this, 0); }
88 2668 : iterator end() { return const_noconst_iterator<false>(*this, _nentries); }
89 : const_iterator begin() const { return const_noconst_iterator<true>(*this, 0); }
90 : const_iterator end() const { return const_noconst_iterator<true>(*this, _nentries); }
91 : ///@}
92 :
93 : /// compute the flat index for a size_type index
94 : unsigned int flatIndex(const size_type & indices) const;
95 :
96 : protected:
97 : /// given a flat index computes the vector of indices i0, i1, i2, ...
98 : void findIndexVector(unsigned int flat_index, size_type & indices) const;
99 :
100 : /// the size along each index
101 : size_type _shape;
102 :
103 : /// the number of dimensions TODO: get rid of this -> _shape.size()
104 : unsigned int _dim;
105 :
106 : /// the number of entries TODO: get rid of this -> _data.size()
107 : unsigned int _nentries;
108 :
109 : /// stride for each index, e.g. if you know {i, j, k} -> flat_index, {i, j + 1, k} = flat_index + stride[1]
110 : size_type _stride;
111 :
112 : /// the data unrolled into a vector
113 : std::vector<T> _data;
114 :
115 : private:
116 : /// build accumulated shape vector for flat index calculation
117 : void buildAccumulatedShape();
118 :
119 : /// change the container shape and reset meta data
120 : void reshape(const size_type & shape);
121 : };
122 :
123 : /**
124 : * Nested iterator class for MultiIndex containers
125 : */
126 : template <class T>
127 : template <bool is_const>
128 : class MultiIndex<T>::const_noconst_iterator
129 : {
130 : public:
131 : typedef typename std::conditional<is_const, const MultiIndex<T> &, MultiIndex<T> &>::type
132 : reference_type;
133 :
134 5346 : const_noconst_iterator(reference_type multi_index, unsigned int position)
135 5346 : : _multi_index(multi_index), _flat_index(position), _shape(multi_index.size())
136 : {
137 5346 : _multi_index.findIndexVector(position, _indices);
138 5346 : }
139 :
140 : // Simple data getters
141 24457 : unsigned int flatIndex() const { return _flat_index; }
142 2673 : reference_type getMultiIndexObject() const { return _multi_index; }
143 :
144 : // assignment =
145 4 : const_noconst_iterator & operator=(const const_noconst_iterator & other)
146 : {
147 4 : _multi_index = other.getMultiIndexObject();
148 4 : _flat_index = other.flatIndex();
149 4 : _multi_index.findIndexVector(_flat_index, _indices);
150 4 : return *this;
151 : }
152 :
153 : // Compiler generated copy constructor
154 : const_noconst_iterator(const const_noconst_iterator &) = default;
155 :
156 : // prefix ++
157 22288 : const_noconst_iterator & operator++()
158 : {
159 22288 : ++_flat_index;
160 : // increment indices
161 41209 : for (unsigned int j = 0; j < _indices.size(); ++j)
162 : {
163 38530 : _indices[_indices.size() - j - 1] =
164 38530 : (_indices[_indices.size() - j - 1] + 1) % _shape[_indices.size() - j - 1];
165 38530 : if (_indices[_indices.size() - j - 1] != 0)
166 19609 : break;
167 : }
168 22288 : return *this;
169 : }
170 :
171 : // postfix ++
172 : const_noconst_iterator & operator++(int)
173 : {
174 : const_noconst_iterator clone(*this);
175 : ++_flat_index;
176 : // increment indices
177 : for (unsigned int j = 0; j < _indices.size(); ++j)
178 : {
179 : _indices[_indices.size() - j - 1] =
180 : (_indices[_indices.size() - j - 1] + 1) % _shape[_indices.size() - j - 1];
181 : if (_indices[_indices.size() - j - 1] != 0)
182 : break;
183 : }
184 : return clone;
185 : }
186 :
187 : // prefix --
188 24 : const_noconst_iterator & operator--()
189 : {
190 24 : --_flat_index;
191 : // decrement indices
192 34 : for (unsigned int j = 0; j < _indices.size(); ++j)
193 : {
194 33 : if (_indices[_indices.size() - j - 1] == 0)
195 10 : _indices[_indices.size() - j - 1] = _shape[_indices.size() - j - 1] - 1;
196 : else
197 : {
198 23 : --_indices[_indices.size() - j - 1];
199 23 : break;
200 : }
201 : }
202 24 : return *this;
203 : }
204 :
205 : // postfix --
206 : const_noconst_iterator & operator--(int)
207 : {
208 : const_noconst_iterator clone(*this);
209 : --_flat_index;
210 : // decrement indices
211 : for (unsigned int j = 0; j < _indices.size(); ++j)
212 : {
213 : if (_indices[_indices.size() - j - 1] == 0)
214 : _indices[_indices.size() - j - 1] = _shape[_indices.size() - j - 1] - 1;
215 : else
216 : {
217 : --_indices[_indices.size() - j - 1];
218 : break;
219 : }
220 : }
221 : return clone;
222 : }
223 :
224 : /// to be equal both iterators must hold a reference to the same MultiIndexObject and be at the same _flat_index
225 24453 : bool operator==(const const_noconst_iterator & other) const
226 : {
227 24453 : return _flat_index == other.flatIndex() && &_multi_index == &other.getMultiIndexObject();
228 : }
229 24453 : bool operator!=(const const_noconst_iterator & other) const { return !(*this == other); }
230 :
231 : /// dereferencing operator
232 22840 : std::pair<const size_type &, T &> operator*()
233 : {
234 22840 : return std::pair<const size_type &, T &>(_indices, _multi_index._data[_flat_index]);
235 : }
236 : std::pair<const size_type &, const T &> operator*() const
237 : {
238 : return std::pair<const size_type &, const T &>(_indices, _multi_index._data[_flat_index]);
239 : }
240 :
241 : protected:
242 : reference_type _multi_index;
243 : unsigned int _flat_index;
244 : size_type _shape;
245 : size_type _indices;
246 : };
247 :
248 : template <class T>
249 2684 : MultiIndex<T>::MultiIndex(const size_type & shape)
250 : {
251 2684 : reshape(shape);
252 2684 : _data.resize(_nentries);
253 2684 : }
254 :
255 : template <class T>
256 2 : MultiIndex<T>::MultiIndex(const size_type & shape, const std::vector<T> & data)
257 : {
258 2 : reshape(shape);
259 :
260 2 : if (data.size() != _nentries)
261 0 : mooseError("shape and data arguments' sizes are inconsistent.");
262 2 : _data = data;
263 2 : }
264 :
265 : template <class T>
266 : T &
267 690 : MultiIndex<T>::operator()(const size_type & indices)
268 : {
269 690 : return _data[flatIndex(indices)];
270 : }
271 :
272 : template <class T>
273 : const T &
274 : MultiIndex<T>::operator()(const size_type & indices) const
275 : {
276 : return _data[flatIndex(indices)];
277 : }
278 :
279 : template <class T>
280 : MultiIndex<T>
281 1 : MultiIndex<T>::slice(unsigned int dimension, unsigned int index) const
282 : {
283 1 : size_type dim(1, dimension);
284 1 : size_type ind(1, index);
285 2 : return slice(dim, ind);
286 1 : }
287 :
288 : template <class T>
289 : MultiIndex<T>
290 8 : MultiIndex<T>::slice(size_type dimension, size_type index) const
291 : {
292 : // input checks
293 8 : if (dimension.size() != index.size())
294 0 : mooseError("dimension and index must have the same size.");
295 8 : if (dimension.size() > _dim - 1)
296 0 : mooseError("The result of slice must be at least of dimension 1.");
297 :
298 : #if DEBUG
299 : for (unsigned int d = 0; d < dimension.size(); ++d)
300 : {
301 : if (dimension[d] >= _dim)
302 : mooseError("dimension is set to ", dimension[d], " which is larger than _dim ", _dim);
303 : if (index[d] >= _shape[dimension[d]])
304 : mooseError("index= ",
305 : index[d],
306 : " at dimension=",
307 : dimension[d],
308 : " is larger than ",
309 : _shape[dimension[d]]);
310 : }
311 : #endif // DEBUG
312 :
313 : // create a MultiIndex object with new dim = dim - dimension.size()
314 8 : size_type new_shape;
315 37 : for (unsigned int j = 0; j < _dim; ++j)
316 29 : if (std::find(dimension.begin(), dimension.end(), j) == dimension.end())
317 23 : new_shape.push_back(_shape[j]);
318 8 : MultiIndex<T> multi_index = MultiIndex(new_shape);
319 :
320 : // copy the data
321 8 : MultiIndex<T>::iterator it = multi_index.begin();
322 1388 : for (unsigned int n = 0; n < _nentries; ++n)
323 : {
324 1380 : size_type indices;
325 1380 : findIndexVector(n, indices);
326 1380 : bool addTo = true;
327 1728 : for (unsigned int d = 0; d < dimension.size(); ++d)
328 1272 : if (indices[dimension[d]] != index[d])
329 : {
330 924 : addTo = false;
331 924 : break;
332 : }
333 :
334 1380 : if (addTo)
335 : {
336 456 : (*it).second = _data[n];
337 456 : ++it;
338 : }
339 : }
340 16 : return multi_index;
341 8 : }
342 :
343 : template <class T>
344 : void
345 1 : MultiIndex<T>::resize(const size_type & shape)
346 : {
347 1 : if (shape.size() != _shape.size())
348 0 : mooseError("resize cannot change the dimensionality of MultiIndex.");
349 :
350 : // first copy the old shape and data
351 1 : size_type old_shape = _shape;
352 1 : size_type old_stride = _stride;
353 1 : std::vector<T> old_data = _data;
354 :
355 : // reset _shape & recompute meta data
356 1 : reshape(shape);
357 :
358 : // fill in _data to all possible indices
359 1 : _data.assign(_nentries, T(0));
360 31 : for (unsigned int j = 0; j < _nentries; ++j)
361 : {
362 30 : size_type indices;
363 30 : findIndexVector(j, indices);
364 :
365 : // check if indices existed in the old version of _data
366 30 : bool existed_in_old = true;
367 96 : for (unsigned int d = 0; d < _dim; ++d)
368 80 : if (indices[d] >= old_shape[d])
369 : {
370 14 : existed_in_old = false;
371 14 : break;
372 : }
373 :
374 : // find the corresponding old_j
375 30 : if (existed_in_old)
376 : {
377 16 : unsigned int old_j = 0;
378 64 : for (unsigned int d = 0; d < _dim; ++d)
379 48 : old_j += indices[d] * old_stride[d];
380 :
381 : // finally set the data entry
382 16 : _data[j] = old_data[old_j];
383 : }
384 : }
385 1 : }
386 :
387 : template <class T>
388 : unsigned int
389 3 : MultiIndex<T>::stride(unsigned int j) const
390 : {
391 : mooseAssert(j < _dim, "Dimension is" << _dim << ", stride(j) called with j = " << j);
392 3 : return _stride[j];
393 : }
394 :
395 : template <class T>
396 : void
397 1 : MultiIndex<T>::assign(const size_type & shape, T value)
398 : {
399 1 : reshape(shape);
400 1 : _data.assign(_nentries, value);
401 1 : }
402 :
403 : template <class T>
404 : void
405 1 : MultiIndex<T>::dataStore(std::ostream & stream, void * context)
406 : {
407 1 : ::dataStore(stream, _shape, context);
408 1 : ::dataStore(stream, _data, context);
409 1 : }
410 :
411 : template <class T>
412 : void
413 1 : MultiIndex<T>::dataLoad(std::istream & stream, void * context)
414 : {
415 1 : ::dataLoad(stream, _shape, context);
416 1 : ::dataLoad(stream, _data, context);
417 1 : _dim = _shape.size();
418 1 : _nentries = _data.size();
419 1 : buildAccumulatedShape();
420 1 : }
421 :
422 : template <class T>
423 : void
424 6760 : MultiIndex<T>::findIndexVector(unsigned int flat_index, MultiIndex<T>::size_type & indices) const
425 : {
426 6760 : indices.resize(_dim);
427 28241 : for (unsigned int d = 0; d < _dim; ++d)
428 : {
429 21481 : unsigned int i = flat_index / _stride[d];
430 21481 : indices[d] = i;
431 21481 : flat_index -= i * _stride[d];
432 : }
433 6760 : }
434 :
435 : // compute the accumulated shapes as:
436 : // as[0] = I_1 * I_2 ...* I_{M}, as[1] = I_2 * I_3 ...* I_{M} ...
437 : template <class T>
438 : void
439 2684 : MultiIndex<T>::buildAccumulatedShape()
440 : {
441 : // TODO: simplify - this is needlessly complicated - can be done in a single loop
442 2684 : _stride.resize(_dim);
443 10741 : for (unsigned int d = 0; d < _dim; ++d)
444 : {
445 8057 : unsigned int k = 1;
446 16127 : for (unsigned int j = d + 1; j < _dim; ++j)
447 8070 : k *= _shape[j];
448 8057 : _stride[d] = k;
449 : }
450 2684 : }
451 :
452 : template <class T>
453 : void
454 2688 : MultiIndex<T>::reshape(const size_type & shape)
455 : {
456 2688 : if (shape.size() == 0)
457 : {
458 5 : _shape = {};
459 5 : _dim = 0;
460 5 : _stride = {};
461 5 : _nentries = 0;
462 5 : return;
463 : }
464 :
465 2683 : _shape = shape;
466 2683 : _dim = shape.size();
467 :
468 2683 : _nentries = 1;
469 10737 : for (unsigned int d = 0; d < _dim; ++d)
470 8054 : _nentries *= _shape[d];
471 :
472 2683 : buildAccumulatedShape();
473 : }
474 :
475 : template <class T>
476 : unsigned int
477 3354 : MultiIndex<T>::flatIndex(const size_type & indices) const
478 : {
479 : mooseAssert(indices.size() == _dim,
480 : "Indices vector has wrong size. size=" << indices.size() << " vs. dim=" << _dim);
481 : #if DEBUG
482 : for (unsigned int j = 0; j < indices.size(); ++j)
483 : if (indices[j] >= _shape[j])
484 : mooseError("Indices vector at entry ", j, " is ", indices[j], " vs. shape ", _shape[j]);
485 : #endif // DEBUG
486 :
487 : // implement the index
488 : // index = i_M + i_{M-1} * I_M + i_{M-1} * I_M * I_{M-1} ...
489 3354 : unsigned int index = 0;
490 13588 : for (unsigned int d = 0; d < _dim; ++d)
491 10234 : index += indices[d] * _stride[d];
492 :
493 3354 : return index;
494 : }
495 :
496 : template <class T>
497 : void
498 1 : dataStore(std::ostream & stream, MultiIndex<T> & mi, void * context)
499 : {
500 1 : mi.dataStore(stream, context);
501 1 : }
502 :
503 : template <class T>
504 : void
505 1 : dataLoad(std::istream & stream, MultiIndex<T> & mi, void * context)
506 : {
507 1 : mi.dataLoad(stream, context);
508 1 : }
|
