(novembre 2015) Simulation d'un feu de forêt

# /usr/bin/env python3
# coding: utf-8

"""Simulation de percolation dans le cadre des feux de forêt.
Cette version est a deux paramètres :
- la probabilité *p* que le feux se propage
- *l* la taille du système.
Les conditions aux bords sont fermées (_open boundary_) et aucun parallélisme ne sera tenté.
Entrée : aucune ;
Sortie : flottant entre 0 et 1, r, proportion d'arbre restant après l'incendie.
         et une image.
"""

def voisins(case):
    """Renvoie les voisins de la case via un générateur."""
    x=case[0]
    y=case[1]
    yield x+1,y
    yield x-1,y
    yield x,y+1
    yield x,y-1
    
from PIL import Image, ImageDraw
import numpy as np
from random import randint, random
p = 0.55
l = 250

monde = np.zeros((l, l), dtype=int) #Convention : 0 = arbre, 1 = feu, 2 = cendre
#Conditions aux limites: on entoure de cendre
monde[0:,0] = 2 ; monde[0:,l-1] = 2 ; monde[0,0:] = 2 ;monde[l-1,0:] = 2 

monde[randint(1,l-2), randint(1, l-2)] = 1 #1er feu

#On continue tant qu'il reste du feu (1) sur la grille.
while 1 in monde:
    #Boucle principale
    en_feu = []
    #Recherche des cases en feu, de manière pas du tout optimale.
    for i in range(1,l):
        for j in range(1,l):
            if monde[i,j] == 1:
                en_feu.append((i,j))
    #Enflammement des arbres.
    for case in en_feu:
        for voisin in voisins(case):
            #va tester chaque voisin tels que définis dans la fonction.
            if monde[voisin] == 0 and random()<p:
                monde[voisin] = 1
        monde[case] = 2

#Post-traitement
restant = 0
#Création d'une image
img = Image.new('RGBA', (l, l), (255,255,255,0))
d = ImageDraw.Draw(img)
for i in range(0,l):
    for j in range(0,l):
        if monde[i,j] == 0:
            #En passant, on regarde combien d'arbre il reste.
            restant += 1
            d.rectangle((i, j, i+1, j+1), (0, 150, 50, 255))#vert
        else:
            #Le cas "feu" n'est pas possible.
            #Si ce n'est pas un arbre, c'est donc des cendres.
            d.rectangle((i, j, i+1, j+1), (50, 50, 50, 255))#gris
print(restant/(l-2)**2)#Proportion d'arbre restant.
img.save("perco_arbre"+str(p)+".png", "PNG")

import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors

# Pour s'épargner des "RuntimeError"
sys.setrecursionlimit(3000)

# Paramètres
taille_zone = 50,50

# Constantes
CENDRE = 0
ARBRE = 1
EAU = 2

# ALLUUUMEEEEEEER LE FEU !
def allumer_le_feu(foret, position, probabilite_embrasement):
   foret[position] = CENDRE
   decalages = [(-1,-1), ( 0,-1), ( 1,-1), ( 1, 0),
               ( 1, 1), ( 0, 1), (-1, 1), (-1, 0)]
   coord_voisins = [(position[0]+x, position[1]+y) for x, y in decalages]
   for v in coord_voisins:
       if (v[0] > 0 and v[0] < taille_zone[0] -1
      and v[1] > 0 and v[1] < taille_zone[1] -1         
      and foret[v] == ARBRE         
      and np.random.rand(1) < probabilite_embrasement):
           allumer_le_feu(foret,v,probabilite_embrasement)

# Une forêt carrée, bordée de cendres, avec un lac au milieu.
def planter_foret():
   foret = np.zeros(taille_zone,dtype='int_')
   foret.fill(CENDRE)
   foret[1:-1,1:-1].fill(ARBRE)
   foret[20:30,20:-10].fill(EAU)
   return foret


# Calcul

nb_echantillons = 60
proba_embrasement = [i/80 for i in range(80)]
proba_embrasement.extend([i/100 for i in range(20,30)])
res_proba = []
res_ratio = []
for proba in proba_embrasement:
   for i in range(nb_echantillons):
       foret = planter_foret()
       nb_arbres_total = np.count_nonzero(foret == ARBRE)
       position_foyer = np.random.randint(0,taille_zone[0]), np.random.randint(0,taille_zone[1])
       allumer_le_feu(foret, position_foyer, proba)
       nb_arbres_restant = np.count_nonzero(foret == ARBRE)
       res_proba.append(proba)
       res_ratio.append(nb_arbres_restant/nb_arbres_total)

