FCostZones.hpp 13.6 KB
Newer Older
1 2 3 4
// ==== CMAKE ====
// Keep in private GIT
// @SCALFMM_PRIVATE

Quentin Khan's avatar
Quentin Khan committed
5 6 7
#ifndef _COSTZONES_HPP_
#define _COSTZONES_HPP_

8
#include "FCostCell.hpp"
9
#include "FCoordColour.hpp"
Quentin Khan's avatar
Quentin Khan committed
10

11
#include <vector>
12 13
#include <stdexcept>
#include <sstream>
14

Quentin Khan's avatar
Quentin Khan committed
15 16 17 18
/**
 * \brief The costzones algorithm implementation.
 * \author Quentin Khan <quentin.khan@inria.fr>
 *
19 20
 * This class is an implementation of the costzones algorithm described in "A
 * Parallel Adaptive Fast Multipole Method" (1993). The algorithm consists in an
Quentin Khan's avatar
Quentin Khan committed
21 22 23
 * in-order traversal of the octree where cell costs are accumulated. When an
 * accumulation is too big, a new zone is created.
 *
24 25 26
 * It is possible to set the levels of the tree that must be considered to
 * compute the costzones.
 *
Quentin Khan's avatar
Quentin Khan committed
27
 * \tparam OctreeClass The type of the octree to work on.
28
 * \tparam CellClass   The type of the cells to work with.
29 30
 * \parblock
 * This class must provide a typedef named CostType that identifies
31
 * the type used to store data.
32
 * \endparblock
Quentin Khan's avatar
Quentin Khan committed
33
 */
Quentin Khan's avatar
Quentin Khan committed
34
template<typename OctreeClass, typename CellClass>
35
class FCostZones {
36
public:
37
    using CostType = typename CellClass::costtype;
38
    using TreeIterator = typename OctreeClass::Iterator;
Quentin Khan's avatar
Quentin Khan committed
39

40 41 42 43 44
    /**
     * \brief Class used to store the bounds of a zone.
     * The bounds consist in the Morton index of the first node and the number
     * of subsequent nodes.
     */ 
45
    using BoundClass = std::pair<TreeIterator, int>;
46
    /// Initial value for empty bounds.
47
    const BoundClass _boundInit {TreeIterator(),0};
48

49

50
protected:
51 52 53 54 55 56 57 58 59 60 61
    /// Alias to FCoordColour#coord2colour -*- not optimised
    constexpr static int(*coord2colour)(const FTreeCoordinate&)
    = FCoordColour::coord2colour;

    /**
     * \brief Enumeration to specify the children to move to during the in-order
     * traversal.
     */
    enum ChildrenSide {LEFT, RIGHT};


62 63 64 65
    /// The iterator to move through the tree.
    typename OctreeClass::Iterator _it;
    /// The number of zones to create.
    int _nbZones;
66
    /// The tree height.
67
    int _treeHeight;
68 69 70 71 72
    /// The highest level in the tree in which to consider cells (inclusive)
    int _topMostLevel = 0;
    /// The lowest level in the tree in which to consider cells (exclusive)
    int _bottomMostLevel = 1;

73

Quentin Khan's avatar
Quentin Khan committed
74
    /// The current cumulative cost of visited cells.
75 76 77 78 79 80 81 82 83 84
    CostType _internalCurrentCost = 0;
    /// The total cost of internal cell.
    CostType _internalTotalCost = 0;
    /**
     * \brief The vector containing the boundaries of the internal node zones
     *
     * Sorted by zone > level > bounds.
     * \details This means there are _treeHeight elements in the inner vectors.
     */
    std::vector< std::vector< BoundClass > > _internalZoneBounds;
85

Quentin Khan's avatar
Quentin Khan committed
86

87 88 89 90
    /// The current cumulative cost of the visited leaves.
    std::vector<CostType> _leafCurrentCost;
    /// The total cost of the leaves.
    std::vector<CostType> _leafTotalCost;
91
    /**
92 93 94 95
     * \brief The vector containing the boundaries of the leaf zones
     *
     * Sorted by zone > colour > bounds.
     * \details This means there are FCoordColour::range elements in the inner
96 97
     * vectors.
     */
98 99 100 101
    std::vector< std::vector< BoundClass > > _leafZoneBounds;

