final version (20230208)

source/DECISION.m

@@ -0,0 +1,163 @@
+function new_position_people = DECISION(position_to_be_updated,position_people,levels,walls,exit)
+%INPUT ... %position_to_be_updated ... náhodně vybraný index buňky, jejíž okolí budeme
+%aktualizovat v celé mřížce position_people...... NESMI TO BYT ZADNA Z
+%BUNEK NA OKRAJI
+%position_people ... tabulka 0 a 1 vyjadřující polohy lidí na
+%mřížce
+%levels ... tabulka potenciálů
+%walls ... tabulka stěn a sloupů
+%exit ... zvlast vypichnute souradnice vychodu
+%OUTPUT ... new_position ... aktualizovaná mřížka position_people
+
+new_position_people = position_people; %inicializace
+
+%v patem sloupci promenne neighborhood jsou hodnoty potenciálů, případně
+%NAN pro stěny, nebo sloupce
+neighborhood = GET_OKOLI(position_to_be_updated,position_people,levels,walls);
+
+if position_people(position_to_be_updated(1),position_to_be_updated(2)) == 1 %%%%%%nejdriv resim, pokud je ta vybrana bunka plna
+    min_pot = min(neighborhood(:,5));
+
+    %nalezeni radku s minimální hodnotou potencialu
+    indeces_of_min_pot = find(neighborhood(:,5)==min_pot);
+    %pocet radku s tim minimalnim potencialem
+    [n_min_pots,~] = size(indeces_of_min_pot);
+    %nejriv je snažíme tlacit na pole s minimalnim potencialem
+    is_there_a_man = 0;
+    is_there_a_wall = 0;
+    for ind = 1:n_min_pots
+        is_there_a_man(ind) = neighborhood(indeces_of_min_pot(ind),4);
+        is_there_a_wall(ind) = neighborhood(indeces_of_min_pot(ind),3);
+    end
+    if sum(is_there_a_man) < n_min_pots && sum(is_there_a_wall) < n_min_pots %pokud bude mozny vubec nejaky pohyb do vrstvy s nizsim potencialem
+        %pocet volnych bunek s nejnizsim potencialem
+        free_spots = n_min_pots - sum(is_there_a_man);
+        %nahodne vybere jednu v dostupnych volnych pozic
+        r = randi(free_spots);
+        %ted vyselektuju, kterou volnou pozici jsme vlastne vybrali
+        free_spot_position = find(is_there_a_man == 0);
+        selected_free_spot = free_spot_position(r);
+        selected_neigh_row = neighborhood(indeces_of_min_pot(selected_free_spot),:);
+
+        new_position_people(position_to_be_updated(1),position_to_be_updated(2)) = 0;
+        new_position_people(selected_neigh_row(1),selected_neigh_row(2)) = 1;
+
+
+    elseif  sum(is_there_a_man) == n_min_pots || sum(is_there_a_wall) == n_min_pots %pokud se nelze pohnout z duvodu preplnenosti nebo sten do vrstvy s nejnizsim potencialem
+        sec_min_pot = min(neighborhood(:,5))+1;
+        %nalezeni radku s druhou minimální hodnotou potencialu
+        indeces_of_sec_min_pot = find(neighborhood(:,5)==sec_min_pot);
+        %pocet radku s timto potencialem
+        [n_sec_min_pots,~] = size(indeces_of_sec_min_pot);
+        %nejdriv je snažíme tlacit na pole s minimalnim potencialem
+        is_there_a_man = 0;
+        is_there_a_wall = 0;
+        for ind = 1:n_sec_min_pots
+            is_there_a_man(ind) = neighborhood(indeces_of_sec_min_pot(ind),4);
+            is_there_a_wall(ind) = neighborhood(indeces_of_sec_min_pot(ind),3);
+        end
+        if sum(is_there_a_man) < n_sec_min_pots && sum(is_there_a_wall) < n_sec_min_pots %pokud bude mozny vubec nejaky pohyb
+            %pocet volnych bunek s nejnizsim potencialem
+            free_spots = n_sec_min_pots - sum(is_there_a_man);
+            %nahodne vybere jednu v dostupnych volnych pozic
+            r = randi(free_spots);
+            %ted vyselektuju, kterou volnou pozici jsme vlastne vybrali
+            free_spot_position = find(is_there_a_man == 0);
+            selected_free_spot = free_spot_position(r);
+            selected_neigh_row = neighborhood(indeces_of_sec_min_pot(selected_free_spot),:);
+
+            new_position_people(position_to_be_updated(1),position_to_be_updated(2)) = 0;
+            new_position_people(selected_neigh_row(1),selected_neigh_row(2)) = 1;
+        end
+
+    end
+else %aktualizace prazdne bunky je zalozena na tom, ze se tam prednostne tlaci spis lidi s vetsim potencialem a az potom ti na stejne urovni
+    if isnan(neighborhood(5,5)) %5. radek pole neighborhood odpovida bodu, ktery prave aktualizujeme, 5. sloupec pokud ma hodnotu NAN znamena, ze je tam zed/sloup
+        ;
+    else
+        max_pot = max(neighborhood(:,5));
+        %nalezeni radku s maximalni hodnotou potencialu
+        indeces_of_max_pot = find(neighborhood(:,5)==max_pot);
+        %pocet radku s tim maximalnim potencialem
+        [n_max_pots,~] = size(indeces_of_max_pot);
+        is_there_a_man = 0;
+        for ind = 1:n_max_pots
+            is_there_a_man(ind) = neighborhood(indeces_of_max_pot(ind),4);
+        end
+        if sum(is_there_a_man) > 0 %pokud bude mozny vubec nejaky pohyb do vrstvy s nizsim potencialem
+            %pocet plnych bunek s nejvyssim potencialem
+            full_spots = sum(is_there_a_man);
+            %nahodne vybere jednu z dostupnych plnych pozic
+            r = randi(full_spots);
+            %ted vyselektuju, kterou plnou pozici jsme vlastne vybrali
+            full_spot_position = find(is_there_a_man == 1);
+            selected_full_spot = full_spot_position(r);
+            selected_neigh_row = neighborhood(indeces_of_max_pot(selected_full_spot),:);
+
+            new_position_people(position_to_be_updated(1),position_to_be_updated(2)) = 1;
+            new_position_people(selected_neigh_row(1),selected_neigh_row(2)) = 0;
+        elseif sum(is_there_a_man) == 0 %z vrstvy s nejvyssim potencialem nema kdo sejit niz
+            sec_max_pot = max(neighborhood(:,5))-1;
+            %nalezeni radku s druhou nejvyssi hodnotou potencialu
+            indeces_of_sec_max_pot = find(neighborhood(:,5)==sec_max_pot);
+            %pocet radku s timto potencialem
+            [n_sec_max_pots,~] = size(indeces_of_sec_max_pot);
+            %nejdriv je snažíme tlacit na pole s minimalnim potencialem
+            is_there_a_man = 0;
+            for ind = 1:n_sec_max_pots
+                is_there_a_man(ind) = neighborhood(indeces_of_sec_max_pot(ind),4);
+            end
+            if sum(is_there_a_man) > 0 %pokud bude mozny vubec nejaky pohyb
+                %pocet plnych bunek s nejvyssim potencialem
+                full_spots = sum(is_there_a_man);
+                %nahodne vybere jednu z dostupnych plnych pozic
+                r = randi(full_spots);
+                %ted vyselektuju, kterou plnou pozici jsme vlastne vybrali
+                full_spot_position = find(is_there_a_man == 1);
+                selected_full_spot = full_spot_position(r);
+                selected_neigh_row = neighborhood(indeces_of_sec_max_pot(selected_full_spot),:);
+
+                new_position_people(position_to_be_updated(1),position_to_be_updated(2)) = 1;
+                new_position_people(selected_neigh_row(1),selected_neigh_row(2)) = 0;
+            elseif sum(is_there_a_man) == 0 %pokud ani z druheho nejvyssiho potencialu neni mozne sejit niz
+                thrd_max_pot = max(neighborhood(:,5))-2;
+                if thrd_max_pot >= neighborhood(5,5) %vezmu jeste treti mozny potencial, nesmi ale byt nizsi nez ten, ktery prislusi prazdne bunce, kterou se snazim aktualizovat
+                    %nalezeni radku s treti nejvyssi hodnotou potencialu
+                    indeces_of_thrd_max_pot = find(neighborhood(:,5)==thrd_max_pot);
+                    %pocet radku s timto potencialem
+                    [n_thrd_max_pots,~] = size(indeces_of_thrd_max_pot);
+                    %nejdriv je snažíme tlacit na pole s minimalnim potencialem
+                    is_there_a_man = 0;
+                    for ind = 1:n_thrd_max_pots
+                        is_there_a_man(ind) = neighborhood(indeces_of_thrd_max_pot(ind),4);
+                    end
+                    if sum(is_there_a_man) > 0 %pokud bude mozny vubec nejaky pohyb
+                        %pocet plnych bunek s nejvyssim potencialem
+                        full_spots = sum(is_there_a_man);
+                        %nahodne vybere jednu z dostupnych plnych pozic
+                        r = randi(full_spots);
+                        %ted vyselektuju, kterou plnou pozici jsme vlastne vybrali
+                        full_spot_position = find(is_there_a_man == 1);
+                        selected_full_spot = full_spot_position(r);
+                        selected_neigh_row = neighborhood(indeces_of_thrd_max_pot(selected_full_spot),:);
+
+                        new_position_people(position_to_be_updated(1),position_to_be_updated(2)) = 1;
+                        new_position_people(selected_neigh_row(1),selected_neigh_row(2)) = 0;
+                    end
+                end
+            end
+        end
+
+
+
+    end
+
+end
+%pokud bude plna bunka se souradnicemi "exit", pred vykreslenim ji vzdy
+%vyprazdnim
+
+if new_position_people(exit(1),exit(2)) == 1
+    new_position_people(exit(1),exit(2)) = 0;
+end
+
+end

