Line data Source code
1 : // This file is part of the Mamico project. For conditions of distribution
2 : // and use, please see the copyright notice in Mamico's main folder
3 :
4 : #pragma once
5 :
6 : #include <memory>
7 :
8 : namespace coupling {
9 :
10 : /** Modular filtering and data analytics system. Includes inferfaces and
11 : * services that manage the filtering subsystem, as well as several filter
12 : * implementations.
13 : */
14 : namespace filtering {
15 :
16 : /***
17 : * Quantities<dim>[0] := mass
18 : * Quantities<dim>[1..dim] := momentum
19 : */
20 : template <unsigned int dim> using Quantities = tarch::la::Vector<dim + 1, double>;
21 :
22 : /***
23 : * Spacetime window of data, i.e. 4D field of T
24 : */
25 : template <unsigned int dim, typename T> class Field;
26 :
27 : /***
28 : * Spacetime window of flow quantities
29 : */
30 : template <unsigned int dim> using Flowfield = Field<dim, Quantities<dim>>;
31 :
32 : /***
33 : * includes local flowfield, and mean and standard deviation of quantities
34 : */
35 : template <unsigned int dim> class Patch;
36 :
37 : /***
38 : * similar to patch, but
39 : * does not store local flowfield, but accesses basefield for distance
40 : * computation
41 : */
42 : template <unsigned int dim> class PatchView;
43 :
44 : /***
45 : * Spacetime window of patches
46 : */
47 : template <unsigned int dim> using Patchfield = Field<dim, Patch<dim>>;
48 : // using Patchfield = Field<dim, PatchView<dim>>;
49 :
50 : } // namespace filtering
51 : } // namespace coupling
52 :
53 : /***
54 : * Spacetime window of data, i.e. 4D field of T
55 : */
56 : template <unsigned int dim, typename T> class coupling::filtering::Field {
57 : public:
58 0 : Field(const tarch::la::Vector<dim, unsigned int>& spatialSize, const unsigned int& temporalSize)
59 0 : : _spatialSize(spatialSize), _temporalSize(temporalSize), _scalarSize(computeScalarSize(spatialSize, temporalSize)),
60 0 : _data(std::allocator_traits<std::allocator<T>>::allocate(allo, _scalarSize)) {}
61 :
62 0 : template <class... Args> void construct(const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t, Args&&... args) {
63 0 : std::allocator_traits<std::allocator<T>>::construct(allo, _data + idx(pos, t), std::forward<Args>(args)...);
64 0 : }
65 :
66 0 : void destroy(const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t) {
67 0 : std::allocator_traits<std::allocator<T>>::destroy(allo, _data + idx(pos, t));
68 0 : }
69 :
70 0 : T& operator()(const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t) { return _data[idx(pos, t)]; }
71 :
72 0 : const T& operator()(const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t) const { return _data[idx(pos, t)]; }
73 :
74 0 : T& operator[](unsigned int pos) {
75 : #ifdef NLM_DEBUG
76 : if (pos < 0 || pos > _scalarSize - 1) {
77 : std::cout << "ERROR Field T& operator[](int pos): pos out of range!" << std::endl;
78 : std::cout << "pos=" << pos << ", _scalarSize=" << _scalarSize << std::endl;
79 : exit(EXIT_FAILURE);
80 : }
81 : #endif
82 0 : return _data[pos];
83 : }
84 :
85 : const T& operator[](unsigned int pos) const {
86 : #ifdef NLM_DEBUG
87 : if (pos < 0 || pos > _scalarSize - 1) {
88 : std::cout << "ERROR Field T& operator[](int pos): pos out of range!" << std::endl;
89 : std::cout << "pos=" << pos << ", _scalarSize=" << _scalarSize << std::endl;
90 : exit(EXIT_FAILURE);
91 : }
92 : #endif
93 : return _data[pos];
94 : }
95 :
96 0 : ~Field() { std::allocator_traits<std::allocator<T>>::deallocate(allo, _data, _scalarSize); }
97 :
98 0 : unsigned int getScalarSize() const { return _scalarSize; }
99 :
100 0 : const tarch::la::Vector<dim, unsigned int>& getSpatialSize() const { return _spatialSize; }
101 :
102 0 : unsigned int getTemporalSize() const { return _temporalSize; }
103 :
104 : void print() {
105 : for (unsigned int i = 0; i < _scalarSize; i++) {
106 : std::cout << "Field elem " << i << " = " << _data[i] << std::endl;
107 : }
108 : }
109 :
110 : private:
111 0 : unsigned int computeScalarSize(const tarch::la::Vector<dim, unsigned int>& spatialSize, const unsigned int& temporalSize) const {
112 0 : unsigned int res = spatialSize[0];
113 0 : for (unsigned int d = 1; d < dim; d++) {
114 0 : res *= (spatialSize[d]);
115 : }
116 0 : res *= temporalSize;
117 : return res;
118 : }
119 :
120 0 : unsigned int idx(const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t) const {
121 : #ifdef NLM_DEBUG
122 : for (unsigned int d = 0; d < dim; d++) {
123 : if (pos[d] < 0 || pos[d] > _spatialSize[d] - 1) {
124 : std::cout << "ERROR Field idx(): pos out of range!" << std::endl;
125 : std::cout << "pos=" << pos << ", _spatialSize=" << _spatialSize << std::endl;
126 : exit(EXIT_FAILURE);
127 : }
128 : }
129 : if (t < 0 || t > _temporalSize - 1) {
130 : std::cout << "ERROR Field idx(): t out of range!" << std::endl;
131 : std::cout << "t=" << t << ", _temporalSize=" << _temporalSize << std::endl;
132 : exit(EXIT_FAILURE);
133 : }
134 : #endif
135 0 : unsigned int idx = 0, step = 1;
136 0 : for (unsigned int d = 0; d < dim; d++) {
137 0 : idx += pos[d] * step;
138 0 : step *= _spatialSize[d];
139 : }
140 0 : idx += t * step;
141 0 : return idx;
142 : }
143 :
144 : static std::allocator<T> allo;
145 : const tarch::la::Vector<dim, unsigned int> _spatialSize;
146 : const unsigned int _temporalSize;
147 : const unsigned int _scalarSize;
148 : T* const _data;
149 :
150 : friend class coupling::filtering::Patch<dim>; // patches need direct access to
151 : // their field's _data for
152 : // memory efficiency
153 : };
154 :
155 : template <unsigned int dim, typename T> std::allocator<T> coupling::filtering::Field<dim, T>::allo;
156 :
157 : /***
158 : * includes local flowfield, and mean and standard deviation of quantities
159 : */
160 : template <unsigned int dim> class coupling::filtering::Patch {
161 : public:
162 0 : Patch(const tarch::la::Vector<dim, unsigned int>& spatialSize, const unsigned int& temporalSize, const Flowfield<dim>& basefield,
163 : const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t)
164 0 : : _flowfield(spatialSize, temporalSize), _localMean(0.0), _localStandardDeviation(0.0) {
165 0 : fillFromBasefield(basefield, pos, t);
166 0 : computeLocalMean();
167 0 : computeLocalStandardDeviation();
168 0 : }
169 :
170 0 : double distance(const coupling::filtering::Patch<dim>& other) {
171 0 : unsigned int size = _flowfield.getScalarSize() * (dim + 1);
172 :
173 0 : double* const my_data = reinterpret_cast<double* const>(_flowfield._data);
174 0 : double* const other_data = reinterpret_cast<double* const>(other._flowfield._data);
175 :
176 0 : double res = 0;
177 0 : for (unsigned int i = 0; i < size; i += 1) {
178 0 : double diff(my_data[i] - other_data[i]);
179 0 : res += diff * diff;
180 : }
181 0 : return res;
182 : }
183 :
184 0 : ~Patch() {}
185 :
186 : void print() { _flowfield.print(); }
187 :
188 : private:
189 0 : inline unsigned int posmod(int i, int n) { return (i % n + n) % n; }
190 :
191 0 : void fillFromBasefield(const Flowfield<dim>& basefield, const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t) {
192 : static_assert(dim == 2 || dim == 3, "ERROR filtering::Patch::fillFromBasefield only implemented "
193 : "for 2D and 3D");
194 :
195 0 : tarch::la::Vector<dim, unsigned int> center;
196 0 : for (unsigned int d = 0; d < dim; d++) {
197 0 : center[d] = _flowfield.getSpatialSize()[d] / 2;
198 : }
199 :
200 0 : tarch::la::Vector<dim, unsigned int> local_pos(0);
201 0 : unsigned int local_t(0);
202 0 : tarch::la::Vector<dim, unsigned int> base_pos(0);
203 0 : unsigned int base_t(0);
204 :
205 0 : for (int offset_t = 0; offset_t > -(int)_flowfield.getTemporalSize(); offset_t--) {
206 0 : local_t = posmod(offset_t, _flowfield.getTemporalSize());
207 0 : base_t = posmod(t + offset_t, basefield.getTemporalSize());
208 0 : for (int offset_x = -(int)(_flowfield.getSpatialSize()[0]) / 2; offset_x <= (int)(_flowfield.getSpatialSize()[0]) / 2; offset_x++) {
209 0 : local_pos[0] = center[0] + offset_x;
210 0 : base_pos[0] = pos[0] + offset_x;
211 0 : for (int offset_y = -(int)(_flowfield.getSpatialSize()[1]) / 2; offset_y <= (int)(_flowfield.getSpatialSize()[1]) / 2; offset_y++) {
212 0 : local_pos[1] = center[1] + offset_y;
213 0 : base_pos[1] = pos[1] + offset_y;
214 : if constexpr (dim == 3) {
215 0 : for (int offset_z = -(int)(_flowfield.getSpatialSize()[2]) / 2; offset_z <= (int)(_flowfield.getSpatialSize()[2]) / 2; offset_z++) {
216 0 : local_pos[2] = center[2] + offset_z;
217 0 : base_pos[2] = pos[2] + offset_z;
218 0 : _flowfield(local_pos, local_t) = basefield(base_pos, base_t);
219 : }
220 : } else {
221 : _flowfield(local_pos, local_t) = basefield(base_pos, base_t);
222 : }
223 : }
224 : }
225 : }
226 0 : }
227 :
228 0 : void computeLocalMean() {
229 0 : for (unsigned int i = 0; i < _flowfield.getScalarSize(); i++)
230 0 : _localMean += _flowfield[i];
231 0 : _localMean = _localMean * (1.0 / _flowfield.getScalarSize());
232 0 : }
233 :
234 0 : void computeLocalStandardDeviation() {
235 0 : for (unsigned int i = 0; i < _flowfield.getScalarSize(); i++) {
236 0 : Quantities<dim> diff = _localMean - _flowfield[i];
237 0 : _localStandardDeviation += tarch::la::dot(diff, diff);
238 : }
239 0 : _localStandardDeviation = sqrt(_localStandardDeviation / _flowfield.getScalarSize());
240 0 : }
241 :
242 : Flowfield<dim> _flowfield;
243 : Quantities<dim> _localMean;
244 : double _localStandardDeviation;
245 : };
246 :
247 : /***
248 : * similar to patch, but does not store local flowfield,
249 : * accesses basefield for distance computation
250 : */
251 : template <unsigned int dim> class coupling::filtering::PatchView {
252 : public:
253 : PatchView(const tarch::la::Vector<dim, unsigned int>& spatialSize, const unsigned int& temporalSize, const Flowfield<dim>& basefield,
254 : const tarch::la::Vector<dim, unsigned int>& pos, const unsigned int& t)
255 : : _spatialSize(spatialSize), _temporalSize(temporalSize), _basefield(basefield), _pos(pos), _t(t) {}
256 :
257 : double distance(const coupling::filtering::PatchView<dim>& other) const {
258 : double res = 0;
259 :
260 : tarch::la::Vector<dim, unsigned int> base_pos(0);
261 : unsigned int base_t(0);
262 :
263 : for (int offset_t = 0; offset_t > -(int)_temporalSize; offset_t--) {
264 : base_t = posmod(_t + offset_t, _basefield.getTemporalSize());
265 : for (int offset_x = -(int)(_spatialSize[0]) / 2; offset_x <= (int)(_spatialSize[0]) / 2; offset_x++) {
266 : base_pos[0] = _pos[0] + offset_x;
267 : for (int offset_y = -(int)(_spatialSize[1]) / 2; offset_y <= (int)(_spatialSize[1]) / 2; offset_y++) {
268 : base_pos[1] = _pos[1] + offset_y;
269 : if constexpr (dim == 3) {
270 : for (int offset_z = -(int)(_spatialSize[2]) / 2; offset_z <= (int)(_spatialSize[2]) / 2; offset_z++) {
271 : base_pos[2] = _pos[2] + offset_z;
272 : // TODO: this is not correct. fix me
273 : auto diff(_basefield(base_pos, base_t) - _basefield(base_pos + other._pos - _pos, posmod(base_t + other._t - _t, _basefield.getTemporalSize())));
274 : res += diff[0] * diff[0];
275 : }
276 : } else {
277 : auto diff(_basefield(base_pos, base_t) - other._basefield(base_pos, base_t));
278 : res += tarch::la::dot(diff, diff);
279 : }
280 : }
281 : }
282 : }
283 : return res;
284 : }
285 :
286 : private:
287 : inline unsigned int posmod(int i, int n) const { return (i % n + n) % n; }
288 :
289 : const tarch::la::Vector<dim, unsigned int> _spatialSize;
290 : const unsigned int _temporalSize;
291 : const Flowfield<dim>& _basefield;
292 : const tarch::la::Vector<dim, unsigned int> _pos;
293 : const unsigned int _t;
294 : };
|
