From c5a183ada30fbd1b8587121e0d313b4cda9390c0 Mon Sep 17 00:00:00 2001 From: Richard Barnes Date: Sun, 2 Jul 2023 17:56:18 -0400 Subject: [PATCH] Adds CachingOutletChecker. In tests on the Steele County 3m dataset this decreases wall-time for marginal volume calculations from ~14s to ~7s. --- .../depressions/depression_hierarchy.hpp | 86 ++++++++++++++++++- 1 file changed, 82 insertions(+), 4 deletions(-) diff --git a/include/richdem/depressions/depression_hierarchy.hpp b/include/richdem/depressions/depression_hierarchy.hpp index 72e8d129..101fc9db 100644 --- a/include/richdem/depressions/depression_hierarchy.hpp +++ b/include/richdem/depressions/depression_hierarchy.hpp @@ -724,6 +724,71 @@ GetDepressionHierarchy(const Array2D& dem, Array2D& label, A return depressions; } +// Accelerates marginal volume calculations by caching previous DH lookups +template +struct CachingOutletChecker { + // The previous cell we looked at was part of this leaf depression + dh_label_t previous_leaf_label = NO_DEP; + + // This is a record of the elevation of each outlet we passed and the + // metadepression associated with it during the previous cell's marginal + // volume calculation. If the cell we're currently looking at was part of the + // same leaf depression as the previous cell then the same set of outlets will + // apply. We cache them here to avoid having to revisit them in the DH. + // + // This is faster for two reasons: + // 1. This outlets vector is more cache friendly. + // 2. We can start our search for an appropriate outlet deeper in the DH; we + // don't have to start back at the leaf depression. + std::vector> outlets; + + void reset(const dh_label_t leaf_label){ + previous_leaf_label = leaf_label; + // Clearing allows us to reuse memory, which is faster than rebuilding the + // object + outlets.clear(); + } + + /// Remember that this metadepression had an outlet at this elevation + void operator()(const dh_label_t this_label, const elev_t this_elev){ + outlets.emplace_back(this_elev, this_label); + } + + // But it contributed to the marginal volume of this depression + // dh_label_t previous_mv_label = NO_DEP; + // This is the outlet just below the MV depression. Maybe the next cell will + // be contained in it if it's lower than the previous cell. + // elev_t previous_lower_outlet = std::numeric_limits::max(); + // This is the outlet of the MV depression. If the next cell is higher than + // the outlet we need to keep looking. + // elev_t previous_upper_outlet = std::numeric_limits::max(); + + + /// Which metadepression is the best one to start looking for an outlet in + /// given that the cell we're currently looking at has this \p leaf_label and + /// \p this_elev ? + dh_label_t find_start(const dh_label_t leaf_label, const elev_t this_elev){ + // Cells are part of different leaf depressions, so we have to search the + // whole DH to find where the marginal volume (MV) should go + if(previous_leaf_label != leaf_label){ + reset(leaf_label); + return leaf_label; + } + + // Search until we find a depression that contains this cell + while (!outlets.empty() && this_elev <= outlets.back().first) { + outlets.pop_back(); + } + + if(!outlets.empty()){ + return outlets.back().second; + } + + reset(leaf_label); + return leaf_label; + } +}; + template void CalculateMarginalVolumes( DepressionHierarchy& deps, @@ -739,6 +804,7 @@ void CalculateMarginalVolumes( { std::vector cell_counts(deps.size(), 0); std::vector total_elevations(deps.size(), 0); + CachingOutletChecker coc; #pragma omp for for (unsigned int i = 0; i < dem.size(); i++) { @@ -746,6 +812,10 @@ void CalculateMarginalVolumes( const auto my_elev = dem(i); auto clabel = label(i); + // Choose a good place in the DH to start looking for our containing + // metadepression + clabel = coc.find_start(clabel, my_elev); + // This cell contributes volume to some depression or meta-depression. // Let's find out which. To do so, we start with the leaf depression the // cell is part of. If it's not contained by that depression, we look at @@ -754,11 +824,19 @@ void CalculateMarginalVolumes( // elevation. If at any point we pass through an oceanlink, our search // terminates: any water falling on the cell will ultimately end up in the // ocean and the cell is not contained by any depression. - while (clabel != OCEAN && my_elev > deps.at(clabel).out_elev) { + while (clabel != OCEAN) { + const auto this_out_elev = deps.at(clabel).out_elev; + coc(clabel, this_out_elev); + + // Stop: we found a depression that contains us + if (my_elev <= this_out_elev){ + break; + } + + // There is no meta-depression: the next depression in the hierarchy is + // downhill of the one we started in. Therefore, our marginal volume is + // contributed to the ocean. if (deps[clabel].ocean_parent){ - // There is no meta-depression: the next depression in the hierarchy - // is downhill of the one we started in. Therefore, our marginal - // volume is contributed to the ocean. clabel = OCEAN; break; }