Add read progress to GeocodeCov (#761)

gshiroma · GitHub Enterprise · commit b65ea6325994 · 2021-05-03T08:58:28.000-07:00
* add read progress in % to GeocodeCov

* improve getNBlocks() function name to getBlocksNumberAndLength() to be more descriptive and update associated dosctrings

* update getBlocksNumberAndLength() max_block_size type from int to long long

* update areaProjGetNBlocks() parameter max_block_size type from int to long long

* simplify declaration of variables in  getBlocksNumberAndLength()

* remove unnecessary line

* add comments to the GeocodeCov function that reads the RSLC rasters

* substitute c-type cast by static_cast in the GeocodeCov module

* substitute c-type cast by static_cast in the RTC module
diff --git a/cxx/isce3/geocode/GeocodeCov.cpp b/cxx/isce3/geocode/GeocodeCov.cpp
@@ -86,6 +86,57 @@ void Geocode<T>::geoGrid(double geoGridStartX, double geoGridStartY,
     _epsgOut = epsgcode;
 }
 
+
+void getBlocksNumberAndLength(const int array_length, const int array_width,
+                              const int nbands, const size_t type_size,
+                              pyre::journal::info_t* channel, 
+                              int* block_length, int* nblocks_y, 
+                              const long long max_block_size) {
+
+    const int max_block_length = max_block_size / 
+        (nbands * array_width * type_size);
+
+    int _block_length = std::min(array_length, max_block_length);
+
+    // update nblocks
+    int _nblocks_y = std::ceil((static_cast<float>(array_length)) / _block_length);
+
+    if (nblocks_y != nullptr)
+        *nblocks_y = _nblocks_y;
+
+    if (channel != nullptr) {
+        *channel << "array length: " << array_length << pyre::journal::newline;
+        *channel << "array width: " << array_width << pyre::journal::newline;
+        *channel << "number of block(s): " << _nblocks_y
+                 << pyre::journal::newline;
+    }
+
+    if (block_length != nullptr) {
+        *block_length = _block_length;
+        if (channel != nullptr) {
+            *channel << "block length: " << *block_length
+                     << pyre::journal::newline;
+            *channel << "block width: " << array_width
+                     << pyre::journal::newline;
+        }
+    }
+
+    if (channel != nullptr) {
+
+        long long block_size_bytes =
+                ((static_cast<long long>(_block_length)) *
+                 array_width * nbands * type_size);
+
+        std::string block_size_bytes_str = geometry::getNbytesStr(block_size_bytes);
+        if (nbands > 1)
+            block_size_bytes_str += " (" + std::to_string(nbands) + " bands)";
+
+        *channel << "block size: " << block_size_bytes_str
+                 << pyre::journal::endl;
+    }
+}
+
+
 template<class T>
 void Geocode<T>::geocode(
         const isce3::product::RadarGridParameters& radar_grid,
@@ -597,23 +648,91 @@ void _getUpsampledBlock(
 {
     int nbands = input_raster.numBands();
     rdrData.reserve(nbands);
+    bool flag_parallel_radargrid_read = geocode_memory_mode ==
+            geocodeMemoryMode::BLOCKS_GEOGRID_AND_RADARGRID;
     for (int band = 0; band < nbands; ++band) {
-        if (geocode_memory_mode !=
-            geocodeMemoryMode::BLOCKS_GEOGRID_AND_RADARGRID) {
+        if (!flag_parallel_radargrid_read) {
             info << "reading input raster band: " << band + 1
                  << pyre::journal::endl;
         }
+
         rdrData.emplace_back(
                 std::make_unique<isce3::core::Matrix<T_out>>(size_y, size_x));
-        if (!flag_upsample_radar_grid && std::is_same<T, T_out>::value) {
-#pragma omp critical
+
+        if (!flag_upsample_radar_grid && std::is_same<T, T_out>::value &&
+                !flag_parallel_radargrid_read) {
+            /*
+            Enter here if:
+                1. No upsampling is required;
+                2. No type convertion is required (input and output have same
+                   types);
+                3. Not parallel (which allows messages to be printed to stdout).
+            */
+
+            int radargrid_nblocks, radar_block_length;
+            const int max_block_size = 1 << 28; // 256MB
+            isce3::geocode::getBlocksNumberAndLength(size_y, size_x, 1, sizeof(T),
+                   nullptr, &radar_block_length, &radargrid_nblocks,
+                   max_block_size);
+
+            for (size_t block = 0; block < (size_t) radargrid_nblocks;
+                 ++block) {
+
+                int this_radar_block_length = radar_block_length;
+                if ((block + 1) * radar_block_length > size_y) {
+                    this_radar_block_length = size_y % radar_block_length;
+                }
+                if (radargrid_nblocks > 1) {
+                    std::cout << "reading band " << band + 1 << " progress: " 
+                              << static_cast<int>((100.0 * block) / radargrid_nblocks) 
+                              << "% \r";
+                    std::cout.flush();
+                }
+                auto ptr = rdrData[band]->data();
+                input_raster.getBlock(ptr + 
+                                      block * radar_block_length * size_x, 
+                                      xidx, block * radar_block_length + yidx, size_x,
+                                      this_radar_block_length, band + 1);
+            }
+
+            if (radargrid_nblocks > 1) {
+                std::cout << "reading band " << band + 1 
+                          << " progress: 100%" << std::endl;
+            }
+        } 
+        else if (!flag_upsample_radar_grid && std::is_same<T, T_out>::value) {
+            /*
+            Enter here if:
+                1. No upsampling is required;
+                2. No type convertion is required (input and output have same 
+                   types);
+                3. Is parallel (which does not allow messages to be printed to
+                   stdout).
+            */
+            _Pragma("omp critical")
+            {
             input_raster.getBlock(rdrData[band]->data(), xidx, yidx, size_x,
                                   size_y, band + 1);
-        } else if (!flag_upsample_radar_grid) {
+            }    
+        } 
+        else if (!flag_upsample_radar_grid) {
+            /*
+            Enter here if:
+                1. No upsampling is required;
+                2. Type convertion is required (input and output have different 
+                   types).
+            */
             isce3::core::Matrix<T> radar_data_in(size_y, size_x);
-#pragma omp critical
-            input_raster.getBlock(radar_data_in.data(), xidx, yidx, size_x,
-                                  size_y, band + 1);
+            if (flag_parallel_radargrid_read) {
+                _Pragma("omp critical")
+                {
+                input_raster.getBlock(radar_data_in.data(), xidx, yidx, size_x,
+                                       size_y, band + 1);
+                }
+            } else {
+                input_raster.getBlock(radar_data_in.data(), xidx, yidx, size_x,
+                                       size_y, band + 1);
+            }
 
             /*
             Iteratively converts input pixel (ptr_1) to output pixel (ptr_2).
@@ -630,10 +749,20 @@ void _getUpsampledBlock(
                 _convertToOutputType(*ptr_1++, *ptr_2++);
             }
         } else if (flag_upsample_radar_grid && !isce3::is_complex<T>()) {
+            /*
+            Enter here if:
+                1. Upsampling is required;
+                2. Input is not complex.
+            */
             std::string error_msg = "radar-grid upsampling is only available";
             error_msg += " for complex inputs";
             throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
         } else {
+            /*
+            Enter here if:
+                1. Upsampling is required;
+                2. Input is complex.
+            */
             int radargrid_nblocks, radar_block_size;
             if (geocode_memory_mode == geocodeMemoryMode::SINGLE_BLOCK ||
                 geocode_memory_mode ==
diff --git a/cxx/isce3/geocode/GeocodeCov.h b/cxx/isce3/geocode/GeocodeCov.h
@@ -241,6 +241,29 @@ class Geocode {
     void updateGeoGrid(const isce3::product::RadarGridParameters& radar_grid,
                        isce3::io::Raster& dem_raster);
 
+    constexpr static long long DEFAULT_MAX_BLOCK_SIZE = 1 << 29; // 512MB
+    /** Compute the number of blocks and their length for block processing
+     * in the length direction
+     *
+     * @param[in]  array_length        Length of the data to be processed
+     * @param[in]  array_width         Width of the data to be processed
+     * @param[in]  nbands              Number of the bands to be processed
+     * @param[in]  type_size           Type size of the data to be processed
+     * @param[in]  channel             Pyre info channel
+     * @param[out] block_length        Block length
+     * @param[out] nblock_y            Number of blocks in the Y direction
+     * @param[in]  max_block_size      Maximum block size in Bytes (per thread)
+     */
+    void getBlocksNumberAndLength(const int array_length, 
+                                  const int array_width,
+                                  const int nbands = 1,
+                                  const size_t type_size = 4, // Float32
+                                  pyre::journal::info_t* channel = nullptr, 
+                                  int* block_length = nullptr, 
+                                  int* nblock_y = nullptr,
+                                  const long long max_block_size = 
+                                  DEFAULT_MAX_BLOCK_SIZE);
+
     // Get/set data interpolator
     isce3::core::dataInterpMethod dataInterpolator() const 
     { 
diff --git a/cxx/isce3/geometry/RTC.cpp b/cxx/isce3/geometry/RTC.cpp
@@ -426,7 +426,8 @@ void areaProjGetNBlocks(const int array_length, const int array_width,
                         int* block_length_with_upsampling, int* block_length,
                         int* nblocks_y, int* block_width_with_upsampling,
                         int* block_width, int* nblocks_x,
-                        const int min_block_size, const int max_block_size,
+                        const int min_block_size, 
+                        const long long max_block_size,
                         const int nblocks_per_thread)
 {
 
@@ -458,7 +459,7 @@ void areaProjGetNBlocks(const int array_length, const int array_width,
     // set nblocks and block size (Y-axis)
     _nblocks_y = std::max(_nblocks_y, 1);
     _block_length_with_upsampling =
-            std::ceil(((float) array_length) / _nblocks_y);
+            std::ceil((static_cast<float>(array_length)) / _nblocks_y);
     _block_length_with_upsampling =
             std::max(_block_length_with_upsampling, min_block_length);
     _block_length_with_upsampling =
@@ -470,7 +471,7 @@ void areaProjGetNBlocks(const int array_length, const int array_width,
     if (flag_2d) {
         _nblocks_x = std::max(_nblocks_x, 1);
         _block_width_with_upsampling =
-                std::ceil(((float) array_width) / _nblocks_x);
+                std::ceil((static_cast<float>(array_width)) / _nblocks_x);
         _block_width_with_upsampling =
                 std::max(_block_width_with_upsampling, min_block_width);
         _block_width_with_upsampling =
diff --git a/cxx/isce3/geometry/RTC.h b/cxx/isce3/geometry/RTC.h
@@ -38,7 +38,7 @@ enum rtcMemoryMode {
 };
 
 constexpr static int AP_DEFAULT_MIN_BLOCK_SIZE = 1 << 22; // 4MB
-constexpr static int AP_DEFAULT_MAX_BLOCK_SIZE = 1 << 30; // 1GB
+constexpr static long long AP_DEFAULT_MAX_BLOCK_SIZE = 1 << 30; // 1GB
 
 /**Enumeration type to indicate RTC area mode (AREA or AREA_FACTOR) */
 enum rtcAreaMode { AREA = 0, AREA_FACTOR = 1 };
@@ -378,7 +378,7 @@ void areaProjGetNBlocks(const int array_length, const int array_width,
         int* block_width_with_upsampling = nullptr, int* block_width = nullptr,
         int* nblock_x = nullptr,
         const int min_block_size = AP_DEFAULT_MIN_BLOCK_SIZE,
-        const int max_block_size = AP_DEFAULT_MAX_BLOCK_SIZE,
+        const long long max_block_size = AP_DEFAULT_MAX_BLOCK_SIZE,
         const int nblocks_per_thread = 4);
 
 double