    /// The vector containing the costzone cells sorted by zone > level > cells.
    std::vector< std::vector< std::pair<int, CellClass*> > > _zones;
Quentin Khan's avatar
Quentin Khan committed
102 103


104

Quentin Khan's avatar
Quentin Khan committed
105 106
public:

Quentin Khan's avatar
Quentin Khan committed
107 108 109 110 111
    /**
     * \brief Constructor
     * \param tree    The tree to work on.
     * \param nbZones The number of zones to create.
     */
112
    FCostZones(OctreeClass* tree, int nbZones) : 
113 114 115
        _it( tree ),
        _nbZones( nbZones ),
        _treeHeight( tree->getHeight() ),
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
        _bottomMostLevel(_treeHeight),
        _internalZoneBounds(
            _nbZones,
            std::vector< BoundClass >(_treeHeight, _boundInit )),
        _leafCurrentCost( FCoordColour::range, 0),
        _leafTotalCost( FCoordColour::range, 0),
        _leafZoneBounds(
            _nbZones,
            std::vector< BoundClass >(FCoordColour::range, _boundInit)),
        _zones(
            _nbZones,
            std::vector< std::pair< int, CellClass*> >( ))
        {
            _it.gotoBottomLeft();            
            typename OctreeClass::Iterator ittest(_it);
            ittest.gotoBottomLeft();
        }
Quentin Khan's avatar
Quentin Khan committed
133

Quentin Khan's avatar
Quentin Khan committed
134
    /**
135 136
     * \brief Gets the computed zones.
     *
137
     * See #_zones.
Quentin Khan's avatar
Quentin Khan committed
138 139
     * \return The computed zones.
     */
140
    const std::vector< std::vector< std::pair<int, CellClass*> > >& getZones() const {
Quentin Khan's avatar
Quentin Khan committed
141 142
        return _zones;
    }
Quentin Khan's avatar
Quentin Khan committed
143

Quentin Khan's avatar
Quentin Khan committed
144
    /**
145
     * \brief Gets the computed internal zone bounds.
146
     *
147
     * See #_internalZoneBounds.
148 149
     * \return The computed zone bounds.
     */
Quentin Khan's avatar
Quentin Khan committed
150
    const std::vector< std::vector< BoundClass > >& getZoneBounds() const {
151 152 153 154 155 156 157 158 159 160 161
        return _internalZoneBounds;
    }

    /**
     * \brief Gets the computed leaf zone bounds.
     *
     * See #_leafZoneBounds.
     * \return The computed zone bounds.
     */
    const std::vector< std::vector< BoundClass > >& getLeafZoneBounds() const {
        return _leafZoneBounds;
162 163
    }

164
    /// Gets the tree topmost level used. 
165 166 167 168
    int getTopMostLevel() const {
        return _topMostLevel;
    }

169
    /// Sets the tree topmost level used. 
170 171 172 173 174
    void setTopMostLevel(unsigned int level) {
        if( level > _treeHeight-1 ) {
            std::stringstream msgstream;
            msgstream << __FUNCTION__ << ": level is to deep. level=" << level
                      << " tree height=" << _treeHeight;
175
            throw std::out_of_range(msgstream.str());
176 177 178 179 180
        }
            
        _topMostLevel = level;
    }

181
    /// Gets the tree bottom most level that we use. 
182 183 184 185
    int getBottomMostLevel() const {
        return _bottomMostLevel;
    }

186
    /// Sets the tree bottom most level that we use. 
187 188 189 190 191
    void setBottomMostLevel(unsigned int level) {
        if( level > _treeHeight-1 ) {
            std::stringstream msgstream;
            msgstream << __FUNCTION__ << ": level is to deep. level=" << level
                      << " tree height=" << _treeHeight;
192
            throw std::out_of_range(msgstream.str());
193 194 195 196
        }
            
        _bottomMostLevel = level;
    }
197

198

199 200
    /**
     * \brief Runs the costzones algorithm.
201 202 203
     *
     * Since ScalFMM's implementation of the octree doesn't have a real root, we
     * run the algorithm on each of its children.
Quentin Khan's avatar
Quentin Khan committed
204
     */
Quentin Khan's avatar
Quentin Khan committed
205
    void run() {
206 207 208 209 210 211 212 213

        // Compute tree leaves total cost;
        computeLeavesCost();
        // Compute tree internal nodes total cost;
        computeInternalCost();

        // Count the root's children (the root is not stored in the tree)
        _it.gotoTop();
Quentin Khan's avatar
Quentin Khan committed
214
        _it.gotoLeft();
215
        do {
216
            this->costzones();
217
        } while( _it.moveRight() );
218

Quentin Khan's avatar
Quentin Khan committed
219 220
    }

221
protected:
Quentin Khan's avatar
Quentin Khan committed
222 223