# Graphe
plt.plot(res_proba,res_ratio,".")
plt.show()


# Affichage de la forêt initiale

cmap = colors.ListedColormap(['black', 'green','blue'])
bounds=[0,1,2,3]
norm = colors.BoundaryNorm(bounds, cmap.N)
plt.imshow(planter_foret(),interpolation="nearest",cmap=cmap, norm=norm)
plt.show()

#include <array>
#include <chrono>
#include <iostream>
#include <random>
#include <thread>

using namespace std::chrono_literals;

void clear() {
  std::cout << "\x1B[2J\x1B[H";
}

enum Cell { Live, Burning, Ashes };

template<int N>
using State = std::array<std::array<Cell, N>, N>;

template<int N>
void print(const State<N>& s)
{
  clear();
  for(const auto &l : s) {
    for(Cell c : l) {
      std::cout <<
        (c == Live ? "T":
         c == Burning ? "B":
         ".");
    }
    std::cout << std::endl;
  }
}

template<int N>
bool update(State<N>& s, double p) {
  State<N> newS{s};
  std::random_device rd;
  std::mt19937 g(rd());
  std::uniform_real_distribution<double> dis(0, 1);
  bool updated = false;

  for(int i=0; i<N; ++i) {
    for(int j=0; j<N; ++j) {
      if(s[i][j] == Burning) {
        if(i > 0 && newS[i-1][j] == Live && dis(g) < p) {
          newS[i-1][j] = Burning;
          updated = true;
        }
        if(i < N-1 && newS[i+1][j] == Live && dis(g) < p) {
          newS[i+1][j] = Burning;
          updated = true;
        }
        if(j > 0 && newS[i][j-1] == Live && dis(g) < p) {
          newS[i][j-1] = Burning;
          updated = true;
        }
        if(j < N-1 && newS[i][j+1] == Live && dis(g) < p) {
          newS[i][j+1] = Burning;
          updated = true;
        }
        newS[i][j] = Ashes;
      }
    }
  }

  std::swap(s, newS);
  return updated;
}

int main() {
  const int size = 60;
  const double p = 0.55;
  const int max_iter = 100;
  int total_burned = 0;

  for(int iter=0; iter<max_iter; ++iter) {
    State<size> state{};
    state[size/2][size/2] = Burning;
    do {
      print<size>(state);
      std::this_thread::sleep_for(20ms);
    } while(update<size>(state, p));
    print<size>(state);
    int burned = 0;
    for(const auto &l : state) {
      for(Cell c : l) {
        if(c == Ashes) {
          burned++;
        }
      }
    }
    total_burned += burned;
  }
  std::cout << "Probability " << p << std::endl;
  std::cout << "Average burning ratio " << (total_burned*100)/(size*size*max_iter) << "%" << std::endl;

  return 0;
}

import random
from itertools import product
import pygame
from pygame.locals import *


# variables globales
C_SZ = 2                            # taille des cases
WIDTH, HEIGHT = 512, 256            # nombre de cases (hauteur / largeur)
COLORS = 'green', 'gray', 'red'     # couleurs possibles des cases (sauf fond)
STATES = 'arbre', 'cendre', 'brule' # etats des cases
PROBA = 0.60                        # probabilite qu'une case soit occupee par un arbre
DELAY = 10                          # delai entre 2 frames (ms)
SURFACES = {}                       # dictionnaire de surfaces (color => sprite)
S2C = dict(zip(STATES, COLORS))     # correspondance (state => color)


# initialise les surfaces à blitter
def init_surfaces():
    for color in COLORS:
        SURFACES[color] = pygame.Surface((C_SZ, C_SZ))
        SURFACES[color].fill(pygame.Color(color))

# affichage
def display(screen, dct):
    for (i,j), content in dct.items():
        screen.blit(SURFACES[S2C[content]], (j * C_SZ, i * C_SZ))
    pygame.display.flip()

# avance d'un tour
def step(mat, dct):
    new_dct = {}
    for (i,j), content in dct.items():
        if content == 'brule':
            new_dct[i,j] = 'cendre'
            for x, y in [(-1,0), (1,0), (0,-1), (0,1)]:
                pos = i + x, j + y
                if mat.get(pos) == 'arbre':
                    new_dct[pos] = 'brule'
    return new_dct


