From 350ccab6f4d30186078bbf9e016b161c37474d13 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 26 Jun 2026 09:32:36 -0700 Subject: [PATCH] feat: add configurable spread neighborhood radius (SpreadRad) Port of the SpreadRad feature from cell2fire/Cell2Fire#136 and fire2a/C2FK#3, adapted to C2F-W. A new --SpreadRad option controls how far a burning cell can send spread messages: 1 (default) is the legacy immediate neighborhood; n includes (2n+1)^2 - 1 candidate neighbors. adjacentCells() takes a radius and returns the box of in-bounds neighbors. Because a radius greater than 1 makes several neighbors share an integer bearing, the per-cell spread bookkeeping (ROSAngleDir, angleDict, fireProgress, distToCenter) is rekeyed from angle to neighbor id, the angleToNb map is removed, and a new selectHeadCell() picks the neighbor closest to the wind azimuth for the slope calculation. The cleaning step in Cell2Fire.cpp erases the burnt cell from each neighbor's maps by id. SpreadRad 1 reproduces legacy results exactly; SpreadRad >= 2 spreads further. Adds test/spreadrad_test.sh and a CI step in build-debian.yml. Co-Authored-By: Claude Opus 4.8 (1M context) --- .github/workflows/build-debian.yml | 7 ++ Cell2Fire/Cell2Fire.cpp | 18 +-- Cell2Fire/Cells.cpp | 169 ++++++++++++++++++----------- Cell2Fire/Cells.h | 11 +- Cell2Fire/ReadArgs.cpp | 17 +++ Cell2Fire/ReadArgs.h | 2 +- README.md | 9 ++ test/spreadrad_test.sh | 83 ++++++++++++++ 8 files changed, 234 insertions(+), 82 deletions(-) create mode 100755 test/spreadrad_test.sh diff --git a/.github/workflows/build-debian.yml b/.github/workflows/build-debian.yml index faa0fea9..ae30ff49 100644 --- a/.github/workflows/build-debian.yml +++ b/.github/workflows/build-debian.yml @@ -156,3 +156,10 @@ jobs: chmod +x "test/test.sh" cd test ./test.sh "${{ env.SUFFIX }}" + + - name: Run SpreadRad test + run: |- + chmod +x "Cell2Fire/Cell2Fire${{ env.SUFFIX }}" + chmod +x "test/spreadrad_test.sh" + cd test + ./spreadrad_test.sh "${{ env.SUFFIX }}" diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index 68a6b8a2..a910d04c 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -603,7 +603,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) @@ -1520,18 +1520,18 @@ Cell2Fire::GetMessages(const std::unordered_map>& sendMess burntList.insert(it->second.realId); // Cleaning step - // Fire can't be propagated back - int cellNum = it->second.realId - 1; - for (auto& angle : it->second.angleToNb) + // Fire can't be propagated back: remove the burnt cell from + // every neighbor's spread bookkeeping (all keyed by cell id). + 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 45213b97..b4b33de0 100644 --- a/Cell2Fire/Cells.cpp +++ b/Cell2Fire/Cells.cpp @@ -8,6 +8,7 @@ #include "ReadCSV.h" #include "Spotting.h" // Include libraries +#include #include #include #include @@ -100,7 +101,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(); } /** @@ -125,9 +125,16 @@ 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) { @@ -138,41 +145,20 @@ Cells::initializeFireFields(std::vector>& coordCells, 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) - { - if (b >= 0) - angle = 270; - else - angle = 90; - } - else if (b == 0) - { - if (a >= 0) - angle = 180; - else - angle = 0; - } - else + // Angle (degrees, w.r.t. E-W, East positive, non-negative). atan2 + // keeps neighbors at the same compass bearing distinct across rings, + // which the old integer-bucketed scheme could not. + double angle = std::atan2(-1.0 * b, -1.0 * a) * radToDeg; + if (angle < 0.0) { - // 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; + // ROSAngleDir is keyed by neighbor id; -1 means "not yet set" + this->ROSAngleDir[nb] = -1; } - this->angleToNb[angle] = nb; this->fireProgress[nb] = 0.0; this->distToCenter[nb] = std::sqrt(a * a + b * b) * this->_ctr2ctrdist; } @@ -197,26 +183,82 @@ Cells::initializeFireFields(std::vector>& coordCells, * -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; } +/** + * @brief Selects the neighbor whose bearing is closest to the wind azimuth. + * + * Replaces the old angleToNb[waz] lookup, which assumed a unique neighbor at + * each integer angle. With a spread radius greater than 1, several neighbors + * can share a bearing, so the closest-by-angle (then closest-by-distance) + * neighbor is chosen. Returns this cell's own id when there are no neighbors. + */ +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; +} + /* New functions for calculating the ROS based on the fire angles Distribute the rate of spread (ROS,ros) to the axes given in the @@ -241,9 +283,10 @@ void Cells::ros_distr_old(double thetafire, double forward, double flank, double back) { // WORKING CHECK OK - for (auto& angle : this->ROSAngleDir) + // ROSAngleDir is keyed by neighbor id; look up the bearing via angleDict. + for (auto& nbRos : this->ROSAngleDir) { - double offset = std::abs(angle.first - thetafire); + double offset = std::abs(this->angleDict[nbRos.first] - thetafire); double base = ((int)(offset)) / 90 * 90; double result; @@ -265,7 +308,7 @@ Cells::ros_distr_old(double thetafire, double forward, double flank, double back { result = this->allocate(offset, 270, flank, forward); } - this->ROSAngleDir[angle.first] = result; + this->ROSAngleDir[nbRos.first] = result; } } @@ -321,10 +364,11 @@ void Cells::ros_distr_V2(double thetafire, double a, double b, double c, double EFactor) { - // Ros allocation for each angle inside the dictionary - for (auto& angle : this->ROSAngleDir) + // Ros allocation for each neighbor inside the dictionary + // ROSAngleDir is keyed by neighbor id; look up the bearing via angleDict. + for (auto& nbRos : this->ROSAngleDir) { - double offset = angle.first - thetafire; + double offset = this->angleDict[nbRos.first] - thetafire; if (offset < 0) { @@ -334,7 +378,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; } } @@ -448,7 +492,7 @@ Cells::manageFire(int period, head_angle = std::round(head_angle / 45.0) * 45.0; - int head_cell = angleToNb[head_angle]; // head cell for slope calculation + int head_cell = selectHeadCell(head_angle); // 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 @@ -623,8 +667,8 @@ Cells::manageFire(int period, // 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]; float ros = (1 + args->ROSCV * ROSRV) * _angle.second; float roundedRos = static_cast(std::ceil(ros * 100.0) / 100.0); @@ -805,7 +849,7 @@ Cells::manageFireBBO(int period, fire_struc headstruct, backstruct, flankstruct, metrics2; // Populate inputs - 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 @@ -976,8 +1020,8 @@ Cells::manageFireBBO(int period, // 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) @@ -1128,7 +1172,7 @@ Cells::get_burned(int period, fire_struc headstruct, backstruct, flankstruct; // Compute main angle and ROSs: forward, flanks and back - 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 @@ -1298,7 +1342,7 @@ Cells::ignition(int period, // << " bui: " << wdf_ptr->bui << std::endl; // Populate inputs - 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 @@ -1436,13 +1480,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 2ecb4723..c88676e9 100644 --- a/Cell2Fire/Cells.h +++ b/Cell2Fire/Cells.h @@ -12,7 +12,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 */ @@ -95,11 +95,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 @@ -116,7 +113,8 @@ class Cells void initializeFireFields(std::vector>& coordCells, std::unordered_set& availSet, int cols, - int rows); // TODO: need TYPE + int rows, + int spreadRadius); // TODO: need TYPE void ros_distr_old(double thetafire, double forward, double flank, double back); static 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: static double allocate(double offset, double base, double ros1, double ros2); static 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 184da890..707bb71d 100644 --- a/Cell2Fire/ReadArgs.cpp +++ b/Cell2Fire/ReadArgs.cpp @@ -249,6 +249,7 @@ parseArgs(int argc, char* argv[], arguments* args_ptr) int dmax_fire_periods = -1; int dseed = 123; int diradius = 0; + int dspreadradius = 1; int dnthreads = 1; int dfmc = 100; int dscen = 3; @@ -372,6 +373,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) @@ -640,6 +655,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; @@ -672,6 +688,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 00a824aa..09043dd4 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, UseWeatherWeights; 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, scenario; std::unordered_set HCells, BCells; } arguments; diff --git a/README.md b/README.md index c2ef54bd..64f5ff7b 100644 --- a/README.md +++ b/README.md @@ -105,6 +105,15 @@ Cell2Fire --final-grid --output-messages --out-ros --sim S --nsims 2 --seed 123 # check the results: to convert to tiff or see the results in QGIS, use the plugin ``` +#### Spread neighborhood radius (`--SpreadRad`) +`--SpreadRad` controls how far a burning cell can send spread messages in one neighborhood-layer definition: + +- `--SpreadRad 1` (default): immediate neighborhood (legacy behavior) +- `--SpreadRad 2`: includes one additional ring of neighbors +- in general, interior cells have `(2n+1)^2 - 1` candidate neighbors for radius `n` + +The spread model keeps the original assumption that outgoing spread from a source cell `i` uses the ROS computed at source cell `i` for all candidate destination cells. + ### Containerized TL;DR: ```bash diff --git a/test/spreadrad_test.sh b/test/spreadrad_test.sh new file mode 100755 index 00000000..a5edbf5d --- /dev/null +++ b/test/spreadrad_test.sh @@ -0,0 +1,83 @@ +#!/bin/bash +set -euo pipefail + +# Focused regression check for the spread neighborhood radius (--SpreadRad). +# $1 is the binary suffix (same convention as test/test.sh). +# +# Invariants: +# * default run == "--SpreadRad 1" (radius 1 is the legacy default) +# * "--SpreadRad 2" is reported and burns at least as many cells as radius 1 + +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" +} + +# Per-simulation burned-cell counts, as a comma-separated sequence. +# C2F-W prints a per-sim summary table with a "Burnt N X%" row. +burned_sequence () { + local logfile="$1" + grep -E '^[[:space:]]+Burnt[[:space:]]' "$logfile" | awk '{print $2}' | paste -sd, +} + +sum_burned () { + local logfile="$1" + grep -E '^[[:space:]]+Burnt[[:space:]]' "$logfile" | awk '{s += $2} END {print s + 0}' +} + +run_case "default" "" +run_case "spread1" "--SpreadRad 1" +run_case "spread2" "--SpreadRad 2" + +SEQ_DEFAULT="$(burned_sequence "$WORKDIR/default.log")" +SEQ_SPREAD1="$(burned_sequence "$WORKDIR/spread1.log")" +S_SPREAD1="$(sum_burned "$WORKDIR/spread1.log")" +S_SPREAD2="$(sum_burned "$WORKDIR/spread2.log")" + +if [ -z "$SEQ_DEFAULT" ] || [ -z "$SEQ_SPREAD1" ]; then + echo "Failed to parse burned-cell counts from logs" + exit 1 +fi + +if [ "$SEQ_DEFAULT" != "$SEQ_SPREAD1" ]; then + echo "SpreadRad regression failed: default and --SpreadRad 1 differ" + echo " default : $SEQ_DEFAULT" + echo " spread1 : $SEQ_SPREAD1" + exit 1 +fi + +if ! grep -q "SpreadRadius: 2" "$WORKDIR/spread2.log"; then + echo "SpreadRad=2 run did not report the spread radius in output" + exit 1 +fi + +# Usually radius 2 burns at least as many cells as radius 1 for this fixed case. +if [ "$S_SPREAD2" -lt "$S_SPREAD1" ]; then + echo "Unexpected spread result: spread2 total burned=$S_SPREAD2 < spread1 total burned=$S_SPREAD1" + exit 1 +fi + +echo "SpreadRad test passed: default==spread1 (total=$S_SPREAD1), spread2 total=$S_SPREAD2" +rm -rf "$WORKDIR"