Add configurable spread neighborhood radius (SpreadRad)

.github/workflows/build-debian-stable.yml

@@ -115,4 +115,9 @@ jobs:
       - name: Bash Test
         run: |
           cd test
-          bash test.sh ${{ env.SUFFIX }}
+          bash test.sh ${{ env.SUFFIX }}
+
+      - name: SpreadRad Test
+        run: |
+          cd test
+          bash spreadrad_test.sh ${{ env.SUFFIX }}

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<int, std::vector<int>> 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);
                         }
                     }
                 }

Cell2Fire/Cells.cpp

@@ -96,7 +96,6 @@ Cells::Cells(int _id,
     this->angleDict = std::unordered_map<int, double>();
     this->ROSAngleDir = std::unordered_map<int, double>();
     this->distToCenter = std::unordered_map<int, double>();
-    this->angleToNb = std::unordered_map<int, int>();
 }
 
 /**
@@ -120,54 +119,34 @@ void
 Cells::initializeFireFields(std::vector<std::vector<int>>& coordCells,  // TODO: should probably make a coordinate type
                             std::unordered_set<int>& availSet,
                             int cols,
-                            int rows)  // WORKING CHECK OK
+                            int rows,
+                            int spreadRadius)  // WORKING CHECK OK
 {
-    std::vector<int> adj = adjacentCells(this->realId, rows, cols);
+    this->angleDict.clear();
+    this->ROSAngleDir.clear();
+    this->distToCenter.clear();
+    this->fireProgress.clear();
+
+    std::vector<int> 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<std::vector<int>>& coordCells,  // TODO:
  * -1.
  */
 std::vector<int>
-adjacentCells(int cell, int nrows, int ncols)
+adjacentCells(int cell, int nrows, int ncols, int radius)
 {
     if (cell <= 0 || cell > nrows * ncols)
     {
-        std::vector<int> adjacents(8, -1);
+        std::vector<int> 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<int> 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<int> 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<double> toPop = vector<double>();
 
         // 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<double> toPop = vector<double>();
 
         // 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)
     {

Cell2Fire/Cells.h

@@ -13,7 +13,7 @@
 
 using namespace std;
 
-std::vector<int> adjacentCells(int cell, int nrows, int ncols);
+std::vector<int> adjacentCells(int cell, int nrows, int ncols, int radius = 1);
 /*
  *   Weather structure
  */
@@ -96,11 +96,8 @@ class Cells
     std::unordered_map<int, std::vector<int>> gMsgListSeason;
     std::unordered_map<int, double> fireProgress;  // CP: dictionary {int: double}
     std::unordered_map<int, double> angleDict;     // CP: dictionary {int: double}
-    std::unordered_map<int, double> 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<int, double> ROSAngleDir;   // dictionary {neighbor id: ROS from this source cell}
     std::unordered_map<int, double> distToCenter;  // CP: dictionary {int: double}
-    std::unordered_map<int, int> angleToNb;        // CP: dictionary {double: int}
 
     // TODO: reference to shared object
 
@@ -117,7 +114,8 @@ class Cells
     void initializeFireFields(std::vector<std::vector<int>>& coordCells,
                               std::unordered_set<int>& 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

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;