/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2012 University of Bonn. All rights reserved. * Copyright (C) 2013 Frederik Heber. All rights reserved. * * * This file is part of MoleCuilder. * * MoleCuilder is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * MoleCuilder is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MoleCuilder. If not, see . */ /* * SamplingGrid.cpp * * Created on: 25.07.2012 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif // include headers that implement a archive in simple text format // otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect #include #include #include "CodePatterns/MemDebug.hpp" #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp" #include #include #include #include #include "CodePatterns/Assert.hpp" #include "CodePatterns/Log.hpp" // static instances const double SamplingGrid::zeroOffset[NDIM] = { 0., 0., 0. }; SamplingGrid::SamplingGrid() : SamplingGridProperties() { setWindowSize(zeroOffset, zeroOffset); ASSERT( getWindowGridPoints() == (size_t)0, "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::SamplingGrid(const double _begin[NDIM], const double _end[NDIM], const int _level) : SamplingGridProperties(_begin, _end, _level) { setWindowSize(zeroOffset, zeroOffset); ASSERT( getWindowGridPoints() == (size_t)0, "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::SamplingGrid(const double _begin[NDIM], const double _end[NDIM], const int _level, const sampledvalues_t &_sampled_grid) : SamplingGridProperties(_begin, _end, _level), sampled_grid(_sampled_grid) { setWindowSize(_begin, _end); ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(), "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::SamplingGrid(const SamplingGrid &_grid) : SamplingGridProperties(_grid), sampled_grid(_grid.sampled_grid) { setWindowSize(_grid.begin_window, _grid.end_window); ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(), "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) : SamplingGridProperties(_props) { setWindowSize(zeroOffset, zeroOffset); ASSERT( getWindowGridPoints() == (size_t)0, "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::SamplingGrid( const SamplingGridProperties &_props, const sampledvalues_t &_sampled_grid) : SamplingGridProperties(_props), sampled_grid(_sampled_grid) { setWindowSize(_props.begin, _props.end); ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(), "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window."); } SamplingGrid::~SamplingGrid() {} bool SamplingGrid::isCongruent(const SamplingGrid &_props) const { bool status = true; status &= (static_cast(*this) == static_cast(_props)); for(size_t i = 0; i(*this) = other; setWindowSize(other.begin_window, other.end_window); sampled_grid = other.sampled_grid; } return *this; } static void multiplyElements( double &dest, const double &source, const double prefactor) { dest *= prefactor*(source); } SamplingGrid& SamplingGrid::operator*=(const double _value) { std::transform( sampled_grid.begin(), sampled_grid.end(), sampled_grid.begin(), boost::bind(std::multiplies(), _1, _value)); return *this; } SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other) { // check that grids are compatible if (isEquivalent(other)) { /// get minimum of window double min_begin_window[NDIM]; double min_end_window[NDIM]; bool doShrink = false; for (size_t index=0; index= other.begin_window[index]) { max_begin_window[index] = other.begin_window[index]; doExtend = true; } else { max_begin_window[index] = begin_window[index]; } if (end_window[index] >= other.end_window[index]) { max_end_window[index] = end_window[index]; } else { max_end_window[index] = other.end_window[index]; doExtend = true; } } LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << "."); LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << "."); if (doExtend) extendWindow(max_begin_window, max_end_window); /// and copy other into larger window, too addOntoWindow(other.begin_window, other.end_window, other.sampled_grid, prefactor); } else { ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards."); } } const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const { static const double round_offset( (std::numeric_limits::round_style == std::round_toward_zero) ? 0.5 : 0.); // need offset to get to round_toward_nearest behavior // const double total = getTotalLengthPerAxis(axis); const double delta = getDeltaPerAxis(axis); if (delta == 0) return 0; const double length = getWindowLengthPerAxis(axis); if (length == 0) return 0; return (size_t)(length/delta+round_offset); } double SamplingGrid::integral() const { const double volume_element = getVolume()/(double)getTotalGridPoints(); double int_value = 0.; for (sampledvalues_t::const_iterator iter = sampled_grid.begin(); iter != sampled_grid.end(); ++iter) int_value += *iter; int_value *= volume_element; LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << "."); return int_value; } double SamplingGrid::integral(const SamplingGrid &weight) const { if (isEquivalent(weight)) { const double volume_element = getVolume()/(double)getTotalGridPoints(); double int_value = 0.; sampledvalues_t::const_iterator iter = sampled_grid.begin(); sampledvalues_t::const_iterator weightiter = weight.sampled_grid.begin(); for (;iter != sampled_grid.end();++iter,++weightiter) int_value += (*weightiter) * (*iter); int_value *= volume_element; //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << "."); return int_value; } else return 0.; } void SamplingGrid::setWindowSize( const double _begin_window[NDIM], const double _end_window[NDIM]) { for (size_t index=0;index= begin[index], "SamplingGrid::setWindowSize() - window starts earlier than domain on " +toString(index)+"th component."); end_window[index] = getNearestHigherGridPoint(_end_window[index], index); ASSERT( end_window[index] <= end[index], "SamplingGrid::setWindowSize() - window ends later than domain on " +toString(index)+"th component."); } } void SamplingGrid::setWindow( const double _begin_window[NDIM], const double _end_window[NDIM]) { setWindowSize(_begin_window, _end_window); const size_t gridpoints_window = getWindowGridPoints(); sampled_grid.clear(); sampled_grid.resize(gridpoints_window, 0.); } void SamplingGrid::setDomain( const double _begin[NDIM], const double _end[NDIM]) { setDomainSize(_begin, _end); setWindowSize(_begin, _end); const size_t gridpoints = getTotalGridPoints(); sampled_grid.resize(gridpoints, 0.); } void SamplingGrid::extendWindow( const double _begin_window[NDIM], const double _end_window[NDIM]) { #ifndef NDEBUG for(size_t index=0;index < NDIM; ++index) { // check that we truly have to extend the window ASSERT ( begin_window[index] >= _begin_window[index], "SamplingGrid::extendWindow() - component "+toString(index)+ " of window start is greater than old value."); ASSERT ( end_window[index] <= _end_window[index], "SamplingGrid::extendWindow() - component "+toString(index)+ " of window end is less than old value."); // check that we are still less than domain ASSERT ( _begin_window[index] >= begin[index], "SamplingGrid::extendWindow() - component "+toString(index)+ " of window start is less than domain start."); ASSERT ( _end_window[index] <= end[index], "SamplingGrid::extendWindow() - component "+toString(index)+ " of window end is greater than domain end."); } #endif // copy old window size and values double old_begin_window[NDIM]; double old_end_window[NDIM]; for(size_t index=0;index= _end_window[index], "SamplingGrid::shrinkWindow() - component "+toString(index)+ " of window end is greater than old value."); // check that we are still less than domain ASSERT ( _begin_window[index] >= begin[index], "SamplingGrid::shrinkWindow() - component "+toString(index)+ " of window start is less than domain start."); ASSERT ( _end_window[index] <= end[index], "SamplingGrid::shrinkWindow() - component "+toString(index)+ " of window end is greater than domain end."); } #endif // copy old window size and values double old_begin_window[NDIM]; double old_end_window[NDIM]; for(size_t index=0;index::round_style == std::round_toward_zero) ? 0.5 : 0.; // need offset to get to round_toward_nearest behavior for(size_t index=0;index std::numeric_limits::epsilon()*1e4) { // we refrain from using floor/ceil as the window's starts and ends, // the grids have to be compatible (equal level), should always be on // discrete grid point locations. const double delta = getDeltaPerAxis(index); // delta is conversion factor from box length to discrete length, i.e. number of points _wbegin[index] = (begin_window[index] - begin[index])/delta+round_offset; _wlength[index] = (end_window[index] - begin_window[index])/delta+round_offset; _wend[index] = (end_window[index] - begin[index])/delta+round_offset; } else { _wbegin[index] = 0; _wlength[index] = 0; _wend[index] = 0; } // total is used as safe-guard against loss due to discrete conversion ASSERT( (_wend[index] - _wbegin[index]) == _wlength[index], "SamplingGrid::getDiscreteWindowCopyIndices() - end - begin is not equal to length for " +toString(index)+"th component."); } } void SamplingGrid::getDiscreteWindowOffsets( size_t _pre_offset[NDIM], size_t _post_offset[NDIM], size_t _length[NDIM], size_t _total[NDIM]) const { const double round_offset = (std::numeric_limits::round_style == std::round_toward_zero) ? 0.5 : 0.; // need offset to get to round_toward_nearest behavior for(size_t index=0;index std::numeric_limits::epsilon()*1e4) { // we refrain from using floor/ceil as the window's starts and ends, // the grids have to be compatible (equal level), should always be on // discrete grid point locations. const double delta = getDeltaPerAxis(index); // delta is conversion factor from box length to discrete length, i.e. number of points _pre_offset[index] = (begin_window[index] - begin[index])/delta+round_offset; _post_offset[index] = (end[index] - end_window[index])/delta+round_offset; _length[index] = (end_window[index] - begin_window[index])/delta+round_offset; _total[index] = (end[index] - begin[index])/delta+round_offset; } else { _pre_offset[index] = 0; _post_offset[index] = 0; _length[index] = 0; _total[index] = 0; } // total is used as safe-guard against loss due to discrete conversion ASSERT( (_pre_offset[index] + _post_offset[index]) + _length[index] == _total[index], "SamplingGrid::getDiscreteWindowCopyIndices() - pre, length, post are not equal to total for " +toString(index)+"th component."); } } void SamplingGrid::getDiscreteWindowCopyIndices( const double *larger_wbegin, const double *larger_wend, const double *smaller_wbegin, const double *smaller_wend, size_t *pre_offset, size_t *post_offset, size_t *length, size_t *total) const { const double round_offset = (std::numeric_limits::round_style == std::round_toward_zero) ? 0.5 : 0.; // need offset to get to round_toward_nearest behavior for(size_t index=0;index std::numeric_limits::epsilon()*1e4) { // we refrain from using floor/ceil as the window's starts and ends, // the grids have to be compatible (equal level), should always be on // discrete grid point locations. const double delta = getDeltaPerAxis(index); // delta is conversion factor from box length to discrete length, i.e. number of points pre_offset[index] = (smaller_wbegin[index] - larger_wbegin[index])/delta+round_offset; length[index] = (smaller_wend[index] - smaller_wbegin[index])/delta+round_offset; post_offset[index] = (larger_wend[index] - smaller_wend[index])/delta+round_offset; total[index] = (larger_wend[index] - larger_wbegin[index])/delta+round_offset; } else { pre_offset[index] = 0; length[index] = 0; post_offset[index] = 0; total[index] = 0; } // total is used as safe-guard against loss due to discrete conversion ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index], "SamplingGrid::getDiscreteWindowCopyIndices() - pre, post, and length don't sum up to total for " +toString(index)+"th component."); } } void SamplingGrid::addWindowOntoWindow( const double larger_wbegin[NDIM], const double larger_wend[NDIM], const double smaller_wbegin[NDIM], const double smaller_wend[NDIM], sampledvalues_t &dest_sampled_grid, const sampledvalues_t &source_sampled_grid, boost::function op, enum eLargerWindow larger_window) { #ifndef NDEBUG for(size_t index=0;index= larger_wbegin[index], "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component " +toString(index)+"."); ASSERT( smaller_wend[index] <= larger_wend[index], "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component " +toString(index)+"."); } #endif // the only issue are indices size_t pre_offset[NDIM]; size_t post_offset[NDIM]; size_t length[NDIM]; size_t total[NDIM]; getDiscreteWindowCopyIndices( larger_wbegin, larger_wend, smaller_wbegin, smaller_wend, pre_offset, post_offset, length, total ); // assert that calculated lengths match with given vector sizes #ifndef NDEBUG const size_t calculated_size = length[0]*length[1]*length[2]; if (larger_window == destwindow) { ASSERT( calculated_size == source_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: " +toString(calculated_size)+" != "+toString(source_sampled_grid.size())+"."); ASSERT( calculated_size <= dest_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - not enough sampled values available: " +toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+"."); } else { ASSERT( calculated_size == dest_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: " +toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+"."); ASSERT( calculated_size <= source_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: " +toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+"."); } const size_t total_size = total[0]*total[1]*total[2]; if (larger_window == destwindow) { ASSERT( total_size == dest_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: " +toString(total_size)+" != "+toString(dest_sampled_grid.size())+"."); } else { ASSERT( total_size == source_sampled_grid.size(), "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: " +toString(total_size)+" != "+toString(source_sampled_grid.size())+"."); } #endif size_t N[NDIM]; // size_t counter = 0; sampledvalues_t::iterator destiter = dest_sampled_grid.begin(); sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin(); if (larger_window == destwindow) std::advance(destiter, pre_offset[0]*total[1]*total[2]); else std::advance(sourceiter, pre_offset[0]*total[1]*total[2]); for(N[0]=0; N[0] < length[0]; ++N[0]) { if (larger_window == destwindow) std::advance(destiter, pre_offset[1]*total[2]); else std::advance(sourceiter, pre_offset[1]*total[2]); for(N[1]=0; N[1] < length[1]; ++N[1]) { if (larger_window == destwindow) std::advance(destiter, pre_offset[2]); else std::advance(sourceiter, pre_offset[2]); for(N[2]=0; N[2] < length[2]; ++N[2]) { ASSERT( destiter != dest_sampled_grid.end(), "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window."); ASSERT( sourceiter != source_sampled_grid.end(), "SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window."); op(*destiter, *sourceiter); ++destiter; ++sourceiter; } if (larger_window == destwindow) std::advance(destiter, post_offset[2]); else std::advance(sourceiter, post_offset[2]); } if (larger_window == destwindow) std::advance(destiter, post_offset[1]*total[2]); else std::advance(sourceiter, post_offset[1]*total[2]); } #ifndef NDEBUG if (larger_window == destwindow) std::advance(destiter, post_offset[0]*total[1]*total[2]); else std::advance(sourceiter, post_offset[0]*total[1]*total[2]); ASSERT( destiter == dest_sampled_grid.end(), "SamplingGrid::addWindowOntoWindow() - destiter is not at end of window."); ASSERT( sourceiter == source_sampled_grid.end(), "SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window."); #endif LOG(8, "DEBUG: Grid after adding other is " << dest_sampled_grid << "."); } bool SamplingGrid::operator==(const SamplingGrid &other) const { bool status = static_cast(*this) == static_cast(other); // compare general properties if (status) { // compare windows for (size_t i=0; i displacement_t; displacement_t displacement; double weight; }; static void getLengthsOfWindow( int _total[NDIM], const SamplingGrid &_grid) { const size_t gridpoints_axis = _grid.getGridPointsPerAxis(); static const double round_offset = (std::numeric_limits::round_style == std::round_toward_zero) ? 0.5 : 0.; // need offset to get to round_toward_nearest behavior for (size_t index=0; index std::numeric_limits::epsilon()*1e4) { const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]); _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset; } else _total[index] = 0; // we can only assert that its atmost the maximum number of grid points ASSERT (_total[index] <= ::pow(2, _grid.level), "SamplingGrid::downsample() - total "+toString(_total[index]) +" is not equal or less than 2^level: "+toString(_grid.level)); } } //!> stencil for full weight restriction, see vmg's stencils.hpp static const std::vector< PointWeight_t > FullWeightNearestNeighbor = boost::assign::list_of ( PointWeight_t( 0, 0, 0, 0.125) ) ( PointWeight_t( 1, 0, 0, 0.0625) ) ( PointWeight_t(-1, 0, 0, 0.0625) ) ( PointWeight_t( 0, 1, 0, 0.0625) ) ( PointWeight_t( 0, -1, 0, 0.0625) ) ( PointWeight_t( 0, 0, 1, 0.0625) ) ( PointWeight_t( 0, 0, -1, 0.0625) ) ( PointWeight_t( 1, 1, 0, 0.03125) ) ( PointWeight_t( 1, -1, 0, 0.03125) ) ( PointWeight_t(-1, 1, 0, 0.03125) ) ( PointWeight_t(-1, -1, 0, 0.03125) ) ( PointWeight_t( 0, 1, 1, 0.03125) ) ( PointWeight_t( 0, 1, -1, 0.03125) ) ( PointWeight_t( 0, -1, 1, 0.03125) ) ( PointWeight_t( 0, -1, -1, 0.03125) ) ( PointWeight_t( 1, 0, 1, 0.03125) ) ( PointWeight_t( 1, 0, -1, 0.03125) ) ( PointWeight_t(-1, 0, 1, 0.03125) ) ( PointWeight_t(-1, 0, -1, 0.03125) ) ( PointWeight_t( 1, 1, 1, 0.015625) ) ( PointWeight_t( 1, 1, -1, 0.015625) ) ( PointWeight_t( 1, -1, 1, 0.015625) ) ( PointWeight_t(-1, 1, 1, 0.015625) ) ( PointWeight_t( 1, -1, -1, 0.015625) ) ( PointWeight_t(-1, 1, -1, 0.015625) ) ( PointWeight_t(-1, -1, 1, 0.015625) ) ( PointWeight_t(-1, -1, -1, 0.015625) ) ; int getValidIndex( const PointWeight_t::displacement_t &_disp, const int N[NDIM], const int length[NDIM]) { int index = 0; // we simply truncate in case of out of bounds access if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2])) index += _disp[2]; if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1])) index += _disp[1]*length[2]; if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0])) index += _disp[0]*length[1]*length[2]; return index; } void restrictFullWeight( SamplingGrid::sampledvalues_t &_coarse_level, const int length_c[NDIM], const SamplingGrid::sampledvalues_t &_fine_level, const int length_f[NDIM]) { int N_c[NDIM]; int N_f[NDIM]; SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin(); 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) { 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) { 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) { const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2]; // go through stencil and add each point relative to displacement with weight for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin(); weightiter != FullWeightNearestNeighbor.end(); ++weightiter) { const PointWeight_t::displacement_t disp = weightiter->displacement; const int index_disp = getValidIndex(disp, N_f, length_f); *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight; } ++coarseiter; } ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]), "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2]) +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2])); } ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]), "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1]) +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1])); } ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]), "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0]) +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0])); ASSERT( coarseiter == _coarse_level.end(), "restrictFullWeight() - coarseiter is not at end of values."); } void SamplingGrid::padWithZerosForEvenNumberedSamples() { size_t wbegin_index[NDIM]; size_t wend_index[NDIM]; size_t wlength_index[NDIM]; getDiscreteWindowIndices(wbegin_index, wlength_index, wend_index); // calculate new window (always extend it such that both indices are even) bool changed = false; size_t wnewbegin_index[NDIM]; size_t wnewend_index[NDIM]; for(size_t i=0;i 0) --wnewbegin_index[i]; else wnewbegin_index[i] = 0; changed = true; } wnewend_index[i] = wend_index[i]; if ((wnewend_index[i] % (size_t)2) != 0) { if (wnewend_index[i] < getGridPointsPerAxis()) ++wnewend_index[i]; else wnewend_index[i] = getGridPointsPerAxis(); changed = true; } ASSERT( (wbegin_index[i] >= 0) && (wend_index[i] <= getGridPointsPerAxis()), "SamplingGrid::padWithZerosForEvenNumberedSamples() - indices " +toString(wbegin_index[i])+" and "+toString(wend_index[i])+" larger than grid " +toString(getGridPointsPerAxis())+"."); } if (changed) { double begin_newwindow[NDIM]; double end_newwindow[NDIM]; for(size_t i=0;i(instance) = other; instance.setWindowSize(other.begin_window, other.end_window); ASSERT( _level <= other.level, "SamplingGrid::downsample() - desired level "+toString(_level) +" larger than level "+toString(other.level)+" of the given values."); if (_level == other.level) { instance.sampled_grid = other.sampled_grid; } else { // if desired level is smaller we need to downsample // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest // neighbor interpolation) and always one grid level at a time till we // have reached the desired one // the reference such that we never have to copy the full grid but only // downsampled ones const sampledvalues_t * sourcevalues = &other.sampled_grid; int length_d[3]; int length_s[3]; getLengthsOfWindow(length_s, other); for (instance.level = other.level-1; instance.level >= _level; --instance.level) { getLengthsOfWindow(length_d, instance); // we always have an eighth of the number of sample points as we stop ASSERT( sourcevalues->size() % 8 == 0, "SamplingGrid::downsample() - at level "+toString( instance.level) +" given grid points "+toString(sourcevalues->size())+" are not even numbered per axis anymore."); sampledvalues_t downsampled(sourcevalues->size()/(size_t)8, 0.); restrictFullWeight(downsampled, length_d, *sourcevalues, length_s); // then copy the downsampled values instance.sampled_grid = downsampled; sourcevalues = &instance.sampled_grid; // and exchange lengths for (size_t i=0;i<3;++i) { length_s[i] = length_d[i]; } } instance.level = _level; // and finally, renormalize downsampled grid to old value // instance *= other.integral()/instance.integral(); } } } std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other) { ost << "SamplingGrid"; ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2]; ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2]; ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2]; ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2]; ost << ", level of " << other.level; ost << " and integrated value of " << other.integral(); return ost; } template<> SamplingGrid ZeroInstance() { SamplingGrid returnvalue; return returnvalue; } // we need to explicitly instantiate the serialization functions as // its is only serialized through its base class FragmentJob BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)