if __name__ == '__main__':
    # initialise aléatoirement la position des arbres sur le terrain
    mat = {pos: 'arbre'
           for pos in product(range(HEIGHT), range(WIDTH))
           if random.random() < PROBA}
    N = len(mat)

    pygame.init()
    init_surfaces()
    screen = pygame.display.set_mode((C_SZ * WIDTH, C_SZ * HEIGHT))
    screen.fill(pygame.Color('yellow'))

    display(screen, mat)

    dct = {random.choice(list(mat.keys())): 'brule'}    # premier arbre en feu
    while dct:
        pygame.time.wait(DELAY)
        display(screen, dct)
        mat.update(dct)
        dct = step(mat, dct)

    K = N - len([v for v in mat.values() if v == 'arbre'])
    print("proba = {}, {:.2%} d'arbres brules ({} / {})".format(PROBA, K / N, K, N))

    pygame.quit()

import random

import pygame
from pygame.locals import *


class Tree:
   BURNED = 0
   BURNING = 1
   ALIVE = 2
   def __init__(self):
       self.state = self.ALIVE
   def __repr__(self):
       s = {self.BURNED:'BURNED', self.BURNING:'BURNING', self.ALIVE:'ALIVE'}[self.state]
       return '<Tree ' + s + '>'

class World:
   def __init__(self, s, p):
       self.map = [list([Tree() for i in range(s)]) for j in range(s)]
       self.burning = []
       self.new_burning = []
       self.burned = []

       self.there_is_fire = False

       self.nb_alive = s ** 2
       self.nb_burned = 0

       self.surface = pygame.Surface((s,s))
       self.surface.fill((0,255,0))

       self.font = pygame.font.Font(None,24)
       self.txt = None

       self.p = p

   def fire_on(self, pos):
       x,y = pos
       if x<0 or x>=len(self.map) or y<0 or y>=len(self.map):
           return
       if self.map[x][y].state is not Tree.ALIVE:
           return
       self.map[x][y].state = Tree.BURNING
       if not self.there_is_fire:
           self.burning = [(x,y)]
           self.there_is_fire = True
           self.nb_burned = 0

       self.nb_burned += 1
       self.new_burning.append((x,y))

   def on_update(self):
       if not self.there_is_fire:
           self.txt = self.font.render("{} % d'arbres brûlés".format(round(self.nb_burned/self.nb_alive * 100, 2)), \
                                       True, (0,0,0), (50,200,255))
           return
       elif self.txt:
           self.txt = None
       self.there_is_fire = False

       self.new_burning = []

       #self.surface.fill((0,255,0))
       burner = lambda : random.random() < self.p

       for x,y in self.burned:
           self.surface.set_at((x,y),(0,0,0))
       self.burned = []

       for x,y in self.burning:
           self.there_is_fire = True
           self.surface.set_at((x,y), (255,0,0))

           self.map[x][y].state = Tree.BURNED

           up = burner()
           down = burner()
           left = burner()
           right = burner()
           if up:
               self.fire_on((x,y+1))
           if down:
               self.fire_on((x,y-1))
           if left:
               self.fire_on((x-1,y))
           if right:
               self.fire_on((x+1,y))

           self.burned.append((x,y))
       self.burning = self.new_burning[:]

   def on_render(self,dst):
       dst.blit(self.surface, (0,0))
       if self.txt:
           dst.blit(self.txt, (0,0))

class App:
   def __init__(self):
       self.window = pygame.display.set_mode((512, 512))
       self.running = False

       self.world = World(512, 0.5)

       self.font = pygame.font.Font(None, 24)
       self.menu = self.font.render("Appuyez sur 'R' pour ré-initialiser la simulation.", \
                                    True, (0,0,0),(50,200,255))

       self.clock= pygame.time.Clock()

   def on_event(self, e):
       if e.type == QUIT:
           self.running = False
       elif e.type == MOUSEBUTTONDOWN:
           self.world.fire_on(e.pos)
       elif e.type == KEYDOWN:
           if e.key == K_r:
               p = self.world.p
               self.world = World(512,p)
           elif e.key == K_LEFT:
               self.world.p -= 0.1
               if self.world.p < 0:
                   self.world.p = 0
           elif e.key == K_RIGHT:
               self.world.p += 0.1
               if self.world.p > 1:
                   self.world.p = 1

   def on_render(self):
       self.world.on_render(self.window)
       if not self.world.there_is_fire:
           prob_txt = "Probabilité : < {} >".format(str(round(self.world.p,2)))
           prob = self.font.render(prob_txt, True, (0,0,0), (50,200,255))
           x = (self.window.get_width() - self.menu.get_width()) // 2
           y = (self.window.get_height() - self.menu.get_height()) // 2
           self.window.blit(self.menu,(x,y))
           x = (self.window.get_width() - prob.get_width()) // 2
           y += 1.2 * self.menu.get_height()
           self.window.blit(prob, (x,y))

   def on_update(self):
       self.world.on_update()

   def on_mainloop(self):
       self.running = True

       while self.running:
           for e in pygame.event.get():
               self.on_event(e)
           self.on_update()
           self.on_render()
           pygame.display.flip()
           pygame.display.set_caption(str(round(self.clock.get_fps(),2))+" FPS")
           self.clock.tick(70)

