From 9c898b504a8149b8999d1a2709d77f13a850aec8 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 19 Mar 2026 16:09:43 -0700 Subject: [PATCH] feat: add configurable spread neighborhood radius (SpreadRad) --- .github/workflows/build-debian-stable.yml | 7 +- Cell2Fire/Cell2Fire.cpp | 15 +-- Cell2Fire/Cells.cpp | 154 ++++++++++++---------- Cell2Fire/Cells.h | 11 +- Cell2Fire/ReadArgs.cpp | 16 +++ Cell2Fire/ReadArgs.h | 2 +- README.md | 10 ++ test/spreadrad_test.sh | 69 ++++++++++ 8 files changed, 202 insertions(+), 82 deletions(-) create mode 100644 test/spreadrad_test.sh diff --git a/.github/workflows/build-debian-stable.yml b/.github/workflows/build-debian-stable.yml index eb39d89..f3f00c9 100644 --- a/.github/workflows/build-debian-stable.yml +++ b/.github/workflows/build-debian-stable.yml @@ -115,4 +115,9 @@ jobs: - name: Bash Test run: | cd test - bash test.sh ${{ env.SUFFIX }} \ No newline at end of file + bash test.sh ${{ env.SUFFIX }} + + - name: SpreadRad Test + run: | + cd test + bash spreadrad_test.sh ${{ env.SUFFIX }} diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index f0a66a2..ffc2c08 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -597,7 +597,7 @@ Cell2Fire::InitCell(int id) it2 = this->Cells_Obj.find(id); // Initialize the fire fields for the selected cel - it2->second.initializeFireFields(this->coordCells, this->availCells, this->cols, this->rows); + it2->second.initializeFireFields(this->coordCells, this->availCells, this->cols, this->rows, this->args.SpreadRadius); // Print info for debugging if (this->args.verbose) @@ -1695,17 +1695,16 @@ Cell2Fire::GetMessages(std::unordered_map> sendMessageList // Cleaning step // Fire can't be propagated back - int cellNum = it->second.realId - 1; - for (auto& angle : it->second.angleToNb) + for (auto& nbAndAngle : it->second.angleDict) { - int origToNew = angle.first; - // Which neighbor am I to the burnt cell - int newToOrig = (origToNew + 180) % 360; - int adjCellNum = angle.second; // Check + int adjCellNum = nbAndAngle.first; auto adjIt = Cells_Obj.find(adjCellNum); if (adjIt != Cells_Obj.end()) { - adjIt->second.ROSAngleDir.erase(newToOrig); + adjIt->second.ROSAngleDir.erase(it->second.realId); + adjIt->second.fireProgress.erase(it->second.realId); + adjIt->second.angleDict.erase(it->second.realId); + adjIt->second.distToCenter.erase(it->second.realId); } } } diff --git a/Cell2Fire/Cells.cpp b/Cell2Fire/Cells.cpp index 9a3e7ca..16ba399 100644 --- a/Cell2Fire/Cells.cpp +++ b/Cell2Fire/Cells.cpp @@ -96,7 +96,6 @@ Cells::Cells(int _id, this->angleDict = std::unordered_map(); this->ROSAngleDir = std::unordered_map(); this->distToCenter = std::unordered_map(); - this->angleToNb = std::unordered_map(); } /** @@ -120,54 +119,34 @@ void Cells::initializeFireFields(std::vector>& coordCells, // TODO: should probably make a coordinate type std::unordered_set& availSet, int cols, - int rows) // WORKING CHECK OK + int rows, + int spreadRadius) // WORKING CHECK OK { - std::vector adj = adjacentCells(this->realId, rows, cols); + this->angleDict.clear(); + this->ROSAngleDir.clear(); + this->distToCenter.clear(); + this->fireProgress.clear(); + + std::vector adj = adjacentCells(this->realId, rows, cols, spreadRadius); + const double radToDeg = 180.0 / M_PI; for (auto& nb : adj) { - // CP Default value is replaced: None = -1 - // std::cout << "DEBUG1: adjacent: " << nb.second << std::endl; if (nb != -1) { int a = -1 * coordCells[nb - 1][0] + coordCells[this->id][0]; int b = -1 * coordCells[nb - 1][1] + coordCells[this->id][1]; - int angle = -1; - if (a == 0) + double angle = std::atan2(-1.0 * b, -1.0 * a) * radToDeg; + if (angle < 0.0) { - if (b >= 0) - angle = 270; - else - angle = 90; - } - else if (b == 0) - { - if (a >= 0) - angle = 180; - else - angle = 0; - } - else - { - // TODO: check this logi - double radToDeg = 180 / M_PI; - // TODO: i think all the negatives and abs cancel out - double temp = std::atan(b * 1.0 / a) * radToDeg; - if (a > 0) - temp += 180; - if (a < 0 && b > 0) - temp += 360; - angle = temp; + angle += 360.0; } this->angleDict[nb] = angle; if (availSet.find(nb) != availSet.end()) { - // TODO: cannot be None, replaced None = -1 and ROSAngleDir has - // a double inside - this->ROSAngleDir[angle] = -1; + this->ROSAngleDir[nb] = -1; } - this->angleToNb[angle] = nb; this->fireProgress[nb] = 0.0; this->distToCenter[nb] = std::sqrt(a * a + b * b) * this->_ctr2ctrdist; } @@ -192,26 +171,74 @@ Cells::initializeFireFields(std::vector>& coordCells, // TODO: * -1. */ std::vector -adjacentCells(int cell, int nrows, int ncols) +adjacentCells(int cell, int nrows, int ncols, int radius) { if (cell <= 0 || cell > nrows * ncols) { - std::vector adjacents(8, -1); + std::vector adjacents; return adjacents; } - int total_cells = nrows * ncols; - int north = cell <= ncols ? -1 : cell - ncols; - int south = cell + ncols > total_cells ? -1 : cell + ncols; - int east = cell % ncols == 0 ? -1 : cell + 1; - int west = cell % ncols == 1 ? -1 : cell - 1; - int northeast = cell < ncols || cell % ncols == 0 ? -1 : cell - ncols + 1; - int southeast = cell + ncols > total_cells || cell % ncols == 0 ? -1 : cell + ncols + 1; - int southwest = cell % ncols == 1 || cell + ncols > total_cells ? -1 : cell + ncols - 1; - int northwest = cell % ncols == 1 || cell < ncols ? -1 : cell - ncols - 1; - std::vector adjacents = { west, east, southwest, southeast, south, northwest, northeast, north }; + if (radius < 1) + { + radius = 1; + } + int row = (cell - 1) / ncols; + int col = (cell - 1) % ncols; + std::vector adjacents; + for (int dr = -radius; dr <= radius; ++dr) + { + for (int dc = -radius; dc <= radius; ++dc) + { + if (dr == 0 && dc == 0) + { + continue; + } + int nr = row + dr; + int nc = col + dc; + if (nr < 0 || nr >= nrows || nc < 0 || nc >= ncols) + { + continue; + } + adjacents.push_back(nr * ncols + nc + 1); + } + } return adjacents; } +int +Cells::selectHeadCell(double windAzimuth) const +{ + if (this->angleDict.empty()) + { + return this->realId; + } + double target = std::fmod(windAzimuth, 360.0); + if (target < 0.0) + { + target += 360.0; + } + + int bestNeighbor = this->realId; + double bestAngleDiff = 1e9; + double bestDistance = 1e18; + for (const auto& nbAndAngle : this->angleDict) + { + const int nb = nbAndAngle.first; + const double angle = nbAndAngle.second; + double diff = std::fabs(angle - target); + diff = std::min(diff, 360.0 - diff); + double dist = this->distToCenter.count(nb) ? this->distToCenter.at(nb) : 1e18; + + if (diff < bestAngleDiff || (std::fabs(diff - bestAngleDiff) < 1e-9 && dist < bestDistance)) + { + bestAngleDiff = diff; + bestDistance = dist; + bestNeighbor = nb; + } + } + return bestNeighbor; +} + /** * @brief Calculates the radial distance for a given angle in an ellipse * defined by its semi-major and semi-minor axes. @@ -265,9 +292,9 @@ Cells::ros_distr_V2(double thetafire, double a, double b, double c, double EFact { // Ros allocation for each angle inside the dictionary - for (auto& angle : this->ROSAngleDir) + for (auto& nbRos : this->ROSAngleDir) { - double offset = angle.first - thetafire; + double offset = this->angleDict[nbRos.first] - thetafire; if (offset < 0) { @@ -277,7 +304,7 @@ Cells::ros_distr_V2(double thetafire, double a, double b, double c, double EFact { offset -= 360; } - this->ROSAngleDir[angle.first] = rhoTheta(offset, a, b) * EFactor; + this->ROSAngleDir[nbRos.first] = rhoTheta(offset, a, b) * EFactor; } } @@ -393,7 +420,7 @@ Cells::manageFire(int period, df_ptr[this->realId - 1].bui = wdf_ptr->bui; df_ptr[this->realId - 1].ffmc = wdf_ptr->ffmc; - int head_cell = angleToNb[wdf_ptr->waz]; // head cell for slope calculation + int head_cell = selectHeadCell(wdf_ptr->waz); // head cell for slope calculation if (head_cell <= 0) // solve boundaries case { head_cell = this->realId; // as it is used only for slope calculation, if @@ -527,14 +554,15 @@ Cells::manageFire(int period, args->EFactor); // std::cout << "Sale de Ros Dist" << std::endl; - // Fire progress using ROS from burning cell, not the neighbors // + // Fire progress uses the source-cell ROS for all outgoing arcs. + // This intentionally retains the source-cell spread assumption. // vector toPop = vector(); // this is a iterator through the keyset of a dictionary for (auto& _angle : this->ROSAngleDir) { - double angle = _angle.first; - int nb = angleToNb[angle]; + int nb = _angle.first; + double angle = this->angleDict[nb]; double ros = (1 + args->ROSCV * ROSRV) * _angle.second; if (std::isnan(ros)) @@ -699,7 +727,7 @@ Cells::manageFireBBO(int period, df_ptr->ffmc = wdf_ptr->ffmc; df_ptr->tmp = wdf_ptr->tmp; df_ptr->rh = wdf_ptr->rh; - int head_cell = angleToNb[wdf_ptr->waz]; // head cell for slope calculation + int head_cell = selectHeadCell(wdf_ptr->waz); // head cell for slope calculation if (head_cell <= 0) // solve boundaries case { head_cell = this->realId; // as it is used only for slope calculation, if @@ -821,14 +849,15 @@ Cells::manageFireBBO(int period, EllipseFactors[3]); // std::cout << "Sale de Ros Dist" << std::endl; - // Fire progress using ROS from burning cell, not the neighbors // + // Fire progress uses the source-cell ROS for all outgoing arcs. + // This intentionally retains the source-cell spread assumption. // vector toPop = vector(); // this is a iterator through the keyset of a dictionary for (auto& _angle : this->ROSAngleDir) { - double angle = _angle.first; - int nb = angleToNb[angle]; + int nb = _angle.first; + double angle = this->angleDict[nb]; double ros = (1 + args->ROSCV * ROSRV) * _angle.second; if (args->verbose) @@ -969,7 +998,7 @@ Cells::get_burned(int period, df[this->id].ws = wdf_ptr->ws; df[this->id].tmp = wdf_ptr->tmp; df[this->id].rh = wdf_ptr->rh; - int head_cell = angleToNb[wdf_ptr->waz]; // head cell for slope calculation + int head_cell = selectHeadCell(wdf_ptr->waz); // head cell for slope calculation if (head_cell <= 0) // solve boundaries case { head_cell = this->realId; // as it is used only for slope calculation, if @@ -1108,7 +1137,7 @@ Cells::ignition(int period, df_ptr->ws = wdf_ptr->ws; df_ptr->bui = wdf_ptr->bui; df_ptr->ffmc = wdf_ptr->ffmc; - int head_cell = angleToNb[wdf_ptr->waz]; // head cell for slope calculation + int head_cell = selectHeadCell(wdf_ptr->waz); // head cell for slope calculation if (head_cell <= 0) // solve boundaries case { head_cell = this->realId; // as it is used only for slope calculation, if @@ -1206,13 +1235,6 @@ Cells::print_info() } std::cout << std::endl; - printf("angleToNb Dict: "); - for (auto& nb : this->angleToNb) - { - std::cout << " " << nb.first << " : " << nb.second; - } - std::cout << std::endl; - printf("fireProgress Dict: "); for (auto& nb : this->fireProgress) { diff --git a/Cell2Fire/Cells.h b/Cell2Fire/Cells.h index f443cc5..44dce82 100644 --- a/Cell2Fire/Cells.h +++ b/Cell2Fire/Cells.h @@ -13,7 +13,7 @@ using namespace std; -std::vector adjacentCells(int cell, int nrows, int ncols); +std::vector adjacentCells(int cell, int nrows, int ncols, int radius = 1); /* * Weather structure */ @@ -96,11 +96,8 @@ class Cells std::unordered_map> gMsgListSeason; std::unordered_map fireProgress; // CP: dictionary {int: double} std::unordered_map angleDict; // CP: dictionary {int: double} - std::unordered_map ROSAngleDir; // CP: dictionary {int: double|None} Instead of None we - // can use a determined number like -9999 = None TODO: - // maybe int : double + std::unordered_map ROSAngleDir; // dictionary {neighbor id: ROS from this source cell} std::unordered_map distToCenter; // CP: dictionary {int: double} - std::unordered_map angleToNb; // CP: dictionary {double: int} // TODO: reference to shared object @@ -117,7 +114,8 @@ class Cells void initializeFireFields(std::vector>& coordCells, std::unordered_set& availSet, int cols, - int rows); + int rows, + int spreadRadius); double rhoTheta(double theta, double a, double b); void ros_distr_V2(double thetafire, double a, double b, double c, double EFactor); @@ -198,6 +196,7 @@ class Cells private: double allocate(double offset, double base, double ros1, double ros2); float slope_effect(float elev_i, float elev_j, int cellsize); + int selectHeadCell(double windAzimuth) const; }; #endif diff --git a/Cell2Fire/ReadArgs.cpp b/Cell2Fire/ReadArgs.cpp index 52ec65e..7f559cb 100644 --- a/Cell2Fire/ReadArgs.cpp +++ b/Cell2Fire/ReadArgs.cpp @@ -231,6 +231,7 @@ parseArgs(int argc, char* argv[], arguments* args_ptr) int dmax_fire_periods = 1000; int dseed = 123; int diradius = 0; + int dspreadradius = 1; int dnthreads = 1; int dfmc = 100; float dROS_Threshold = 0.1; @@ -329,6 +330,20 @@ parseArgs(int argc, char* argv[], arguments* args_ptr) else args_ptr->IgnitionRadius = diradius; + //--SpreadRad + char* input_sprad = getCmdOption(argv, argv + argc, "--SpreadRad"); + if (input_sprad) + { + printf("SpreadRadius: %s \n", input_sprad); + args_ptr->SpreadRadius = std::stoi(input_sprad, &sz); + } + else + args_ptr->SpreadRadius = dspreadradius; + if (args_ptr->SpreadRadius < 1) + { + args_ptr->SpreadRadius = 1; + } + //--fmc char* input_fmc = getCmdOption(argv, argv + argc, "--fmc"); if (input_fmc) @@ -588,6 +603,7 @@ printArgs(arguments args) std::cout << "FirePeriodLen: " << args.FirePeriodLen << std::endl; std::cout << "Ignitions: " << args.Ignitions << std::endl; std::cout << "IgnitionRad: " << args.IgnitionRadius << std::endl; + std::cout << "SpreadRad: " << args.SpreadRadius << std::endl; std::cout << "OutputGrid: " << args.OutputGrids << std::endl; std::cout << "FinalGrid: " << args.FinalGrid << std::endl; std::cout << "PromTuned: " << args.PromTuned << std::endl; diff --git a/Cell2Fire/ReadArgs.h b/Cell2Fire/ReadArgs.h index 537f654..42f6cf7 100644 --- a/Cell2Fire/ReadArgs.h +++ b/Cell2Fire/ReadArgs.h @@ -19,7 +19,7 @@ typedef struct NoOutput, verbose, IgnitionsLog, Ignitions, OutputGrids, FinalGrid, PromTuned, Stats, BBOTuning, AllowCROS; float ROSCV, ROSThreshold, CROSThreshold, HFIThreshold, HFactor, FFactor, BFactor, EFactor, FirePeriodLen; float CBDFactor, CCFFactor, ROS10Factor, CROSActThreshold; - int MinutesPerWP, MaxFirePeriods, TotalYears, TotalSims, NWeatherFiles, IgnitionRadius, seed, nthreads, FMC; + int MinutesPerWP, MaxFirePeriods, TotalYears, TotalSims, NWeatherFiles, IgnitionRadius, SpreadRadius, seed, nthreads, FMC; std::unordered_set HCells, BCells; } arguments; diff --git a/README.md b/README.md index 4552018..5574f26 100644 --- a/README.md +++ b/README.md @@ -39,6 +39,7 @@ There are many command line options available to configure the simulation. | nthreads | int | Defines the number of threads to use for simulations. Default is 1. | | seed | int | Seed for the random generator. | | IgnitionRad | int | If greater than 0, picks a point on the circumference of radius = IgnitionRadius from ignition point. | +| SpreadRad | int | Fire spread neighborhood radius. `1` uses immediate neighbors; larger values include additional tiers. | | nweathers | int | Number of weather files. Used to randomly select one. | | CCFFactor | float | Weighs by the cell’s CCF when calculating crown ROS in S&B. Default is 0. | | EFactor | float | Weighs radial distance from the ellipse center. Default is 1.0. | @@ -64,6 +65,15 @@ There are many command line options available to configure the simulation. | out-fl | bool | Saves flame length results. In S&B, saves surface, crown, and max. | | output-messages | bool | Saves file with messages (origin, destination, time, ROS). | +### Spread Neighborhood Radius (`--SpreadRad`) +`--SpreadRad` controls how far a burning cell can send spread messages in one neighborhood layer definition: + +- `--SpreadRad 1`: immediate neighborhood (legacy behavior) +- `--SpreadRad 2`: includes one additional ring of neighbors +- in general, interior cells have `(2n+1)^2 - 1` candidates for radius `n` + +The spread model keeps the original assumption that outgoing spread from source cell `i` uses the ROS computed at source cell `i` for all candidate destination cells. + # Output examples Most outputs are ASCII raster layers, of the same shape as the original fuels raster. diff --git a/test/spreadrad_test.sh b/test/spreadrad_test.sh new file mode 100644 index 0000000..a385d59 --- /dev/null +++ b/test/spreadrad_test.sh @@ -0,0 +1,69 @@ +#!/bin/bash +set -euo pipefail + +# Run a focused regression check for spread neighborhood radius. +# $1 is the binary suffix (same convention as test/test.sh). + +BIN="../Cell2Fire/Cell2Fire$1" +if [ ! -x "$BIN" ]; then + echo "Binary not found or not executable: $BIN" + exit 1 +fi + +WORKDIR="spreadrad_test_results" +rm -rf "$WORKDIR" +mkdir -p "$WORKDIR" + +run_case () { + local label="$1" + local spread_arg="$2" + local out="$WORKDIR/$label" + mkdir -p "$out" + # Keep this short but deterministic. + "$BIN" \ + --input-instance-folder model/kitral-asc \ + --output-folder "$out" \ + --nsims 30 \ + --max-fire-periods 60 \ + --sim K \ + --seed 123 \ + --ignitionsLog \ + $spread_arg > "$WORKDIR/$label.log" +} + +extract_burned () { + local logfile="$1" + grep "Total Burnt Cells:" "$logfile" | tail -n1 | sed -E 's/.*Total Burnt Cells:[[:space:]]*([0-9]+).*/\1/' +} + +run_case "default" "" +run_case "spread1" "--SpreadRad 1" +run_case "spread2" "--SpreadRad 2" + +B_DEFAULT="$(extract_burned "$WORKDIR/default.log")" +B_SPREAD1="$(extract_burned "$WORKDIR/spread1.log")" +B_SPREAD2="$(extract_burned "$WORKDIR/spread2.log")" + +if [ -z "$B_DEFAULT" ] || [ -z "$B_SPREAD1" ] || [ -z "$B_SPREAD2" ]; then + echo "Failed to parse burned-cell counts from logs" + exit 1 +fi + +if [ "$B_DEFAULT" != "$B_SPREAD1" ]; then + echo "SpreadRad regression failed: default burned=$B_DEFAULT, spread1 burned=$B_SPREAD1" + exit 1 +fi + +if ! grep -q "SpreadRadius: 2" "$WORKDIR/spread2.log"; then + echo "SpreadRad=2 run did not report SpreadRadius in output" + exit 1 +fi + +# Usually radius 2 burns at least as many cells as radius 1 for this fixed case. +if [ "$B_SPREAD2" -lt "$B_SPREAD1" ]; then + echo "Unexpected spread result: spread2 burned=$B_SPREAD2 < spread1 burned=$B_SPREAD1" + exit 1 +fi + +echo "SpreadRad test passed: default=$B_DEFAULT spread1=$B_SPREAD1 spread2=$B_SPREAD2" +rm -rf "$WORKDIR"