Skip to content

(rejection!=weighted sampling) make the probability map ignition generation algorithm more robust #183

Description

@fdobad

The ignition point generated when a probability map is provided is not what is expected.

Suppose you have three cells with probabilities [0.1, 0.5, 0.9].

Weighted sampling:
Cell 1: 0.1 / (0.1+0.5+0.9) ≈ 0.07
Cell 2: 0.5 / (0.1+0.5+0.9) ≈ 0.36
Cell 3: 0.9 / (0.1+0.5+0.9) ≈ 0.64
0.07 + 0.36 + 0.64 = 1.07
0.07 / 1.07 ≈ 0.065 ...etc.

Rejection sampling:   
Uniformly pick a cell (1/3 chance each).
Accept with probability equal to its ignition probability.
So, overall probability for cell 1: (1/3) * 0.1 ≈ 0.033
Cell 2: (1/3) * 0.5 ≈ 0.167
Cell 3: (1/3) * 0.9 ≈ 0.3
0.033 + 0.167 + 0.3 = 0.5 ... etc.

  Original	Normalized Weighted(1st)	Normalized Rejection(2nd)
  0.07	0.065	0.066
  0.36	0.336	0.334
  0.64	0.598	0.6
These are similar... but the implementation could be a lot better

On Cell2Fire.cpp :: Cell2Fire::RunIgnition when args.Ignitions == 0, the following occurs: Two nested while True loops attempt to chose a cell uniformly at random, then if another (0,1)random is bigger than the cell probabilityMap number, the cell is selected given availability and burned conditions; if this two conditions aren't met for both TotalCells*100 loop cycles, the ignition is marked as non realized.

A better algorithm: Weighted Sampling

  1. calculate the true available cells (in a array), and sum those probabilities
  2. weight sample on that bound: a uniform random number from 0 to that sum,
  3. cumulatively calculate which cell was selected
if (this->args.Ignitions == 0)
{

// 1. Build a vector of ignitable cell indices and their probabilities
std::vector<int> ignitable_cells;
std::vector<float> ignitable_probs;
float total_prob = 0.0f;

for (int i = 0; i < this->nCells; ++i) {
    // Filter: available, not burnt, not harvested, burnable type
    if (this->statusCells[i] < 3 &&
        this->burntCells.find(i + 1) == this->burntCells.end() &&
        this->fTypeCells[i] != 0) {
        ignitable_cells.push_back(i + 1);
        ignitable_probs.push_back(this->ignProb[i]);
        total_prob += this->ignProb[i];
    }
}

// 2. If no ignitable cells, set noIgnition and exit
if (ignitable_cells.empty()) {
    this->noIgnition = true;
    return true;
}

// 3. Weighted random sampling
boost::random::uniform_real_distribution<float> float_dist(0.0f, total_prob);
float rd_number = float_dist(generator2);

// Find which cell is selected
float cumulative = 0.0f;
int selected_cell = -1;
for (size_t idx = 0; idx < ignitable_cells.size(); ++idx) {
    cumulative += ignitable_probs[idx];
    if (rd_number <= cumulative) {
        selected_cell = ignitable_cells[idx];
        break;
    }
}

// 4. If no cell selected (should not happen), set noIgnition and exit
if (selected_cell == -1) {
    this->noIgnition = true;
    return true;
}

// 5. Initialize and ignite the selected cell
if (this->Cells_Obj.find(selected_cell) == this->Cells_Obj.end()) {
    InitCell(selected_cell);
}
it = this->Cells_Obj.find(selected_cell);

if (it->second.getStatus() == "Available" && it->second.fType != 0) {
    IgnitionHistory[sim] = selected_cell;
    std::vector<int> ignPts = { selected_cell };
    if (it->second.ignition(this->fire_period[year - 1],
                            this->year,
                            ignPts,
                            &df[selected_cell - 1],
                            this->coef_ptr,
                            this->args_ptr,
                            &(this->wdf[this->weatherPeriod]),
                            this->activeCrown,
                            this->perimeterCells)) {
        // Printing info about ignitions
        if (this->args.verbose) {
            std::cout << "Cell " << it->second.realId << " Ignites" << std::endl;
            std::cout << "Cell " << it->second.realId << " Status: " << it->second.getStatus() << std::endl;
        }
        this->statusCells[it->second.realId - 1] = 1;
        if (this->args.OutputGrids) {
            this->outputGrid();
        }
        this->noIgnition = false;
        return false;
    }
}

// If ignition fails, set noIgnition
this->noIgnition = true;
return true;

}

Metadata

Metadata

Labels

enhancementNew feature or request

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions