123 change statistics format

Cell2Fire/Cell2Fire.cpp

@@ -49,6 +49,11 @@ std::unordered_map<int, std::vector<int>> HarvestedCells;
 std::vector<int> NFTypesCells;
 std::unordered_map<int, int> IgnitionHistory;
 std::unordered_map<int, std::string> WeatherHistory;
+std::unordered_map<int, std::vector<float>> StatisticsPerCell;
+std::unordered_map<int, float> meanSurfaceFlameLength;
+std::unordered_map<int, float> meanCrownFlameLength;
+std::unordered_map<int, float> meanMaxFlameLength;
+std::unordered_map<int, std::vector<float>> StatisticsPerSim;
 
 /******************************************************************************
                                                                                                                                 Utils
@@ -176,9 +181,6 @@ Cell2Fire::Cell2Fire(arguments _args) : CSVForest(_args.InFolder + "fuels", " ")
     this->args = _args;
     this->args_ptr = &this->args;
 
-    // this->wdf_ptr = new weatherDF;
-    // this->wdf[100000];
-
     /********************************************************************
      *
      *															Initialization
@@ -341,7 +343,6 @@ Cell2Fire::Cell2Fire(arguments _args) : CSVForest(_args.InFolder + "fuels", " ")
             this->harvestCells.insert(i + 1);
     }
 
-    // POTENCIALMENTE AQUI ESTA MALO
     /* Weather DataFrame */
     // this->WeatherData = this->CSVWeather.getData();
     // std::cout << "\nWeather DataFrame from instance " << this->CSVWeather.fileName << std::endl;
@@ -673,7 +674,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1)
         this->surfaceIntensityFolder = Cell2Fire::createOutputFolder("SurfaceIntensity");
     }
     // Crown Byram Intensity Folder
-    if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->crownIntensityFolder = Cell2Fire::createOutputFolder("CrownIntensity");
     }
@@ -683,15 +684,19 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1)
         this->surfaceFlameLengthFolder = Cell2Fire::createOutputFolder("SurfaceFlameLength");
     }
     // Crown Flame Length Folder
-    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->crownFlameLengthFolder = Cell2Fire::createOutputFolder("CrownFlameLength");
     }
     // max Flame Length Folder
-    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->maxFlameLengthFolder = Cell2Fire::createOutputFolder("MaxFlameLength");
     }
+    if (this->args.Stats)
+    {
+        this->statsFolder = Cell2Fire::createOutputFolder("Statistics");
+    }
     // Crown Folder
     if (this->args.OutCrown && this->args.AllowCROS)
     {
@@ -1749,7 +1754,7 @@ Cell2Fire::Results()
     }
 
     // Crown Intensity
-    if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->crownIntensityFolder = this->args.OutFolder + "CrownIntensity" + separator();
         std::string intensityName;
@@ -1785,7 +1790,7 @@ Cell2Fire::Results()
     }
 
     // Crown Flame length
-    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator();
         std::string fileName;
@@ -1803,7 +1808,7 @@ Cell2Fire::Results()
     }
 
     // Max Flame length
-    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S"))
+    if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C"))
     {
         this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator();
         std::string fileName;
@@ -1820,6 +1825,58 @@ Cell2Fire::Results()
             this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths);
     }
 