    /**
224
     * \brief Main costzone algorithm.
225
     *
226
     * Moves through the tree in-order and assigns each cell to a zone. When a
227
     * zone's cumulative cost is too high, the new cells are inserted in the
228
     * next one.
Quentin Khan's avatar
Quentin Khan committed
229
     */
230 231 232 233 234 235
    void costzones() {
        
        std::pair<int,int> childrenCount;
        const int level = _it.level();
        const bool progressDown = _it.canProgressToDown()
            && (level < _bottomMostLevel);
236
        // Is this cell within the level range we consider
237 238 239 240 241 242 243
        const bool useCell = (level < _bottomMostLevel)
            && (level >= _topMostLevel);

        // When not on a leaf, apply to left children first
        if ( progressDown ) {
            childrenCount = countLeftRightChildren();
            callCostZonesOnChildren(LEFT, childrenCount);
Quentin Khan's avatar
Quentin Khan committed
244
        }
245
        
246
        // if the current cell is within the levels we consider, we add it
247 248 249 250 251 252
        if( useCell )
            addCurrentCell();
        
        // When not on a leaf, apply to right children
        if ( progressDown ) {
            callCostZonesOnChildren(RIGHT, childrenCount);
253 254
        }

Quentin Khan's avatar
Quentin Khan committed
255 256
    }

257
    
258
    /**
Quentin Khan's avatar
Quentin Khan committed
259 260 261 262 263
     * \brief Applies costzones to the left or right children of the current cell.
     *
     * The current cell is the one currently pointed at by the iterator _it.
     *
     * \warning You must check by yourself whether the cell is a leaf or not.
264 265 266
     *
     * \param side The children side we want to visit.
     * \param childrenCount The children count as returned by
Quentin Khan's avatar
Quentin Khan committed
267
     *                      countLeftRightChildren.
268 269
     */
    void callCostZonesOnChildren(const ChildrenSide side, const std::pair<int, int>& childrenCount) {
Quentin Khan's avatar
Quentin Khan committed
270
        
271
        const int& nbChildren = (side == LEFT ? childrenCount.first : childrenCount.second);
Quentin Khan's avatar
Quentin Khan committed
272

273 274 275 276 277
        // Don't move if there are no children on the right when we want to
        // visit them. We test this before moving in case one day moving in the
        // tree becomes expensive.
        if ( side == RIGHT && childrenCount.second == 0)
            return;
Quentin Khan's avatar
Quentin Khan committed
278

279
        // move down to the children level
Quentin Khan's avatar
Quentin Khan committed
280 281 282
        _it.moveDown();
        
        if ( side == RIGHT ) {
283 284 285 286
            // move to the first right child
            for ( int childIdx = 0; childIdx < childrenCount.first; childIdx++) {
                _it.moveRight();
            }
Quentin Khan's avatar
Quentin Khan committed
287 288
        }

289 290
        // Call costzones
        for ( int childIdx = 0; childIdx < nbChildren; childIdx++ ) {
291
            this->costzones();
292 293 294
            if(childIdx < nbChildren -1) // nbChildren-1 to avoid changing tree
                _it.moveRight();
        }
Quentin Khan's avatar
Quentin Khan committed
295

296
        // move up to the previous cell level
Quentin Khan's avatar
Quentin Khan committed
297
        _it.moveUp();
298

Quentin Khan's avatar
Quentin Khan committed
299 300
    }

301 302 303 304