source/Floor_field_model.m

@@ -0,0 +1,104 @@
+%% FLOOR FIELD MODEL - spousteci skript
+
+%% uklid
+close all
+clear variables
+
+%% Nastaveni konstant (pomoci fce 'params.m')
+% lze vytvorit nekolik ruznych sad pocatecnich parametru, treba fce
+% 'params1.m', 'params2.m' atd)
+
+params = parameters(4);
+P = params.P;
+Q = params.Q;
+loc_of_exit = params.exit;
+mtx_loc_of_pillars = params.pillars;
+num_of_peds = params.pedestrians;
+
+%[num_of_peds, P, Q, mtx_loc_of_pillars, loc_of_exit, animation_show] = params();
+%[num_of_peds, P, Q, mtx_loc_of_pillars, loc_of_exit, animation_show] = params_big_room_full();  % pozor: Dijkstra bezi cca 4min na tak velke mape
+%[num_of_peds, P, Q, mtx_loc_of_pillars, loc_of_exit, animation_show] = params_big_room_empty(); % pozor: Dijkstra bezi cca 4min na tak velke mape
+
+animation_show = true;
+%% Nastaveni mapy, rozmisteni sloupu a sten
+walls = double(isnan(get_map (P, Q, loc_of_exit, mtx_loc_of_pillars)));  % AG: melo by to vratit mapu, kde na vsech "nevkrocitelnych" polich je 1, a na volnych 0
+% nevim, kde presne se to pouzije, ale jestli je potreba, tak tu je.
+
+%% -------------------------- simulace FFM
+n_simulations = 30;                         % počet simulácií
+iterations    = ones(n_simulations, 1);     % počet krokov do úplnej evakuácie miestnosti
+rho           = [];                         % prázdne pole pre ukladanie hustôt
+
+rng(1)                                      % pevný seed
+
+for i = 1:n_simulations
+    %% Vypocet potencialu na mape (fce 'get_grad_field.m' od AG)
+    tStart = tic;
+    levels = get_grad_field(P, Q, loc_of_exit, mtx_loc_of_pillars);  % Dijkstra runs inside
+    tEnd = toc(tStart)
+    'log: get_grad_field.m successfully done'
+
+    %% Umisteni chodcu do mistnosti (fce 'get_pedestrians.m' od SM)
+    position_people = get_pedestrians(P, Q, mtx_loc_of_pillars, num_of_peds);
+    'log: get_pedestrians.m successfully done'
+
+    %% Pomocna matice pro vykresleni polohy exitu zelenou barvou
+    exit_matrix = 0*position_people; %stejne rozmery jako matice position people
+    exit_matrix(loc_of_exit(1),loc_of_exit(2)) = 1;
+
+    while (sum(sum(position_people)) ~= 0)
+        [row_num,col_num] = size(position_people);
+        %budu se snazit generovat cisla pouze z te zaplnene casti, abych to
+        %urychlila
+        first_nonempty_col = 1;
+        while sum(position_people(:,first_nonempty_col)) == 0 && (first_nonempty_col + 1) <= (col_num - 1)
+            first_nonempty_col = first_nonempty_col + 1;
+        end
+
+        %ind_row a ind_col budou indexy bunky, kterou prave aktualizuju, nesmi
+        %to byt nic z okraju pole
+        ind_row = randi([2,row_num-1]);
+        ind_col = randi([first_nonempty_col,col_num-1]);
+
+        position_people = DECISION([ind_row,ind_col],position_people,levels,walls,loc_of_exit);
+        R = 255 * (-position_people+1) - 255*exit_matrix;
+        G = 255 * (-position_people+1)-255*walls;
+        B = 255 * (-position_people+1)-255*walls - 255*exit_matrix;
+        res = cat(3,R,G,B);
+
+%         %plot result
+%              if animation_show == true
+%                  imshow(res,'InitialMagnification',20000)
+%                  pause(0.07);
+%              end
+
+        % Hustota chodcov v miestnosti v danej iterácii
+        rho = [rho; stats(P, Q, mtx_loc_of_pillars, position_people)];
+
+        iterations(i) = iterations(i)+1;
+    end
+end
+
+% Stredný počet krokov potrebný pre evakuáciu miestnosti
+[mu, delta] = expectation(iterations, n_simulations);
+histogram(iterations)
+xlabel("number of steps",'interpreter','latex')
+ylabel("frequency",'interpreter','latex')
+set(gca,'TickLabelInterpreter','latex')
+
+'Done.'
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

