source: src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

Candidate_v1.6.1
Last change on this file was 9eb71b3, checked in by Frederik Heber <frederik.heber@…>, 8 years ago

Commented out MemDebug include and Memory::ignore.

  • MemDebug clashes with various allocation operators that use a specific placement in memory. It is so far not possible to wrap new/delete fully. Hence, we stop this effort which so far has forced us to put ever more includes (with clashes) into MemDebug and thereby bloat compilation time.
  • MemDebug does not add that much usefulness which is not also provided by valgrind.
  • Property mode set to 100644
File size: 36.2 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder 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
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * SamplingGrid.cpp
26 *
27 * Created on: 25.07.2012
28 * Author: heber
29 */
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36// include headers that implement a archive in simple text format
37// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
38#include <boost/archive/text_oarchive.hpp>
39#include <boost/archive/text_iarchive.hpp>
40
41//#include "CodePatterns/MemDebug.hpp"
42
43#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
44
45#include <boost/assign.hpp>
46#include <boost/bind.hpp>
47#include <boost/shared_ptr.hpp>
48
49#include <algorithm>
50#include <limits>
51
52#include "CodePatterns/Assert.hpp"
53#include "CodePatterns/Log.hpp"
54
55#include "LinearAlgebra/Vector.hpp"
56
57// static instances
58const double SamplingGrid::zeroOffset[NDIM] = { 0., 0., 0. };
59
60SamplingGrid::SamplingGrid() :
61 SamplingGridProperties()
62{
63 setWindowSize(zeroOffset, zeroOffset);
64 ASSERT( getWindowGridPoints() == (size_t)0,
65 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
66}
67
68SamplingGrid::SamplingGrid(const double _begin[NDIM],
69 const double _end[NDIM],
70 const int _level) :
71 SamplingGridProperties(_begin, _end, _level)
72{
73 setWindowSize(zeroOffset, zeroOffset);
74 ASSERT( getWindowGridPoints() == (size_t)0,
75 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
76}
77
78SamplingGrid::SamplingGrid(const double _begin[NDIM],
79 const double _end[NDIM],
80 const int _level,
81 const sampledvalues_t &_sampled_grid) :
82 SamplingGridProperties(_begin, _end, _level),
83 sampled_grid(_sampled_grid)
84{
85 setWindowSize(_begin, _end);
86 ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
87 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
88}
89
90SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
91 SamplingGridProperties(_grid),
92 sampled_grid(_grid.sampled_grid)
93{
94 setWindowSize(_grid.begin_window, _grid.end_window);
95 ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(),
96 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
97}
98
99SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
100 SamplingGridProperties(_props)
101{
102 setWindowSize(zeroOffset, zeroOffset);
103 ASSERT( getWindowGridPoints() == (size_t)0,
104 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
105}
106
107SamplingGrid::SamplingGrid(
108 const SamplingGridProperties &_props,
109 const sampledvalues_t &_sampled_grid) :
110 SamplingGridProperties(_props),
111 sampled_grid(_sampled_grid)
112{
113 setWindowSize(_props.begin, _props.end);
114 ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
115 "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
116}
117
118SamplingGrid::~SamplingGrid()
119{}
120
121bool SamplingGrid::isCongruent(const SamplingGrid &_props) const
122{
123 bool status = true;
124 status &= (static_cast<const SamplingGridProperties &>(*this) ==
125 static_cast<const SamplingGridProperties &>(_props));
126 for(size_t i = 0; i<NDIM; ++i) {
127 status &= begin_window[i] == _props.begin_window[i];
128 status &= end_window[i] == _props.end_window[i];
129 }
130 return status;
131}
132
133SamplingGrid& SamplingGrid::operator=(const SamplingGrid& other)
134{
135 // check for self-assignment
136 if (this != &other) {
137 static_cast<SamplingGridProperties &>(*this) = other;
138 setWindowSize(other.begin_window, other.end_window);
139 sampled_grid = other.sampled_grid;
140 }
141 return *this;
142}
143
144static void multiplyElements(
145 double &dest,
146 const double &source,
147 const double prefactor)
148{
149 dest *= prefactor*(source);
150}
151
152SamplingGrid& SamplingGrid::operator*=(const double _value)
153{
154 std::transform(
155 sampled_grid.begin(), sampled_grid.end(),
156 sampled_grid.begin(),
157 boost::bind(std::multiplies<double>(), _1, _value));
158
159 return *this;
160}
161
162static void getDownsampledGrid(
163 const SamplingGrid &_reference_grid,
164 const SamplingGrid &_grid,
165 boost::shared_ptr<SamplingGrid> &_weight_downsampled)
166{
167 static const double round_offset(
168 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
169 0.5 : 0.); // need offset to get to round_toward_nearest behavior
170 const int surplus_level = _reference_grid.getSurplusLevel(_grid)+round_offset;
171 // need to downsample first
172 _weight_downsampled.reset(new SamplingGrid); // may use default cstor as downsample does all settings
173 SamplingGrid::downsample(*_weight_downsampled, _grid, _reference_grid.level-surplus_level);
174}
175
176SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
177{
178 // check that grids are compatible
179 ASSERT(isCompatible(other),
180 "SamplingGrid::operator*=() - multiplying incomatible grids is so far not in the cards.");
181 const SamplingGrid *other_grid = &other;
182 boost::shared_ptr<SamplingGrid> other_downsampled;
183 if (!isEquivalent(other)) {
184 getDownsampledGrid(*this, other, other_downsampled);
185 other_grid = other_downsampled.get();
186 }
187 /// get minimum of window
188 double min_begin_window[NDIM];
189 double min_end_window[NDIM];
190 bool doShrink = false;
191 for (size_t index=0; index<NDIM;++index) {
192 if (begin_window[index] <= other.begin_window[index]) {
193 min_begin_window[index] = other.begin_window[index];
194 doShrink = true;
195 } else {
196 min_begin_window[index] = begin_window[index];
197 }
198 if (end_window[index] <= other.end_window[index]) {
199 min_end_window[index] = end_window[index];
200 } else {
201 min_end_window[index] = other.end_window[index];
202 doShrink = true;
203 }
204 }
205 LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
206 LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
207 if (doShrink)
208 shrinkWindow(min_begin_window, min_end_window);
209 addWindowOntoWindow(
210 other_grid->begin_window,
211 other_grid->end_window,
212 begin_window,
213 end_window,
214 sampled_grid,
215 other_grid->sampled_grid,
216 boost::bind(multiplyElements, _1, _2, 1.),
217 sourcewindow);
218
219 return *this;
220}
221
222void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
223{
224 // check that grids are compatible
225 ASSERT(isCompatible(other),
226 "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
227 const SamplingGrid *other_grid = &other;
228 boost::shared_ptr<SamplingGrid> other_downsampled;
229 if (!isEquivalent(other)) {
230 getDownsampledGrid(*this, other, other_downsampled);
231 other_grid = other_downsampled.get();
232 }
233 /// get maximum of window
234 double max_begin_window[NDIM];
235 double max_end_window[NDIM];
236 bool doExtend = false;
237 for (size_t index=0; index<NDIM;++index) {
238 if (begin_window[index] >= other.begin_window[index]) {
239 max_begin_window[index] = other.begin_window[index];
240 doExtend = true;
241 } else {
242 max_begin_window[index] = begin_window[index];
243 }
244 if (end_window[index] >= other.end_window[index]) {
245 max_end_window[index] = end_window[index];
246 } else {
247 max_end_window[index] = other.end_window[index];
248 doExtend = true;
249 }
250 }
251 LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
252 LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
253 if (doExtend)
254 extendWindow(max_begin_window, max_end_window);
255 /// and copy other into larger window, too
256 addOntoWindow(other_grid->begin_window, other_grid->end_window, other_grid->sampled_grid, prefactor);
257}
258
259const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const
260{
261 static const double round_offset(
262 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
263 0.5 : 0.); // need offset to get to round_toward_nearest behavior
264// const double total = getTotalLengthPerAxis(axis);
265 const double delta = getDeltaPerAxis(axis);
266 if (delta == 0)
267 return 0;
268 const double length = getWindowLengthPerAxis(axis);
269 if (length == 0)
270 return 0;
271 return (size_t)(length/delta+round_offset);
272}
273
274double SamplingGrid::integral() const
275{
276 const double volume_element = getVolume()/(double)getTotalGridPoints();
277 double int_value = 0.;
278 for (sampledvalues_t::const_iterator iter = sampled_grid.begin();
279 iter != sampled_grid.end();
280 ++iter)
281 int_value += *iter;
282 int_value *= volume_element;
283 LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
284 return int_value;
285}
286
287double SamplingGrid::integral(const SamplingGrid &weight) const
288{
289 // check that grids are compatible
290 ASSERT(isCompatible(weight),
291 "SamplingGrid::integral() - integrating with weights from incompatible grids is so far not in the cards.");
292 const SamplingGrid *weight_grid = &weight;
293 boost::shared_ptr<SamplingGrid> weight_downsampled;
294 if (!isEquivalent(weight)) {
295 getDownsampledGrid(*this, weight, weight_downsampled);
296 weight_grid = weight_downsampled.get();
297 }
298 const double volume_element = getVolume()/(double)getTotalGridPoints();
299 double int_value = 0.;
300 sampledvalues_t::const_iterator iter = sampled_grid.begin();
301 sampledvalues_t::const_iterator weightiter = weight_grid->sampled_grid.begin();
302 for (;iter != sampled_grid.end();++iter,++weightiter)
303 int_value += (*weightiter) * (*iter);
304 int_value *= volume_element;
305 //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
306 return int_value;
307}
308
309void SamplingGrid::setWindowSize(
310 const double _begin_window[NDIM],
311 const double _end_window[NDIM])
312{
313 for (size_t index=0;index<NDIM;++index) {
314 begin_window[index] = getNearestLowerGridPoint(_begin_window[index], index);
315 ASSERT( begin_window[index] >= begin[index],
316 "SamplingGrid::setWindowSize() - window starts earlier than domain on "
317 +toString(index)+"th component.");
318 end_window[index] = getNearestHigherGridPoint(_end_window[index], index);
319 ASSERT( end_window[index] <= end[index],
320 "SamplingGrid::setWindowSize() - window ends later than domain on "
321 +toString(index)+"th component.");
322 }
323}
324
325void SamplingGrid::setWindow(
326 const double _begin_window[NDIM],
327 const double _end_window[NDIM])
328{
329 setWindowSize(_begin_window, _end_window);
330 const size_t gridpoints_window = getWindowGridPoints();
331 sampled_grid.clear();
332 sampled_grid.resize(gridpoints_window, 0.);
333}
334
335void SamplingGrid::setDomain(
336 const double _begin[NDIM],
337 const double _end[NDIM])
338{
339 setDomainSize(_begin, _end);
340 setWindowSize(_begin, _end);
341 const size_t gridpoints = getTotalGridPoints();
342 sampled_grid.resize(gridpoints, 0.);
343}
344
345void SamplingGrid::extendWindow(
346 const double _begin_window[NDIM],
347 const double _end_window[NDIM])
348{
349#ifndef NDEBUG
350 for(size_t index=0;index < NDIM; ++index) {
351 // check that we truly have to extend the window
352 ASSERT ( begin_window[index] >= _begin_window[index],
353 "SamplingGrid::extendWindow() - component "+toString(index)+
354 " of window start is greater than old value.");
355 ASSERT ( end_window[index] <= _end_window[index],
356 "SamplingGrid::extendWindow() - component "+toString(index)+
357 " of window end is less than old value.");
358
359 // check that we are still less than domain
360 ASSERT ( _begin_window[index] >= begin[index],
361 "SamplingGrid::extendWindow() - component "+toString(index)+
362 " of window start is less than domain start.");
363 ASSERT ( _end_window[index] <= end[index],
364 "SamplingGrid::extendWindow() - component "+toString(index)+
365 " of window end is greater than domain end.");
366 }
367#endif
368 // copy old window size and values
369 double old_begin_window[NDIM];
370 double old_end_window[NDIM];
371 for(size_t index=0;index<NDIM;++index) {
372 old_begin_window[index] = begin_window[index];
373 old_end_window[index] = end_window[index];
374 }
375 sampledvalues_t old_values(sampled_grid);
376 // set new window
377 setWindow(_begin_window,_end_window);
378 // now extend it ...
379 addOntoWindow(old_begin_window, old_end_window, old_values, +1.);
380 LOG(6, "DEBUG: Grid after extension is " << sampled_grid << ".");
381}
382
383void SamplingGrid::shrinkWindow(
384 const double _begin_window[NDIM],
385 const double _end_window[NDIM])
386{
387#ifndef NDEBUG
388 for(size_t index=0;index < NDIM; ++index) {
389 // check that we truly have to shrink the window
390 ASSERT ( begin_window[index] <= _begin_window[index],
391 "SamplingGrid::shrinkWindow() - component "+toString(index)+
392 " of window start is less than old value.");
393 ASSERT ( end_window[index] >= _end_window[index],
394 "SamplingGrid::shrinkWindow() - component "+toString(index)+
395 " of window end is greater than old value.");
396
397 // check that we are still less than domain
398 ASSERT ( _begin_window[index] >= begin[index],
399 "SamplingGrid::shrinkWindow() - component "+toString(index)+
400 " of window start is less than domain start.");
401 ASSERT ( _end_window[index] <= end[index],
402 "SamplingGrid::shrinkWindow() - component "+toString(index)+
403 " of window end is greater than domain end.");
404 }
405#endif
406 // copy old window size and values
407 double old_begin_window[NDIM];
408 double old_end_window[NDIM];
409 for(size_t index=0;index<NDIM;++index) {
410 old_begin_window[index] = begin_window[index];
411 old_end_window[index] = end_window[index];
412 }
413 sampledvalues_t old_values(sampled_grid);
414 // set new window
415 setWindow(_begin_window,_end_window);
416 // now extend it ...
417 addIntoWindow(old_begin_window, old_end_window, old_values, +1.);
418 LOG(6, "DEBUG: Grid after extension is " << sampled_grid << ".");
419}
420
421static void addElements(
422 double &dest,
423 const double &source,
424 const double prefactor)
425{
426 dest += prefactor*(source);
427}
428
429void SamplingGrid::addOntoWindow(
430 const double _begin_window[NDIM],
431 const double _end_window[NDIM],
432 const sampledvalues_t &_sampled_grid,
433 const double prefactor)
434{
435 addWindowOntoWindow(
436 begin_window,
437 end_window,
438 _begin_window,
439 _end_window,
440 sampled_grid,
441 _sampled_grid,
442 boost::bind(addElements, _1, _2, boost::cref(prefactor)),
443 destwindow);
444}
445
446void SamplingGrid::addIntoWindow(
447 const double _begin_window[NDIM],
448 const double _end_window[NDIM],
449 const sampledvalues_t &_sampled_grid,
450 const double prefactor)
451{
452 addWindowOntoWindow(
453 _begin_window,
454 _end_window,
455 begin_window,
456 end_window,
457 sampled_grid,
458 _sampled_grid,
459 boost::bind(addElements, _1, _2, boost::cref(prefactor)),
460 sourcewindow);
461}
462
463void SamplingGrid::getDiscreteWindowIndices(
464 size_t _wbegin[NDIM],
465 size_t _wlength[NDIM],
466 size_t _wend[NDIM]) const
467{
468 const double round_offset =
469 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
470 0.5 : 0.; // need offset to get to round_toward_nearest behavior
471 for(size_t index=0;index<NDIM;++index) {
472 if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
473 // we refrain from using floor/ceil as the window's starts and ends,
474 // the grids have to be compatible (equal level), should always be on
475 // discrete grid point locations.
476 const double delta = getDeltaPerAxis(index);
477 // delta is conversion factor from box length to discrete length, i.e. number of points
478 _wbegin[index] = (begin_window[index] - begin[index])/delta+round_offset;
479 _wlength[index] = (end_window[index] - begin_window[index])/delta+round_offset;
480 _wend[index] = (end_window[index] - begin[index])/delta+round_offset;
481 } else {
482 _wbegin[index] = 0;
483 _wlength[index] = 0;
484 _wend[index] = 0;
485 }
486 // total is used as safe-guard against loss due to discrete conversion
487 ASSERT( (_wend[index] - _wbegin[index]) == _wlength[index],
488 "SamplingGrid::getDiscreteWindowCopyIndices() - end - begin is not equal to length for "
489 +toString(index)+"th component.");
490 }
491}
492
493void SamplingGrid::getDiscreteWindowOffsets(
494 size_t _pre_offset[NDIM],
495 size_t _post_offset[NDIM],
496 size_t _length[NDIM],
497 size_t _total[NDIM]) const
498{
499 const double round_offset =
500 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
501 0.5 : 0.; // need offset to get to round_toward_nearest behavior
502 for(size_t index=0;index<NDIM;++index) {
503 if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
504 // we refrain from using floor/ceil as the window's starts and ends,
505 // the grids have to be compatible (equal level), should always be on
506 // discrete grid point locations.
507 const double delta = getDeltaPerAxis(index);
508 // delta is conversion factor from box length to discrete length, i.e. number of points
509 _pre_offset[index] = (begin_window[index] - begin[index])/delta+round_offset;
510 _post_offset[index] = (end[index] - end_window[index])/delta+round_offset;
511 _length[index] = (end_window[index] - begin_window[index])/delta+round_offset;
512 _total[index] = (end[index] - begin[index])/delta+round_offset;
513 } else {
514 _pre_offset[index] = 0;
515 _post_offset[index] = 0;
516 _length[index] = 0;
517 _total[index] = 0;
518 }
519 // total is used as safe-guard against loss due to discrete conversion
520 ASSERT( (_pre_offset[index] + _post_offset[index]) + _length[index] == _total[index],
521 "SamplingGrid::getDiscreteWindowCopyIndices() - pre, length, post are not equal to total for "
522 +toString(index)+"th component.");
523 }
524}
525
526void SamplingGrid::getDiscreteWindowCopyIndices(
527 const double *larger_wbegin,
528 const double *larger_wend,
529 const double *smaller_wbegin,
530 const double *smaller_wend,
531 size_t *pre_offset,
532 size_t *post_offset,
533 size_t *length,
534 size_t *total) const
535{
536 const double round_offset =
537 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
538 0.5 : 0.; // need offset to get to round_toward_nearest behavior
539 for(size_t index=0;index<NDIM;++index) {
540 if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
541 // we refrain from using floor/ceil as the window's starts and ends,
542 // the grids have to be compatible (equal level), should always be on
543 // discrete grid point locations.
544 const double delta = getDeltaPerAxis(index);
545 // delta is conversion factor from box length to discrete length, i.e. number of points
546 pre_offset[index] = (smaller_wbegin[index] - larger_wbegin[index])/delta+round_offset;
547 length[index] = (smaller_wend[index] - smaller_wbegin[index])/delta+round_offset;
548 post_offset[index] = (larger_wend[index] - smaller_wend[index])/delta+round_offset;
549 total[index] = (larger_wend[index] - larger_wbegin[index])/delta+round_offset;
550 } else {
551 pre_offset[index] = 0;
552 length[index] = 0;
553 post_offset[index] = 0;
554 total[index] = 0;
555 }
556 // total is used as safe-guard against loss due to discrete conversion
557 ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
558 "SamplingGrid::getDiscreteWindowCopyIndices() - pre, post, and length don't sum up to total for "
559 +toString(index)+"th component.");
560 }
561}
562
563void SamplingGrid::addWindowOntoWindow(
564 const double larger_wbegin[NDIM],
565 const double larger_wend[NDIM],
566 const double smaller_wbegin[NDIM],
567 const double smaller_wend[NDIM],
568 sampledvalues_t &dest_sampled_grid,
569 const sampledvalues_t &source_sampled_grid,
570 boost::function<void (double &, const double &)> op,
571 enum eLargerWindow larger_window)
572{
573#ifndef NDEBUG
574 for(size_t index=0;index<NDIM;++index) {
575 ASSERT( smaller_wbegin[index] >= larger_wbegin[index],
576 "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
577 +toString(index)+".");
578 ASSERT( smaller_wend[index] <= larger_wend[index],
579 "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
580 +toString(index)+".");
581 }
582#endif
583 // the only issue are indices
584 size_t pre_offset[NDIM];
585 size_t post_offset[NDIM];
586 size_t length[NDIM];
587 size_t total[NDIM];
588 getDiscreteWindowCopyIndices(
589 larger_wbegin, larger_wend,
590 smaller_wbegin, smaller_wend,
591 pre_offset,
592 post_offset,
593 length,
594 total
595 );
596 // assert that calculated lengths match with given vector sizes
597#ifndef NDEBUG
598 const size_t calculated_size = length[0]*length[1]*length[2];
599 if (larger_window == destwindow) {
600 ASSERT( calculated_size == source_sampled_grid.size(),
601 "SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: "
602 +toString(calculated_size)+" != "+toString(source_sampled_grid.size())+".");
603 ASSERT( calculated_size <= dest_sampled_grid.size(),
604 "SamplingGrid::addWindowOntoWindow() - not enough sampled values available: "
605 +toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+".");
606 } else {
607 ASSERT( calculated_size == dest_sampled_grid.size(),
608 "SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: "
609 +toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+".");
610 ASSERT( calculated_size <= source_sampled_grid.size(),
611 "SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: "
612 +toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+".");
613 }
614 const size_t total_size = total[0]*total[1]*total[2];
615 if (larger_window == destwindow) {
616 ASSERT( total_size == dest_sampled_grid.size(),
617 "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: "
618 +toString(total_size)+" != "+toString(dest_sampled_grid.size())+".");
619 } else {
620 ASSERT( total_size == source_sampled_grid.size(),
621 "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: "
622 +toString(total_size)+" != "+toString(source_sampled_grid.size())+".");
623 }
624#endif
625 size_t N[NDIM];
626// size_t counter = 0;
627 sampledvalues_t::iterator destiter = dest_sampled_grid.begin();
628 sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin();
629 if (larger_window == destwindow)
630 std::advance(destiter, pre_offset[0]*total[1]*total[2]);
631 else
632 std::advance(sourceiter, pre_offset[0]*total[1]*total[2]);
633 for(N[0]=0; N[0] < length[0]; ++N[0]) {
634 if (larger_window == destwindow)
635 std::advance(destiter, pre_offset[1]*total[2]);
636 else
637 std::advance(sourceiter, pre_offset[1]*total[2]);
638 for(N[1]=0; N[1] < length[1]; ++N[1]) {
639 if (larger_window == destwindow)
640 std::advance(destiter, pre_offset[2]);
641 else
642 std::advance(sourceiter, pre_offset[2]);
643 for(N[2]=0; N[2] < length[2]; ++N[2]) {
644 ASSERT( destiter != dest_sampled_grid.end(),
645 "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
646 ASSERT( sourceiter != source_sampled_grid.end(),
647 "SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window.");
648 op(*destiter, *sourceiter);
649 ++destiter;
650 ++sourceiter;
651 }
652 if (larger_window == destwindow)
653 std::advance(destiter, post_offset[2]);
654 else
655 std::advance(sourceiter, post_offset[2]);
656 }
657 if (larger_window == destwindow)
658 std::advance(destiter, post_offset[1]*total[2]);
659 else
660 std::advance(sourceiter, post_offset[1]*total[2]);
661 }
662#ifndef NDEBUG
663 if (larger_window == destwindow)
664 std::advance(destiter, post_offset[0]*total[1]*total[2]);
665 else
666 std::advance(sourceiter, post_offset[0]*total[1]*total[2]);
667 ASSERT( destiter == dest_sampled_grid.end(),
668 "SamplingGrid::addWindowOntoWindow() - destiter is not at end of window.");
669 ASSERT( sourceiter == source_sampled_grid.end(),
670 "SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window.");
671#endif
672 LOG(8, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
673}
674
675
676bool SamplingGrid::operator==(const SamplingGrid &other) const
677{
678 bool status =
679 static_cast<const SamplingGridProperties &>(*this)
680 == static_cast<const SamplingGridProperties &>(other);
681 // compare general properties
682 if (status) {
683 // compare windows
684 for (size_t i=0; i<NDIM; ++i) {
685 status &= begin_window[i] == other.begin_window[i];
686 status &= end_window[i] == other.end_window[i];
687 }
688 // compare grids
689 if (status)
690 status &= sampled_grid == other.sampled_grid;
691 }
692 return status;
693}
694
695/** Struct contains a single point with displacements from the
696 * central point and the weight in the restriction.
697 */
698struct PointWeight_t {
699 PointWeight_t(const int &d1, const int &d2, const int &d3, const double &_weight) :
700 displacement(NDIM),
701 weight(_weight)
702 {
703 displacement[0] = d1; displacement[1] = d2; displacement[2] = d3;
704 }
705 typedef std::vector<int> displacement_t;
706 displacement_t displacement;
707 double weight;
708};
709
710static void getLengthsOfWindow(
711 int _total[NDIM],
712 const SamplingGrid &_grid)
713{
714 const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
715 static const double round_offset =
716 (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
717 0.5 : 0.; // need offset to get to round_toward_nearest behavior
718 for (size_t index=0; index<NDIM; ++index) {
719 if (fabs(_grid.end[index] - _grid.begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
720 const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
721 _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
722 } else
723 _total[index] = 0;
724 // we can only assert that its atmost the maximum number of grid points
725 ASSERT (_total[index] <= ::pow(2, _grid.level),
726 "SamplingGrid::downsample() - total "+toString(_total[index])
727 +" is not equal or less than 2^level: "+toString(_grid.level));
728 }
729}
730
731//!> stencil for full weight restriction, see vmg's stencils.hpp
732static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
733 boost::assign::list_of
734 ( PointWeight_t( 0, 0, 0, 0.125) )
735 ( PointWeight_t( 1, 0, 0, 0.0625) )
736 ( PointWeight_t(-1, 0, 0, 0.0625) )
737 ( PointWeight_t( 0, 1, 0, 0.0625) )
738 ( PointWeight_t( 0, -1, 0, 0.0625) )
739 ( PointWeight_t( 0, 0, 1, 0.0625) )
740 ( PointWeight_t( 0, 0, -1, 0.0625) )
741 ( PointWeight_t( 1, 1, 0, 0.03125) )
742 ( PointWeight_t( 1, -1, 0, 0.03125) )
743 ( PointWeight_t(-1, 1, 0, 0.03125) )
744 ( PointWeight_t(-1, -1, 0, 0.03125) )
745 ( PointWeight_t( 0, 1, 1, 0.03125) )
746 ( PointWeight_t( 0, 1, -1, 0.03125) )
747 ( PointWeight_t( 0, -1, 1, 0.03125) )
748 ( PointWeight_t( 0, -1, -1, 0.03125) )
749 ( PointWeight_t( 1, 0, 1, 0.03125) )
750 ( PointWeight_t( 1, 0, -1, 0.03125) )
751 ( PointWeight_t(-1, 0, 1, 0.03125) )
752 ( PointWeight_t(-1, 0, -1, 0.03125) )
753 ( PointWeight_t( 1, 1, 1, 0.015625) )
754 ( PointWeight_t( 1, 1, -1, 0.015625) )
755 ( PointWeight_t( 1, -1, 1, 0.015625) )
756 ( PointWeight_t(-1, 1, 1, 0.015625) )
757 ( PointWeight_t( 1, -1, -1, 0.015625) )
758 ( PointWeight_t(-1, 1, -1, 0.015625) )
759 ( PointWeight_t(-1, -1, 1, 0.015625) )
760 ( PointWeight_t(-1, -1, -1, 0.015625) )
761;
762
763int getValidIndex(
764 const PointWeight_t::displacement_t &_disp,
765 const int N[NDIM],
766 const int length[NDIM])
767{
768 int index = 0;
769 // we simply truncate in case of out of bounds access
770 if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
771 index += _disp[2];
772 if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
773 index += _disp[1]*length[2];
774 if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
775 index += _disp[0]*length[1]*length[2];
776 return index;
777}
778
779void restrictFullWeight(
780 SamplingGrid::sampledvalues_t &_coarse_level,
781 const int length_c[NDIM],
782 const SamplingGrid::sampledvalues_t &_fine_level,
783 const int length_f[NDIM])
784{
785 int N_c[NDIM];
786 int N_f[NDIM];
787 SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
788 for(N_c[0]=0, N_f[0]=0; (N_c[0] < length_c[0]) && (N_f[0] < length_f[0]); ++N_c[0], N_f[0] +=2) {
789 for(N_c[1]=0, N_f[1]=0; (N_c[1] < length_c[1]) && (N_f[1] < length_f[1]); ++N_c[1], N_f[1] +=2) {
790 for(N_c[2]=0, N_f[2]=0; (N_c[2] < length_c[2]) && (N_f[2] < length_f[2]); ++N_c[2], N_f[2] +=2) {
791 const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
792 // go through stencil and add each point relative to displacement with weight
793 for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
794 weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
795 const PointWeight_t::displacement_t disp = weightiter->displacement;
796 const int index_disp = getValidIndex(disp, N_f, length_f);
797 *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
798 }
799 ++coarseiter;
800 }
801 ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
802 "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
803 +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
804 }
805 ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
806 "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
807 +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
808 }
809 ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
810 "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
811 +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
812 ASSERT( coarseiter == _coarse_level.end(),
813 "restrictFullWeight() - coarseiter is not at end of values.");
814}
815
816void SamplingGrid::padWithZerosForEvenNumberedSamples()
817{
818 size_t wbegin_index[NDIM];
819 size_t wend_index[NDIM];
820 size_t wlength_index[NDIM];
821 getDiscreteWindowIndices(wbegin_index, wlength_index, wend_index);
822
823 // calculate new window (always extend it such that both indices are even)
824 bool changed = false;
825 size_t wnewbegin_index[NDIM];
826 size_t wnewend_index[NDIM];
827 for(size_t i=0;i<NDIM;++i) {
828 LOG(2, "INFO: window[" << i << "] starts at " << wbegin_index[i] << " and ends at "
829 << wend_index[i] << " with length " << wlength_index[i]);
830 // We require begin and end of window on even indices (and inside grid).
831 wnewbegin_index[i] = wbegin_index[i];
832 if ((wnewbegin_index[i] % (size_t)2) != 0) {
833 if (wnewbegin_index[i] > 0)
834 --wnewbegin_index[i];
835 else
836 wnewbegin_index[i] = 0;
837 changed = true;
838 }
839 wnewend_index[i] = wend_index[i];
840 if ((wnewend_index[i] % (size_t)2) != 0) {
841 if (wnewend_index[i] < getGridPointsPerAxis())
842 ++wnewend_index[i];
843 else
844 wnewend_index[i] = getGridPointsPerAxis();
845 changed = true;
846 }
847 ASSERT( (wbegin_index[i] >= 0) && (wend_index[i] <= getGridPointsPerAxis()),
848 "SamplingGrid::padWithZerosForEvenNumberedSamples() - indices "
849 +toString(wbegin_index[i])+" and "+toString(wend_index[i])+" larger than grid "
850 +toString(getGridPointsPerAxis())+".");
851 }
852 if (changed) {
853 double begin_newwindow[NDIM];
854 double end_newwindow[NDIM];
855 for(size_t i=0;i<NDIM;++i) {
856 const double delta = getDeltaPerAxis(i);
857 begin_newwindow[i] = begin[i]+delta*wnewbegin_index[i];
858 end_newwindow[i] = begin[i]+delta*wnewend_index[i];
859 }
860 LOG(2, "INFO: Original window is " << Vector(begin_window) << " <-> " << Vector(end_window));
861 LOG(2, "INFO: Padded window is " << Vector(begin_newwindow) << " <-> " << Vector(end_newwindow));
862 // extend window
863 extendWindow(begin_newwindow, end_newwindow);
864 }
865 ASSERT( getWindowGridPoints() % (size_t)8 == 0,
866 "SamplingGrid::padWithZerosForEvenNumberedSamples() - new size "
867 +toString(getWindowGridPoints())+" is still not divisible by 8.");
868}
869
870void SamplingGrid::downsample(
871 SamplingGrid& instance,
872 const SamplingGrid& other,
873 const int _level)
874{
875 if (&instance != &other) {
876 LOG(2, "INFO: other's window is " << Vector(other.begin_window)
877 << " <-> " << Vector(other.end_window));
878 // take over properties
879 static_cast<SamplingGridProperties &>(instance) = other;
880 instance.setWindowSize(other.begin_window, other.end_window);
881 LOG(2, "INFO: Set instance's window to " << Vector(instance.begin_window)
882 << " <-> " << Vector(instance.end_window));
883 ASSERT( _level <= other.level,
884 "SamplingGrid::downsample() - desired level "+toString(_level)
885 +" larger than level "+toString(other.level)+" of the given values.");
886 if (_level == other.level) {
887 instance.sampled_grid = other.sampled_grid;
888 } else {
889 // if desired level is smaller we need to downsample
890 // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
891 // neighbor interpolation) and always one grid level at a time till we
892 // have reached the desired one
893
894 // we need to copy the other grid because we might need to pad it with zeros anyway
895 SamplingGrid FinerGrid(other);
896 int length_d[NDIM];
897 int length_s[NDIM];
898 for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
899 // pad with zeros for even indices and get length of fine grid window
900 FinerGrid.padWithZerosForEvenNumberedSamples();
901 getLengthsOfWindow(length_s, FinerGrid);
902 ASSERT( FinerGrid.getWindowGridPoints() % 8 == 0,
903 "SamplingGrid::downsample() - at level "+toString( instance.level)
904 +" given grid points "+toString(FinerGrid.getWindowGridPoints())+" are not even numbered per axis anymore.");
905 // re-adjust the window (length), get the corresponding window length and downsample
906 instance.setWindow(FinerGrid.begin_window, FinerGrid.end_window);
907 getLengthsOfWindow(length_d, instance);
908 ASSERT( instance.sampled_grid.size() == FinerGrid.getWindowGridPoints()/(size_t)8,
909 "SamplingGrid::downsample() - at level "+toString( instance.level)
910 +" points on coarser grid "+toString(instance.sampled_grid.size())
911 +" and downsampled number on finer grid "
912 +toString(FinerGrid.getWindowGridPoints()/(size_t)8)+" do not match.");
913 restrictFullWeight(instance.sampled_grid, length_d, FinerGrid.sampled_grid, length_s);
914 // then use as new finer grid for next downsampling (if it's not the last level)
915 if (instance.level > _level)
916 FinerGrid = instance;
917 }
918 // loop stops at _level-1
919 instance.level = _level;
920
921 // and finally, renormalize downsampled grid to old value
922// instance *= other.integral()/instance.integral();
923 }
924 }
925}
926
927std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
928{
929 ost << "SamplingGrid";
930 ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
931 ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
932 ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2];
933 ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2];
934 ost << ", level of " << other.level;
935 ost << " and integrated value of " << other.integral();
936 return ost;
937}
938
939template<> SamplingGrid ZeroInstance<SamplingGrid>()
940{
941 SamplingGrid returnvalue;
942 return returnvalue;
943}
944
945// we need to explicitly instantiate the serialization functions as
946// its is only serialized through its base class FragmentJob
947BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)
Note: See TracBrowser for help on using the repository browser.