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;

}

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

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