source/GET_OKOLI.m

@@ -0,0 +1,68 @@
+function okoli = GET_OKOLI(ped,pozice_lidi,urovne,zdi)
+%ped jsou souradnice bunky, kterou jsem se rozhodla aktualizovat, nesmi to
+%byt z kraje mrizky, tam jsou proste steny
+okoli = nan(9,5);
+okoli(1,1) = ped(1)-1;
+okoli(1,2) = ped(2)-1;
+okoli(1,3) = zdi(okoli(1,1),okoli(1,2));
+okoli(1,4) = pozice_lidi(okoli(1,1),okoli(1,2));
+okoli(1,5) = urovne(okoli(1,1),okoli(1,2));
+
+
+okoli(2,1) = ped(1)-1;
+okoli(2,2) = ped(2);
+okoli(2,3) = zdi(okoli(2,1),okoli(2,2));
+okoli(2,4) = pozice_lidi(okoli(2,1),okoli(2,2));
+okoli(2,5) = urovne(okoli(2,1),okoli(2,2));
+
+
+okoli(3,1) = ped(1)-1;
+okoli(3,2) = ped(2)+1;
+okoli(3,3) = zdi(okoli(3,1),okoli(3,2));
+okoli(3,4) = pozice_lidi(okoli(3,1),okoli(3,2));
+okoli(3,5) = urovne(okoli(3,1),okoli(3,2));
+
+
+okoli(4,1) = ped(1);
+okoli(4,2) = ped(2)-1;
+okoli(4,3) = zdi(okoli(4,1),okoli(4,2));
+okoli(4,4) = pozice_lidi(okoli(4,1),okoli(4,2));
+okoli(4,5) = urovne(okoli(4,1),okoli(4,2));
+
+
+okoli(5,1) = ped(1);
+okoli(5,2) = ped(2);
+okoli(5,3) = zdi(okoli(5,1),okoli(5,2));
+okoli(5,4) = pozice_lidi(okoli(5,1),okoli(5,2));
+okoli(5,5) = urovne(okoli(5,1),okoli(5,2));
+
+
+okoli(6,1) = ped(1);
+okoli(6,2) = ped(2)+1;
+okoli(6,3) = zdi(okoli(6,1),okoli(6,2));
+okoli(6,4) = pozice_lidi(okoli(6,1),okoli(6,2));
+okoli(6,5) = urovne(okoli(6,1),okoli(6,2));
+
+
+okoli(7,1) = ped(1)+1;
+okoli(7,2) = ped(2)-1;
+okoli(7,3) = zdi(okoli(7,1),okoli(7,2));
+okoli(7,4) = pozice_lidi(okoli(7,1),okoli(7,2));
+okoli(7,5) = urovne(okoli(7,1),okoli(7,2));
+
+
+okoli(8,1) = ped(1)+1;
+okoli(8,2) = ped(2);
+okoli(8,3) = zdi(okoli(8,1),okoli(8,2));
+okoli(8,4) = pozice_lidi(okoli(8,1),okoli(8,2));
+okoli(8,5) = urovne(okoli(8,1),okoli(8,2));
+
+
+okoli(9,1) = ped(1)+1;
+okoli(9,2) = ped(2)+1;
+okoli(9,3) = zdi(okoli(9,1),okoli(9,2));
+okoli(9,4) = pozice_lidi(okoli(9,1),okoli(9,2));
+okoli(9,5) = urovne(okoli(9,1),okoli(9,2));
+
+end
+