if __name__ == "__main__":
   pygame.init()
   a = App()
   a.on_mainloop()
   pygame.quit()

! Gestion de l'alétoire.
 module alea
 use iso_fortran_env, only: int64
 implicit none
 contains 
 subroutine init_random_seed()
 integer, allocatable :: seed(:)
 integer :: i, n, dt(8), pid
 integer(int64) :: t
 
 call random_seed(size = n)
 allocate(seed(n))
 call system_clock(t)
 call date_and_time(values=dt)
 t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
 + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
 + dt(3) * 24_int64 * 60 * 60 * 1000 &
 + dt(5) * 60 * 60 * 1000 &
 + dt(6) * 60 * 1000 + dt(7) * 1000 &
 + dt(8)
 pid = getpid()
 t = ieor(t, int(pid, kind(t)))
 do i = 1, n
   seed(i) = lcg(t)
   end do
   call random_seed(put=seed)
   
   contains
   ! This simple PRNG might not be good enough for real work, but is
   ! sufficient for seeding a better PRNG.
   function lcg(s)
   integer :: lcg
   integer(int64) :: s
   s = mod(s, 4294967296_int64)
   s = mod(s * 279470273_int64, 4294967291_int64)
   lcg = int(mod(s, int(huge(0), int64)), kind(0))
   end function lcg
end subroutine init_random_seed
end module alea
   
   
program feu
! programme principale.
use alea
integer, dimension(:,:), allocatable :: monde
real, dimension(:,:), allocatable :: stocha
integer, dimension(:,:), allocatable :: en_feu, en_feu2
integer :: l, i, j, k, k2, nb_en_feu, ok
real :: pa, x, y

!Les nx/ny permettent de facilement passer d'un réseau à un autre.
!Il suffit de modifier une variable pour passer à un réseau 6.
integer, parameter, dimension(6) :: nx=(/1,-1, 0, 0, 1,-1/)
integer, parameter, dimension(6) :: ny=(/0, 0, 1,-1, 1,-1/)
integer, parameter :: res = 4

l = 1000
pa = 1!proba d'un arbre.

!Convention : 0 = arbre, 1 = lac, 2 = feu, 3 = cendre
allocate(monde(l,l)) ; allocate(stocha(l,l)) ; allocate(en_feu(l*l,2))
en_feu = -1
call init_random_seed()
call random_number(stocha)

! Création du monde, avec arbre et lac.
forall (i=1:l)
  forall (j=1:l)
    monde(i,j) = CEILING(stocha(i,j)-pa)
  end forall
end forall
deallocate(stocha)

!Flamme et CL
call random_number(x) ; call random_number(y)
monde(1,:) = 3 ; monde(l,:) = 3 ; monde(:,1) = 3 ; monde(:,l) = 3
monde(int(x*(l-2))+2,int(y*(l-2))+2) = 2
en_feu(1,:) = (/int(x*(l-2))+2, int(y*(l-2))+2/)
nb_en_feu = 1

allocate(en_feu2(l*l,2))
do while (nb_en_feu /= 0)
  
  k2 = nb_en_feu
  nb_en_feu = 0
  do k=1,k2
    do i=1,res
      do j=1,res
        if (monde(en_feu(k, 1)+nx(i), en_feu(k,2)+ny(j)) == 0 ) then
          monde(en_feu(k, 1)+nx(i), en_feu(k,2)+ny(j)) = 2
          nb_en_feu = nb_en_feu + 1
          en_feu2(nb_en_feu,:) = (/ en_feu(k, 1)+nx(i), en_feu(k,2)+ny(j)/)
        end if
      end do
    end do
    monde(en_feu(k, 1), en_feu(k,2)) = 3
  end do
  en_feu(:,:) = 0
  en_feu(:nb_en_feu,:) = en_feu2(:nb_en_feu,:)
end do
write(*,*) "Proportion d'arbre brulé : ", real(count(monde == 3)-4*(l-1))/((l-2.)*(l-2.))
deallocate(monde)
end program