+    // Flame Length StatisticsPerCell
+    // This will be used to store accumulated flame lengths across simulations
+    if (this->args.Stats)
+    {
+        float totalCrown = 0;
+        float totalSurface = 0;
+        float maxPerSim = 0;
+        std::vector<float> simStats;
+        for (int cell : this->burntCells)
+        {
+            std::vector<float> cellFlameLengthMeans;
+
+            float surfaceFlameLength = this->surfaceFlameLengths[cell] / args.TotalSims;
+            meanSurfaceFlameLength[cell] += surfaceFlameLength;
+            cellFlameLengthMeans.push_back(meanSurfaceFlameLength[cell]);
+            totalSurface += this->surfaceFlameLengths[cell];
+            maxPerSim = max(maxPerSim, this->surfaceFlameLengths[cell]);
+
+            if ((this->args.AllowCROS) && (this->args.Simulator != "C"))
+            {
+                float crownFlameLength = this->crownFlameLengths[cell] / args.TotalSims;
+                float maxFlameLength = this->maxFlameLengths[cell] / args.TotalSims;
+                meanCrownFlameLength[cell] += crownFlameLength;
+                meanMaxFlameLength[cell] += maxFlameLength;
+                cellFlameLengthMeans.push_back(meanCrownFlameLength[cell]);
+                cellFlameLengthMeans.push_back(meanMaxFlameLength[cell]);
+                totalCrown += this->maxFlameLengths[cell];
+                maxPerSim = max(maxPerSim, this->crownFlameLengths[cell]);
+            }
+            StatisticsPerCell[cell] = { cellFlameLengthMeans };
+        }
+        simStats.push_back(totalSurface / BCells);
+        if ((this->args.AllowCROS) && (this->args.Simulator != "C"))
+        {
+            simStats.push_back(totalCrown / BCells);
+        }
+        simStats.push_back(maxPerSim);
+
+        StatisticsPerSim[this->sim] = { simStats };
+        std::ostringstream oss;
+        std::string Stats = this->statsFolder + "statisticsPerSim" + oss.str() + ".csv";
+        CSVWriter simStatsFile(Stats);
+        simStatsFile.printStats(StatisticsPerSim, "sim", (this->args.AllowCROS) && (this->args.Simulator != "C"));
+        if (currentSim == args.TotalSims)
+        {
+            std::ostringstream oss;
+            std::string Statsname = this->statsFolder + "statisticsPerCell" + oss.str() + ".csv";
+            CSVWriter statsFile(Statsname);
+            statsFile.printStats(StatisticsPerCell, "cell", (this->args.AllowCROS) && (this->args.Simulator != "C"));
+        }
+    }
+
     // Intensity
     if ((this->args.OutCrownConsumption) && (this->args.AllowCROS))
     {

Cell2Fire/Cell2Fire.h

@@ -89,6 +89,7 @@ class Cell2Fire
     string crownFlameLengthFolder;
     string historyFolder;
     string ignitionsFolder;
+    string statsFolder;
 
     // Vectors
     std::vector<int> fire_period;

Cell2Fire/Cells.cpp

@@ -677,7 +677,7 @@ Cells::manageFire(int period,
                 surfFraction[nb] = metrics.sfc;
                 SurfaceFlameLengths[this->realId - 1] = mainstruct.fl;
                 SurfaceFlameLengths[nb - 1] = metrics.fl;
-                if ((args->AllowCROS) && (args->Simulator == "S"))
+                if ((args->AllowCROS) && (args->Simulator != "C"))
                 {
                     float comp_zero = 0;
                     MaxFlameLengths[this->realId - 1]

Cell2Fire/FuelModelKitral.cpp

@@ -705,6 +705,51 @@ byram_intensity(inputs* data, main_outs* at)
     return ib;  // unidad de medida
 }
 
+/**
+ * Calculates byram fire intensity when there is active crown fire.
+ * In order for this to be calculated, the input folder must contain
+ * files with CBD, CBH and tree height data for each cell.
+ * @param at Structure containing the cell's output data.
+ * @param data Structure containing the cell's input data.
+ * @return Fire intensity.
+ */
+
+float
+crown_byram_intensity_k(main_outs* at, inputs* data)
+{
+    float cbd, cbh;
+    cbd = data->cbd;
+    cbh = data->cbh;
+    // CBH is 0,1012 * Height^(1,4822)
+    float canopy_height = std::pow(cbh / 0.1012, 1 / 1.4822) - cbh;
+    if (canopy_height < 0)
+    {
+        throw std::runtime_error("Tree height is lower than canopy base height, "
+                                 "please provide valid files.");
+    }
+    return std::ceil((hs[data->nftype][0] / 60) * cbd * canopy_height * at->ros_active * 100.0) / 100.0;
+}
+
+/**
+ * @brief Calculates the flame length of a cell when there is crown fire.
+ * @param intensity Byram intensity for crown fires
+ * @return the flame length
+ */
+
+float
+crown_flame_length_k(float intensity)
+{
+    float fl = 0.1 * pow(intensity, 0.5);
+    if (fl < 0.01)
+    {
+        return 0;
+    }
+    else
+    {
+        return std::ceil(fl * 100.0) / 100.0;
+    }
+}
+
 bool
 fire_type(inputs* data, main_outs* at, int FMC)
 {
@@ -935,6 +980,8 @@ calculate_k(inputs* data,
         at->rss = hptr->ros;
         bptr->ros = backfire_ros10_k(hptr, sec);
         fptr->ros = flankfire_ros_k(hptr->ros, bptr->ros, sec->lb);
+        at->crown_intensity = crown_byram_intensity_k(at, data);
+        at->crown_flame_length = crown_flame_length_k(at->crown_intensity);
 
         if (args->verbose)
         {
@@ -957,6 +1004,8 @@ calculate_k(inputs* data,
         at->rss = hptr->ros;
         bptr->ros = backfire_ros10_k(hptr, sec);
         fptr->ros = flankfire_ros_k(hptr->ros, bptr->ros, sec->lb);
+        at->crown_intensity = crown_byram_intensity_k(at, data);
+        at->crown_flame_length = crown_flame_length_k(at->crown_intensity);
 
         if (args->verbose)
         {
@@ -1045,6 +1094,8 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main
         if (crownFire)
         {
             metrics->cfb = crownfractionburn(data, metrics, FMC);
+            metrics->crown_intensity = crown_byram_intensity_k(metrics, data);
+            metrics->crown_flame_length = crown_flame_length_k(metrics->crown_intensity);
         }
         if (args->verbose)
         {

Cell2Fire/FuelModelKitral.h

@@ -32,6 +32,9 @@ float backfire_ros_k(main_outs* at, snd_outs* sec);
 // Flame length [m])
 float flame_length(inputs* data, main_outs* at);
 
+// Crown Flame length [m])
+float crown_flame_length_k(float intensity);
+
 // Angle of the flame w.r.t. horizontal surface (Putnam's)
 float angleFL(inputs* data, main_outs* at);
 
@@ -41,6 +44,9 @@ float flame_height(inputs* data, main_outs* at);
 // byram intensity
 float byram_intensity(inputs* data, main_outs* at);
 
+// crown byram intensity
+float crown_byram_intensity_k(main_outs* at, inputs* data);
+
 // Type of fire (Crown = 1, SROS = 0)
 bool fire_type(inputs* data, main_outs* atr, int FMC);
 // Active crown

Cell2Fire/FuelModelSpain.h

@@ -35,6 +35,8 @@ float backfire_ros_s(main_outs* at, snd_outs* sec);
 
 // Flame length [m])
 float flame_length(inputs* data, fuel_coefs* ptr);
+// Crown Flame length [m])
+float crown_flame_length(float intensity);
 
 // Angle of the flame w.r.t. horizontal surface (Putnam's)
 float angleFL(inputs* data, fuel_coefs* ptr);
@@ -45,6 +47,9 @@ float flame_height(inputs* data, fuel_coefs* ptr);
 // byram intensity
 float byram_intensity(inputs* data, fuel_coefs* ptr);
 
+// crown byram intensity
+float crown_byram_intensity(main_outs* at, inputs* data);
+
 // Type of fire (Crown = 1, SROS = 0)
 bool fire_type(inputs* data, fuel_coefs* ptr);
 

Cell2Fire/ReadArgs.cpp

@@ -206,16 +206,17 @@ parseArgs(int argc, char* argv[], arguments* args_ptr)
     //}
 
     //--statistics
-    // if (cmdOptionExists(argv, argv + argc, "--statistics")) {
-    //	out_stats = true;
-    //	printf("Statistics: %d \n", out_stats);
-    //}
+    if (cmdOptionExists(argv, argv + argc, "--statistics"))
+    {
+        out_stats = true;
+        printf("StatisticsPerCell: %d \n", out_stats);
+    }
 
     //--bbo
     if (cmdOptionExists(argv, argv + argc, "--bbo"))
     {
         bbo_tuning = true;
-        printf("BBOTuning: %s \n", btoa(out_stats));
+        printf("BBOTuning: %s \n", btoa(bbo_tuning));
     }
 
     //--cros
@@ -628,7 +629,7 @@ printArgs(arguments args)
     std::cout << "FinalGrid: " << args.FinalGrid << std::endl;
     std::cout << "PromTuned: " << args.PromTuned << std::endl;
     std::cout << "BBOTuning: " << args.BBOTuning << std::endl;
-    std::cout << "Statistics: " << args.Stats << std::endl;
+    std::cout << "StatisticsPerCell: " << args.Stats << std::endl;
     std::cout << "noOutput: " << args.NoOutput << std::endl;
     std::cout << "verbose: " << args.verbose << std::endl;
     std::cout << "seed: " << args.seed << std::endl;
@@ -660,7 +661,7 @@ printArgs(arguments args)
     std::cout << "FinalGrid: " << args.FinalGrid << std::endl;
     std::cout << "PromTuned: " << args.PromTuned << std::endl;
     std::cout << "BBOTuning: " << args.BBOTuning << std::endl;
-    std::cout << "Statistics: " << args.Stats << std::endl;
+    std::cout << "StatisticsPerCell: " << args.Stats << std::endl;
     std::cout << "noOutput: " << args.NoOutput << std::endl;
     std::cout << "verbose: " << args.verbose << std::endl;
     std::cout << "Ignition Points Log: " << args.IgnitionsLog << std::endl;

Cell2Fire/WriteCSV.cpp

@@ -287,6 +287,48 @@ CSVWriter::printIgnitions(std::unordered_map<int, int> ignitionsHistory,
     ofs.close();
 }
 
+void
+CSVWriter::printStats(std::unordered_map<int, std::vector<float>> statistics, std::string type, bool hasCrown)
+{
+    std::ofstream ofs(fileName, std::ofstream::out);
+
+    // Print column titles
+    ofs << type << ",surfaceFlameLengthMean";
+    if (hasCrown)
+    {
+        ofs << type << ",crownFlameLengthMean,maxFlameLength";
+    }
+    else if (type == "sim")
+    {
+        ofs << ",maxFlameLength";
+    }
+    ofs << std::endl;
+
+    // Iterate through the unordered_map using iterator
+    for (auto it = statistics.begin(); it != statistics.end(); ++it)
+    {
+        int sim = it->first;                                  // The simulation ID (key)
+        const std::vector<float>& flameLengths = it->second;  // The associated flame lengths (value)
+
+        ofs << sim;  // Print the simulation number (key)
+
+        // Print all the flame lengths associated with the current simulation
+        for (size_t j = 0; j < flameLengths.size(); j++)
+        {
+            ofs << "," << flameLengths[j];
+            if (j < flameLengths.size() - 1)
+            {
+                ofs << " ";  // Optional: space between multiple flame lengths
+            }
+        }
+
+        ofs << "\n";  // Move to the next line for each entry
+    }
+
+    // Close the file
+    ofs.close();
+}
+
 void
 CSVWriter::printCSV_V2(int rows, int cols, std::vector<int> statusCells)
 {

Cell2Fire/WriteCSV.h

@@ -45,6 +45,7 @@ class CSVWriter
     void printWeather(std::vector<std::string> weatherHistory);
     void printIgnitions(std::unordered_map<int, int> ignitionsHistory,
                         std::unordered_map<int, std::string> weatherHistory);
+    void printStats(std::unordered_map<int, std::vector<float>> statistics, std::string type, bool hasCrown);
     void printCSV_V2(int rows, int cols, std::vector<int> statusCells);
     // Function to create a directory
     void MakeDir(std::string pathPlot);