    /**
     * \brief Adds the current cell to a zone.
     *
305 306 307 308 309 310 311 312 313 314 315
     * The choice of the zone is made according to the current cost accumulation
     * compared to the mean cost of a zone (_totalCost/_nbZones +1).
     *
     * This method uses the following attributes to choose the zone into which
     * the current cell must be stored :
     *
     *   - #_internalCurrentCost
     *   - #_leafCurrentCost
     *   - #_internalTotalCost
     *   - #_leafTotalCost
     *   - #_nbZones
316
     */
317
    void addCurrentCell() {
318 319

        const int& level = _it.level();
320
        CellClass* cell = _it.getCurrentCell();
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
        CostType cellCost = cell->getCost();
        bool isLeaf = (level == _treeHeight -1);

        if ( 0 == cellCost ) {
            return;
        }

        // find cell zone
        long long int cellZone = 0;

        // Near-field zone /////////////////
        if( isLeaf ) {
            CostType leafCost = cell->getNearCost();
            int colour = coord2colour(_it.getCurrentGlobalCoordinate());

            cellZone = _leafCurrentCost[colour] * _nbZones / (_leafTotalCost[colour]+1);
            _leafCurrentCost[colour] += leafCost;

            if( _leafZoneBounds.at(cellZone)[colour] == _boundInit ) {
340
                _leafZoneBounds.at(cellZone)[colour].first = _it;
341
                _leafZoneBounds.at(cellZone)[colour].second = 1;
342
            } else {
343
                _leafZoneBounds.at(cellZone)[colour].second++;
344
            }
345 346 347 348 349 350 351 352 353
        }
        ////////////////////////////////////

        // Far-field zone //////////////////
        cellZone = _internalCurrentCost * _nbZones / (_internalTotalCost+1);

        _internalCurrentCost += cellCost;

        if( _boundInit == _internalZoneBounds.at(cellZone)[level] ) {
354
            _internalZoneBounds.at(cellZone)[level].first = _it;
355 356 357 358 359 360 361 362
            _internalZoneBounds.at(cellZone)[level].second = 1;
        } else {
            _internalZoneBounds.at(cellZone)[level].second++;                
        }
        ///////////////////////////////

        // add cell to exhaustive zone vector
        (_zones.at(cellZone)).emplace_back(level, cell);
363
    }
Quentin Khan's avatar
Quentin Khan committed
364

Quentin Khan's avatar
Quentin Khan committed
365 366

    /**
367
     * \brief Computes and stores the leaves' total cost.
Quentin Khan's avatar
Quentin Khan committed
368
     *
369
     * The tree iterator (#_it) is moved to the bottom level of the
370 371
     * tree by this method. After the method returns, the iterator is left at
     * the rightmost leaf.
Quentin Khan's avatar
Quentin Khan committed
372
     */
373
    void computeLeavesCost() {
Quentin Khan's avatar
Quentin Khan committed
374

375 376 377
        // Reset colour costs.
        for( CostType& colourCost : _leafTotalCost ) {
            colourCost = 0;
378 379
        }

380 381 382 383 384
        _it.gotoBottomLeft();
        do {
            int leafColour = coord2colour(_it.getCurrentGlobalCoordinate());
            _leafTotalCost[leafColour] += _it.getCurrentCell()->getNearCost();
        } while(_it.moveRight());        
385

386
    }
Quentin Khan's avatar
Quentin Khan committed
387

388 389 390 391
    /**
     * \brief Computes and stores the internal cells' total cost.
     * \warning This method makes use of 
     *
392
     * The tree iterator (#_it) is moved to the bottom level of the
393 394 395 396 397 398 399 400 401
     * tree by this method. After the method returns, the iterator is left at
     * the rightmost leaf.
     */
    void computeInternalCost() {
        _it.gotoBottomLeft();
        //_it.moveUp();

        while( _it.level() >= _bottomMostLevel ) {
            _it.moveUp();
402 403
        }

404 405 406 407 408 409
        do {
            _it.gotoLeft();
            do {
                _internalTotalCost += _it.getCurrentCell()->getCost();
            } while(_it.moveRight());        
        } while(_it.moveUp());
410

411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
    }

    /**
     * \brief Counts the left and right children of the current cell.
     *
     * The current cell is the one currently pointed at by the iterator _it.
     *
     * \warning It must be checked whether the current cell is a leaf or not
     * before calling this method.
     *
     * \return A pair of int containing the count of left (first) and right
     *         (second) children.
     */
    std::pair<int,int> countLeftRightChildren() {
        CellClass** children = _it.getCurrentChildren();
        int nbLeftChildren = 0, nbRightChildren = 0; 
        // Left children
        for ( int childIdx = 0; childIdx < 4; childIdx++) {
            if ( children[childIdx] != nullptr ) {
                ++ nbLeftChildren;
            }
Quentin Khan's avatar
Quentin Khan committed
432
        }
433 434 435 436 437 438 439 440
        // Right children
        for ( int childIdx = 4; childIdx < 8; childIdx++) {
            if ( children[childIdx] != nullptr) {
                ++ nbRightChildren;
            }
        }

        return std::pair<int,int> (nbLeftChildren, nbRightChildren);
441
    }
Quentin Khan's avatar
Quentin Khan committed
442 443 444



445 446
};

Quentin Khan's avatar
Quentin Khan committed
447
#endif