Make a microscopic simulation

In the two following models, individuals are seen as spherical particles. They move at a desired velocity to reach a goal (typically a door). Doing nothing more would lead to many overlaps between individuals. Hence the two families of models below can prevent these overlaps through two different approaches. For the Social Force models, and, for some forces are added to act against overlaps, and for the Granular model the velocities are projected into a permissible velocity space which ensures the absence of overlaps.

Social force model

The Social Force model has been introduced in the 90’s. Pedestrians are identified with inertial particles submitted to a forcing term which implements the individual tendencies and extra forces which account for interactions with other pedestrians (typically the tendency to preserve a certain distance with neighbors).

Reference : [MF2018] Chapter 3.

Some examples can be found in the directory

cromosim/examples/micro/social

and can be launched with

python micro_social.py --json input_room.json
python micro_social.py --json input_event.json
python micro_social.py --json input_stadium.json
python micro_social.py --json input_shibuya_crossing.json
python micro_social.py --json input_stairs.json
Evacuation of a room, visualization of trajectories

Circular track around a stadium

Evacuation of an exhibition hall: two groups of people with the same destination
Evacuation of an exhibition hall: sensors (green lines) results
Evacuation of an exhibition hall : Sensors (green lines) results

Shibuya crossing (Japan) : five groups of people and five different destinations

Two floors of a building: a group of people goes up and another goes down
Floor 1 on the left and 0 on the right
micro_social.py
  1# Authors:
  2#     Sylvain Faure <sylvain.faure@universite-paris-saclay.fr>
  3#     Bertrand Maury <bertrand.maury@universite-paris-saclay.fr>
  4#
  5#     cromosim/examples/micro/granular/micro_granular.py
  6#     python micro_granular.py --json input.json
  7#
  8# License: GPL
  9
 10import sys
 11import os
 12import json
 13import time
 14from optparse import OptionParser
 15import numpy as np
 16import matplotlib.pyplot as plt
 17from matplotlib.patches import Circle, Ellipse, Rectangle, Polygon
 18from matplotlib.lines import Line2D
 19
 20import logging
 21logging.basicConfig(
 22    level=logging.INFO,
 23    # format="%(asctime)s | %(levelname)-8s | %(name)s | %(message)s",
 24    format="%(levelname)-8s | %(name)s | %(message)s",
 25)
 26LOGGER = logging.getLogger(__name__)
 27# LOG_FORMAT = "%(asctime)s | %(levelname)-8s | %(name)s | %(message)s"
 28
 29from cromosim.domain import Destination, Domain
 30from cromosim.micro import (
 31    people_initialization,
 32    add_people_flows,
 33    append_people,
 34    compute_forces,
 35    compute_contacts,
 36    find_duplicate_people,
 37    move_people,
 38    people_update_destination,
 39    projection,
 40)
 41from cromosim.visu import (
 42    plot_sensors
 43)
 44
 45"""
 46    python3 micro_social.py --json input.json
 47"""
 48
 49parser = OptionParser(usage="usage: %prog [options] filename",
 50                      version="%prog 1.0")
 51parser.add_option('--json', dest="jsonfilename", default="input.json",
 52                  type="string", action="store", help="Input json filename")
 53opt, remainder = parser.parse_args()
 54
 55LOGGER.info(f"Loading JSON configuration: {opt.jsonfilename}")
 56
 57with open(opt.jsonfilename) as json_file:
 58    try:
 59        config = json.load(json_file)
 60    except json.JSONDecodeError as exc:
 61        raise ValueError(f"Invalid JSON in {opt.jsonfilename}: {exc}") from exc
 62
 63LOGGER.info(
 64    "Simulation: Tf=%s, dt=%s, projection_method=%s",
 65    config["Tf"],
 66    config["dt"],
 67    config["projection_method"],
 68)
 69LOGGER.info(
 70    "Contacts: dmax=%s, dmin_people=%s, dmin_walls=%s",
 71    config["dmax"],
 72    config["dmin_people"],
 73    config["dmin_walls"],
 74)
 75LOGGER.info(
 76    "Visualization: enabled=%s, backend=%s, people=%s " 
 77    "contacts=%s velocities=%s desired_velocities=%s "
 78    "paths=%s sensors=%s plot_period=%s",
 79    config["with_graphes"],
 80    config["plot_module"],
 81    config["plot_people"],
 82    config["plot_contacts"],
 83    config["plot_velocities"],
 84    config["plot_desired_velocities"],
 85    config["plot_paths"],
 86    config["plot_sensors"],
 87    config["plot_period"],
 88)
 89LOGGER.info(
 90    "Miscellaneous: prefix=%s, seed=%s, dmin_walls=%s",
 91    config["prefix"],
 92    config["seed"],
 93    config["dmin_walls"],
 94)
 95
 96os.makedirs(config["prefix"], exist_ok=True)
 97
 98"""
 99    Build the Domain objects
100"""
101
102domains = {}
103
104domain_configs = config["domains"]
105
106for domain_index, domain_config in enumerate(domain_configs):
107    name = domain_config["name"]
108    LOGGER.info(f"Start building domain {domain_index}: {name}")
109    if domain_config["background"] == "":
110        domain = Domain(
111            name=name,
112            pixel_size=domain_config["px"],
113            width=domain_config["width"],
114            height=domain_config["height"],
115            wall_colors=domain_config["wall_colors"],
116        )
117    else:
118        domain = Domain(
119            name=name,
120            background=domain_config["background"],
121            pixel_size=domain_config["px"],
122            wall_colors=domain_config["wall_colors"],
123        )
124        
125    LOGGER.info(
126        f"Building domain {domain_index}: {name} => "
127        f"Add geometric objects described in the JSON "
128        f"configuration")
129    
130    for shape_config in domain_config["shape_lines"]:
131        shape = Line2D(
132            shape_config["xx"], shape_config["yy"],
133            linewidth=shape_config["linewidth"],
134        )
135        domain.add_shape(
136            shape,
137            outline_color=shape_config["outline_color"],
138            fill_color=shape_config["fill_color"],
139        )
140
141    for shape_config in domain_config["shape_circles"]:
142        shape = Circle(
143            (shape_config["center_x"], shape_config["center_y"]),
144            shape_config["radius"],
145        )
146        domain.add_shape(
147            shape,
148            outline_color=shape_config["outline_color"],
149            fill_color=shape_config["fill_color"],
150        )
151
152    for shape_config in domain_config["shape_ellipses"]:
153        shape = Ellipse(
154            (shape_config["center_x"], shape_config["center_y"]),
155            shape_config["width"], shape_config["height"],
156            shape_config["angle_in_degrees_anti-clockwise"],
157        )
158        domain.add_shape(
159            shape,
160            outline_color=shape_config["outline_color"],
161            fill_color=shape_config["fill_color"],
162        )
163
164    for shape_config in domain_config["shape_rectangles"]:
165        shape = Rectangle(
166            (shape_config["bottom_left_x"], shape_config["bottom_left_y"]),
167            shape_config["width"], shape_config["height"],
168            shape_config["angle_in_degrees_anti-clockwise"],
169        )
170        domain.add_shape(
171            shape,
172            outline_color=shape_config["outline_color"],
173            fill_color=shape_config["fill_color"],
174        )
175
176    for shape_config in domain_config["shape_polygons"]:
177        shape = Polygon(shape_config["xy"])
178        domain.add_shape(
179            shape,
180            outline_color=shape_config["outline_color"],
181            fill_color=shape_config["fill_color"],
182        )
183        
184    LOGGER.info(f"Building domain {domain_index}: {name} => build_domain()")    
185    domain.build_domain()
186    
187    LOGGER.info(
188        f"Building domain {domain_index}: {name} => "
189        f"Add the destinations described in the JSON "
190        f"configuration")
191
192    for destination_index, destination_config in enumerate(
193            domain_config["destinations"]
194        ):
195        desired_velocity_from_color = [
196            np.concatenate((cdv["color"], cdv["desired_velocity"]))
197            for cdv in destination_config["desired_velocity_from_color"]
198        ]
199        destination = Destination(
200            name=destination_config["name"], 
201            colors=destination_config["colors"],
202            excluded_colors=destination_config["excluded_colors"],
203            desired_velocity_from_color=desired_velocity_from_color,
204            velocity_scale=destination_config["velocity_scale"],
205            next_destination=destination_config["next_destination"],
206            next_domain=destination_config["next_domain"],
207            next_destination_weights=destination_config["next_destination_weights"],
208            next_transit_box=destination_config["next_transit_box"],
209            next_transit_box_xshift=destination_config["next_transit_box_xshift"],
210            next_transit_box_yshift=destination_config["next_transit_box_yshift"],
211        )
212        LOGGER.info(
213                f"Building domain {domain_index}: {name} => "
214                f"Add destination {destination_index}: {destination_config["name"]}")
215        domain.add_destination(destination)
216        if (config["with_graphes"]):
217            domain.plot_desired_velocity(
218                destination_config["name"], 
219                id=100*domain_index+10+destination_index, 
220                step=20,
221                title=f"Destination '{destination_config["name"]}' "
222                      f"in '{name}' domain")
223     
224    if (config["with_graphes"]):
225        domain.plot(
226            id=100*domain_index, 
227            title=f"Geometry of the '{name}' domain")
228        domain.plot_wall_dist(
229            id=100*domain_index+1, 
230            step=20,
231            title=f"Distance from walls for the '{name}' domain")
232       
233    domains[name] = domain
234
235
236"""
237    To create the sensors to measure the pedestrian flows
238"""
239
240all_sensors = {domain_name: [] for domain_name in domains}
241
242for sensor_config in config["sensors"]:
243    sensor = dict(sensor_config)
244    sensor["id"] = []
245    sensor["times"] = []
246    sensor["xy"] = []
247    sensor["dir"] = []
248    LOGGER.info(
249        f"Initializing the '{sensor["name"]}' sensor "
250        f"located in the '{sensor["domain"]}' domain")
251    all_sensors[sensor["domain"]].append(sensor)
252
253
254"""
255    Initialize pedestrian populations in all configured domains
256"""
257
258# Current time
259t = 0.0
260counter = 0
261dt = config["dt"]
262Tf = config["Tf"]
263
264all_people = {}
265
266for people_index, people_config in enumerate(config["people_init"]):
267    domain_name = people_config["domain"]
268    domain = domains[domain_name]
269    LOGGER.info(
270        f"Initializing people configuration in the '{domain_name}' domain"
271    )
272
273    people = people_initialization(
274        domain,
275        people_config["groups"],
276        dt,
277        dmin_people=config["dmin_people"],
278        dmin_walls=config["dmin_walls"],
279        seed=config["seed"],
280        itermax=10,
281        projection_method="uzawa",
282        flow_duration=60,
283        log=True,
284    )
285    if people is None:
286        raise RuntimeError(
287            f"No population was created for domain {domain_name!r}."
288        )
289    else:
290        all_people[domain_name] = people
291
292
293"""
294    Visualization: either using Matplotlib or pyvista 
295"""
296
297all_visu = {} # one per domain...
298
299for domain_index, (name, domain) in enumerate(domains.items()):
300
301    LOGGER.info(
302        f"Creating the view vindow for the '{domain_name}' domain "
303        f"(render={config["plot_module"]})"
304    )
305    
306    if not config["with_graphes"]:
307    
308        all_visu[name] = None
309    
310    else:
311
312        people = all_people[name]
313        colors = people["xyrv"][:, 2]
314        if colors.size>0:
315            cmin, cmax = colors.min(), colors.max()
316        else:
317            cmin, cmax = 0, 1
318                
319        II, JJ, Vd = domain.people_desired_velocity(
320            people["xyrv"], 
321            people["destinations"]
322        )
323        people["Vd"] = Vd
324                
325        contacts = compute_contacts(domain, people["xyrv"], config["dmax"])      
326
327        if config["plot_module"] == "matplotlib":
328            
329            from cromosim.visu import matplotlibVisualization
330            
331            visu = matplotlibVisualization(
332                domain,
333                figure_id=100 * (1+domain_index) + 20,
334                cmap="viridis",
335                velocity_reference=config["plot_velocity_reference"],
336                velocity_arrow_length=config["plot_velocity_arrow_length"]
337            )
338            visu.update(
339                people=people,
340                contacts=contacts,
341                sensors=all_sensors[name],
342                time=t,
343                plot_people=config["plot_people"],
344                plot_contacts=config["plot_contacts"],
345                plot_paths=config["plot_paths"],
346                plot_velocities=config["plot_velocities"],
347                plot_desired_velocities=config["plot_desired_velocities"],
348                plot_sensors=config["plot_sensors"],
349                color_scalar_values=colors,
350                color_scalar_name="r",
351                color_minmax=(cmin, cmax),
352                color_cmap="viridis",
353                show_scalar_bar=config["plot_people"],
354                savefig=True,
355                filename=os.path.join(
356                    config["prefix"],
357                    f"{name}_fig_{counter:06d}.png",
358                ),
359            )
360            plt.pause(0.001)
361            visu.show()
362            all_visu[name] = visu
363
364        elif config["plot_module"] == "pyvista":
365        
366            from cromosim.visu import pyvistaVisualization
367        
368            target_width = 1920
369            target_height = round(
370                target_width * domain.height / domain.width
371            )    
372            # scale = max(1, round(1920 / domain.width))
373            movie_path = os.path.join(config["prefix"], f"{name}.mp4")
374            visu = pyvistaVisualization(
375                domain,
376                image_path=domain._Domain__image_filename,
377                pixel_size=domain.pixel_size,
378                movie_path=movie_path,
379                fps=4, # max(1,int(round(1/(config["plot_period"]*dt)))),
380                # window_size=(scale * domain.width, scale * domain.height),
381                window_size=(target_width, target_height),
382                movie_macro_block_size=16,
383                velocity_reference=config["plot_velocity_reference"],
384                velocity_arrow_length=config["plot_velocity_arrow_length"]
385            )
386            visu.update(
387                people=people,
388                contacts=contacts,
389                sensors=all_sensors[name],
390                time=t,
391                plot_people=config["plot_people"],
392                plot_contacts=config["plot_contacts"],
393                plot_velocities=config["plot_velocities"],
394                plot_desired_velocities=config["plot_desired_velocities"],
395                plot_paths=config["plot_paths"],
396                plot_sensors=config["plot_sensors"],
397                color_scalar_values=colors,
398                color_scalar_name="r",
399                color_minmax=(cmin, cmax),
400                color_cmap="viridis",
401                show_scalar_bar=config["plot_people"],
402                write_frame=True,
403            )
404            visu.show()
405            all_visu[name] = visu
406                
407
408# It may be helpful to pause here to check the calculation 
409# domains and the positions of the individuals before 
410# continuing the simulation...
411# sys.exit()
412
413"""
414    Main loop
415"""
416
417cc = 0
418draw = True
419
420while (t < Tf):
421    
422    LOGGER.info( f"-------------------------------------\n" )
423    LOGGER.info( f"===> Time = {t}" )
424
425    # If there are inflows of people arriving in the area, we calculate 
426    # the number of people to be added.
427    
428    for people_index, people_config in enumerate(config["people_init"]):
429        domain_name = people_config["domain"]
430        domain = domains[domain_name]
431        people_to_add = add_people_flows(
432            domain,
433            people_config["groups"],
434            dt,
435            dmin_people=config["dmin_people"],
436            dmin_walls=config["dmin_walls"],
437            rng=all_people[domain_name]["rng"],
438            itermax=10,
439            projection_method="uzawa",
440            flow_duration=config["dt"],
441            log=True,
442        )
443        if people_to_add is not None: 
444            LOGGER.info(
445                f"People flow in the '{domain_name}' domain: "
446                f"number of people added = {people_to_add["xyrv"].shape[0]}"
447            )
448            append_people(all_people[domain_name], people_to_add)
449                            
450 
451     # Compute people desired velocity
452     
453    for domain_index, (name, domain) in enumerate(domains.items()):
454        LOGGER.info(
455            f"Compute people desired velocity in the '{name}' domain"
456            )
457        people = all_people[name]
458        II, JJ, Vd = domain.people_desired_velocity(
459            people["xyrv"],
460            people["destinations"])
461        people["Vd"][:,0] = people["beta"]*Vd[:,0]
462        people["Vd"][:,1] = people["beta"]*Vd[:,1]
463        people["I"] = II
464        people["J"] = JJ
465
466    # Look at if there are people in the transit boxes
467    
468    LOGGER.info(
469                f"Find the people who need to be duplicated when "
470                f"they are in an area where two domains overlap "
471                f"(such as when changing floors)"
472            )
473    virtual_people = find_duplicate_people(all_people, domains)
474 
475    # Social forces
476    
477    for domain_index, (name, domain) in enumerate(domains.items()):
478        LOGGER.info(
479                    f"Start to compute social forces for the '{name}' domain "
480                    f"to prevent individuals from overlapping..."
481                )
482        people = all_people[name]
483
484        if virtual_people[name]["xyrv"] is not None:
485            LOGGER.info(
486                        f"Projection step for the '{name}' domain: "
487                        f"{virtual_people[name]["xyrv"].shape[0]} "
488                        f"virtual guys found"
489                    )
490            xyrv = np.concatenate(
491                (people["xyrv"], virtual_people[name]["xyrv"]))
492            Vd = np.concatenate(
493                (people["Vd"], virtual_people[name]["Vd"]))
494            Uold = np.concatenate(
495                (people["Uold"], virtual_people[name]["Uold"]))
496            Mass = np.concatenate(
497                (people["mass"], virtual_people[name]["mass"]))
498        else:
499            xyrv = people["xyrv"]
500            Vd = people["Vd"]
501            Uold = people["Uold"]
502            Mass = people["mass"]
503
504        if (xyrv.shape[0] > 0):
505
506            if (np.unique(xyrv, axis=0).shape[0] != xyrv.shape[0]):
507                raise ValueError( 
508                                 "There are two identical lines in "
509                                 "the array xyrv used to determine the "
510                                 "contacts between individuals and this "
511                                 "is not normal!"
512                            )
513                         
514            contacts = compute_contacts(domain, xyrv, config["dmax"])
515            LOGGER.info(
516                f"Projection step for the '{name}' domain: "
517                f"{contacts.shape[0]} contacts computed"
518            )
519
520            Forces = compute_forces(
521                config["F"], 
522                config["Fwall"], 
523                xyrv, 
524                contacts, 
525                Uold, 
526                Vd,
527                config["lambda"], 
528                config["delta"], 
529                config["kappa"], 
530                config["eta"]
531            )
532
533            nn = people["xyrv"].shape[0]
534            
535            all_people[name]["U"] = dt * (Vd[:nn, :]-Uold[:nn, :]) / config["tau"] + \
536                            Uold[:nn, :] + dt*Forces[:nn, :]/Mass[:nn, None]
537            
538            virtual_people[name]["U"] = dt * (Vd[nn:, :]-Uold[nn:, :]) / config["tau"] +\
539                            Uold[nn:, :] + dt*Forces[nn:, :]/Mass[nn:, None]
540                            
541            # We retain the last frustration values and compute beta 
542               
543            # Frustration is 1 - <U_i, Vd_i> / <Vd_i, Vd_i>, clipped to [0, 1].
544            
545            satisfactions = np.minimum(
546                np.maximum(
547                    (Vd[:nn, 0] * all_people[name]["U"][:, 0] + Vd[:nn, 1] * all_people[name]["U"][:, 1])
548                    / np.maximum(Vd[:nn, 0] ** 2 + Vd[:nn, 1] ** 2, 1.0e-8),
549                    0.0,
550                ),
551                1.0,
552            )
553            frustrations = 1 - satisfactions
554                 
555            all_people[name]["frustration_last_pos"] = (all_people[name]["frustration_last_pos"] + 1) % \
556                config["frustration_memory_size"]
557            
558            all_people[name]["frustration_memory"      ][ np.arange(nn), all_people[name]["frustration_last_pos"] ] = frustrations[:nn]
559            all_people[name]["frustration_memory_times"][ np.arange(nn), all_people[name]["frustration_last_pos"] ] = t
560                        
561            # The Vd correcting factor, beta, is activated only 
562            # for individual who have neighbors in their cone
563            # of vision and at a distance less than a disc diameter.
564            
565            all_people[name]["beta"] = np.ones((nn,))
566            
567            active = np.zeros((nn,))
568            for ic, cont in enumerate(contacts):
569                i = int(cont[0])
570                j = int(cont[1])
571                dij = cont[2]
572                eij_x = cont[3]
573                eij_y = cont[4]
574                if (j > -1):  # contact person/person
575                    # theta_ij is the angle at which i (assumed to look 
576                    # ahead along their desired velocity) sees j
577                    # Angular dependence
578                    norm_Vdi = np.sqrt(Vd[i, 0]**2+Vd[i, 1]**2)
579                    if (norm_Vdi > 0) and (i<nn): # not virtual people...
580                        theta_ij = np.arccos( max(-1.0, min(1.0, ( (Vd[i,0]*eij_x+Vd[i,1]*eij_y)/norm_Vdi ))))
581                        if theta_ij<=all_people[name]["viewing_angle"][i]/180.0*np.pi:
582                            active[i] = 1
583                                                    
584                    norm_Vdj = np.sqrt(Vd[j, 0]**2+Vd[j, 1]**2)
585                    if (norm_Vdj > 0) and (j<nn): # not virtual people...
586                        theta_ji = np.arccos( max(-1.0, min(1.0, ( -(Vd[j,0]*eij_x+Vd[j,1]*eij_y)/norm_Vdj ))))
587                        if theta_ji<=all_people[name]["viewing_angle"][j]/180.0*np.pi:
588                            active[j] = 1
589
590            ind_frust = np.where( (all_people[name]["with_frustration"]==1) * \
591                (active==1) )[0]
592                                                                
593            all_people[name]["beta"][ind_frust] = 1.0 - \
594                np.mean(all_people[name]["frustration_memory"], axis=1)[ind_frust]
595                                    
596            # Update people positions
597            
598            all_people[name], all_sensors[name] = move_people(
599                t, dt,
600                all_people[name],
601                all_sensors[name])
602
603    # Update people destinations
604    
605    # the distance at which one changes destinations:
606    dist_to_change_dest = 2*dt
607    all_people = people_update_destination(
608        all_people, domains, 
609        dist_to_change_dest)
610
611    # Update previous velocities
612    for idom, name in enumerate(domains):
613        all_people[name]["Uold"] = all_people[name]["U"]
614  
615    t += dt
616    cc += 1
617    counter += 1
618  
619    for domain_index, (name, domain) in enumerate(domains.items()):
620        LOGGER.info(
621            f"--- Number of persons in the '{name}' "
622            f"domain: {all_people[name]["xyrv"].shape[0]}"
623        )    
624
625    if (cc >= config["plot_period"]):
626        draw = True
627        cc = 0
628    else:
629        draw = False
630    
631    # Ntot = sum([people["xyrv"].shape[0] for (name,people) in all_people.items()])
632    # if (Ntot==0):
633    #     LOGGER.info(
634    #         f"No one left in any domain! END..."
635    #     )
636    #     break    
637    
638    """
639    Visualization: either using Matplotlib or pyvista 
640    """
641
642    if draw and config["with_graphes"]:
643
644        for domain_index, (name, domain) in enumerate(domains.items()):
645
646            people = all_people[name]
647            
648            colors = people["xyrv"][:, 2]
649            if colors.size>0:
650                cmin, cmax = colors.min(), colors.max()
651            else:
652                cmin, cmax = 0, 1
653            
654            # Since people have moved, we recalculate Vd and 
655            # the contacts to ensure a correct display
656            II, JJ, Vd = domain.people_desired_velocity(
657                people["xyrv"], 
658                people["destinations"]
659            )
660            people["Vd"][:,0] = people["beta"]*Vd[:,0]
661            people["Vd"][:,1] = people["beta"]*Vd[:,1]        
662            
663            contacts = compute_contacts(domain, people["xyrv"], config["dmax"]) 
664            
665            if config["plot_module"] == "matplotlib":
666            
667                all_visu[name].update(
668                    people=people,
669                    contacts=contacts,
670                    virtual_people=virtual_people[name],
671                    sensors=all_sensors[name],
672                    time=t,
673                    plot_people=config["plot_people"],
674                    plot_contacts=config["plot_contacts"],
675                    plot_paths=config["plot_paths"],
676                    plot_velocities=config["plot_velocities"],
677                    plot_desired_velocities=config["plot_desired_velocities"],
678                    plot_sensors=config["plot_sensors"],
679                    color_scalar_values=colors,
680                    color_scalar_name="r",
681                    color_minmax=(cmin, cmax),
682                    color_cmap="viridis",
683                    show_scalar_bar=config["plot_people"],
684                    savefig=True,
685                    filename=os.path.join(
686                        config["prefix"],
687                        f"{name}_fig_{counter:06d}.png",
688                    ),
689                )
690                plt.pause(0.001)
691                all_visu[name].show()
692
693            elif config["plot_module"] == "pyvista":
694                       
695                all_visu[name].update(
696                    people=people,
697                    contacts=contacts,
698                    virtual_people=virtual_people[name],
699                    sensors=all_sensors[name],
700                    time=t,
701                    plot_people=config["plot_people"],
702                    plot_contacts=config["plot_contacts"],
703                    plot_velocities=config["plot_velocities"],
704                    plot_desired_velocities=config["plot_desired_velocities"],
705                    plot_paths=config["plot_paths"],
706                    plot_sensors=config["plot_sensors"],
707                    color_scalar_values=colors,
708                    color_scalar_name="r",
709                    color_minmax=(cmin, cmax),
710                    color_cmap="viridis",
711                    show_scalar_bar=config["plot_people"],
712                    write_frame=True,
713                )
714                all_visu[name].process_events()
715                all_visu[name].plotter.update()
716                all_visu[name].show()      
717
718
719# End...
720
721LOGGER.info(
722    f"============ It's over! ============\n"
723)    
724
725
726if config["with_graphes"]:
727    
728    for domain_index, (name, domain) in enumerate(domains.items()):
729        
730        if config["plot_module"] == "pyvista":
731            LOGGER.info(
732                f"Close the pyvista visualization window for the '{name}' domain"
733            )    
734            all_visu[name].close()
735        
736        if config["plot_sensors"]:
737            sensors = all_sensors[domain_name]
738            LOGGER.info(
739                f"Plot sensor data for the '{name}' domain "
740                f"({len(sensors)} sensors)"
741            )    
742            plot_sensors(
743                100*domain_index+40, 
744                sensors, 
745                t, 
746                savefig=True,
747                filename=os.path.join(
748                    config["prefix"],
749                    f"sensors_{domain_name}.png"
750                    )
751                )
752            plt.pause(0.001)
753
754plt.show()

Granular model

The Granular model comes from crowd motion models of the granular type : each individual is identified to a hard disk of a prescribed size, subject to a non-overlapping constraint with their neighbors. The approach relies on a desired velocity for each individual (the velocity they would take if they were alone), and the global velocity field shall be defined as the closest to the desired one among all those feasible fields (fields which do not lead to overlapping of disks).

Reference : [MF2018] Chapter 4.

An example can be find in the directory

cromosim/examples/micro/granular

and can be launched with

python micro_granular.py --json input_room.json
python micro_granular.py --json input_event.json
python micro_granular.py --json input_stadium.json
python micro_granular.py --json input_shibuya_crossing.json
python micro_granular.py --json input_stairs.json
Evacuation of a room, visualization of trajectories

Circular track around a stadium

Evacuation of an exhibition hall: two groups of people with the same destination
Evacuation of an exhibition hall: sensors (green lines) results
Evacuation of an exhibition hall : Sensors (green lines) results

Shibuya crossing (Japan) : five groups of people and five different destinations

Two floors of a building: a group of people goes up and another goes down
Floor 1 on the left and 0 on the right
micro_granular.py
  1# Authors:
  2#     Sylvain Faure <sylvain.faure@universite-paris-saclay.fr>
  3#     Bertrand Maury <bertrand.maury@universite-paris-saclay.fr>
  4#
  5#     cromosim/examples/micro/granular/micro_granular.py
  6#     python micro_granular.py --json input.json
  7#
  8# License: GPL
  9
 10import sys
 11import os
 12import json
 13import time
 14from optparse import OptionParser
 15import numpy as np
 16import matplotlib.pyplot as plt
 17from matplotlib.patches import Circle, Ellipse, Rectangle, Polygon
 18from matplotlib.lines import Line2D
 19
 20import logging
 21logging.basicConfig(
 22    level=logging.INFO,
 23    # format="%(asctime)s | %(levelname)-8s | %(name)s | %(message)s",
 24    format="%(levelname)-8s | %(name)s | %(message)s",
 25)
 26LOGGER = logging.getLogger(__name__)
 27# LOG_FORMAT = "%(asctime)s | %(levelname)-8s | %(name)s | %(message)s"
 28
 29from cromosim.domain import Destination, Domain
 30from cromosim.micro import (
 31    people_initialization,
 32    add_people_flows,
 33    append_people,
 34    compute_contacts,
 35    find_duplicate_people,
 36    move_people,
 37    people_update_destination,
 38    projection,
 39)
 40from cromosim.visu import (
 41    plot_sensors
 42)
 43
 44# plt.ion()
 45
 46"""
 47    python3 micro_granular.py --json input.json
 48"""
 49
 50parser = OptionParser(usage="usage: %prog [options] filename",
 51                      version="%prog 1.0")
 52parser.add_option('--json', dest="jsonfilename", default="input.json",
 53                  type="string", action="store", help="Input json filename")
 54opt, remainder = parser.parse_args()
 55
 56LOGGER.info(f"Loading JSON configuration: {opt.jsonfilename}")
 57
 58with open(opt.jsonfilename) as json_file:
 59    try:
 60        config = json.load(json_file)
 61    except json.JSONDecodeError as exc:
 62        raise ValueError(f"Invalid JSON in {opt.jsonfilename}: {exc}") from exc
 63
 64LOGGER.info(
 65    "Simulation: Tf=%s, dt=%s, projection_method=%s",
 66    config["Tf"],
 67    config["dt"],
 68    config["projection_method"],
 69)
 70LOGGER.info(
 71    "Contacts: dmax=%s, dmin_people=%s, dmin_walls=%s",
 72    config["dmax"],
 73    config["dmin_people"],
 74    config["dmin_walls"],
 75)
 76LOGGER.info(
 77    "Visualization: enabled=%s, backend=%s, people=%s " 
 78    "contacts=%s velocities=%s desired_velocities=%s "
 79    "paths=%s sensors=%s plot_period=%s",
 80    config["with_graphes"],
 81    config["plot_module"],
 82    config["plot_people"],
 83    config["plot_contacts"],
 84    config["plot_velocities"],
 85    config["plot_desired_velocities"],
 86    config["plot_paths"],
 87    config["plot_sensors"],
 88    config["plot_period"],
 89)
 90LOGGER.info(
 91    "Miscellaneous: prefix=%s, seed=%s, dmin_walls=%s",
 92    config["prefix"],
 93    config["seed"],
 94    config["dmin_walls"],
 95)
 96
 97os.makedirs(config["prefix"], exist_ok=True)
 98
 99"""
100    Build the Domain objects
101"""
102
103domains = {}
104
105domain_configs = config["domains"]
106
107for domain_index, domain_config in enumerate(domain_configs):
108    name = domain_config["name"]
109    LOGGER.info(f"Start building domain {domain_index}: {name}")
110    if domain_config["background"] == "":
111        domain = Domain(
112            name=name,
113            pixel_size=domain_config["px"],
114            width=domain_config["width"],
115            height=domain_config["height"],
116            wall_colors=domain_config["wall_colors"],
117        )
118    else:
119        domain = Domain(
120            name=name,
121            background=domain_config["background"],
122            pixel_size=domain_config["px"],
123            wall_colors=domain_config["wall_colors"],
124        )
125        
126    LOGGER.info(
127        f"Building domain {domain_index}: {name} => "
128        f"Add geometric objects described in the JSON "
129        f"configuration")
130    
131    for shape_config in domain_config["shape_lines"]:
132        shape = Line2D(
133            shape_config["xx"], shape_config["yy"],
134            linewidth=shape_config["linewidth"],
135        )
136        domain.add_shape(
137            shape,
138            outline_color=shape_config["outline_color"],
139            fill_color=shape_config["fill_color"],
140        )
141
142    for shape_config in domain_config["shape_circles"]:
143        shape = Circle(
144            (shape_config["center_x"], shape_config["center_y"]),
145            shape_config["radius"],
146        )
147        domain.add_shape(
148            shape,
149            outline_color=shape_config["outline_color"],
150            fill_color=shape_config["fill_color"],
151        )
152
153    for shape_config in domain_config["shape_ellipses"]:
154        shape = Ellipse(
155            (shape_config["center_x"], shape_config["center_y"]),
156            shape_config["width"], shape_config["height"],
157            shape_config["angle_in_degrees_anti-clockwise"],
158        )
159        domain.add_shape(
160            shape,
161            outline_color=shape_config["outline_color"],
162            fill_color=shape_config["fill_color"],
163        )
164
165    for shape_config in domain_config["shape_rectangles"]:
166        shape = Rectangle(
167            (shape_config["bottom_left_x"], shape_config["bottom_left_y"]),
168            shape_config["width"], shape_config["height"],
169            shape_config["angle_in_degrees_anti-clockwise"],
170        )
171        domain.add_shape(
172            shape,
173            outline_color=shape_config["outline_color"],
174            fill_color=shape_config["fill_color"],
175        )
176
177    for shape_config in domain_config["shape_polygons"]:
178        shape = Polygon(shape_config["xy"])
179        domain.add_shape(
180            shape,
181            outline_color=shape_config["outline_color"],
182            fill_color=shape_config["fill_color"],
183        )
184        
185    LOGGER.info(f"Building domain {domain_index}: {name} => build_domain()")    
186    domain.build_domain()
187    
188    LOGGER.info(
189        f"Building domain {domain_index}: {name} => "
190        f"Add the destinations described in the JSON "
191        f"configuration")
192
193    for destination_index, destination_config in enumerate(
194            domain_config["destinations"]
195        ):
196        desired_velocity_from_color = [
197            np.concatenate((cdv["color"], cdv["desired_velocity"]))
198            for cdv in destination_config["desired_velocity_from_color"]
199        ]
200        destination = Destination(
201            name=destination_config["name"], 
202            colors=destination_config["colors"],
203            excluded_colors=destination_config["excluded_colors"],
204            desired_velocity_from_color=desired_velocity_from_color,
205            velocity_scale=destination_config["velocity_scale"],
206            next_destination=destination_config["next_destination"],
207            next_domain=destination_config["next_domain"],
208            next_destination_weights=destination_config["next_destination_weights"],
209            next_transit_box=destination_config["next_transit_box"],
210            next_transit_box_xshift=destination_config["next_transit_box_xshift"],
211            next_transit_box_yshift=destination_config["next_transit_box_yshift"],
212        )
213        LOGGER.info(
214                f"Building domain {domain_index}: {name} => "
215                f"Add destination {destination_index}: {destination_config["name"]}")
216        domain.add_destination(destination)
217        if (config["with_graphes"]):
218            domain.plot_desired_velocity(
219                destination_config["name"], 
220                id=100*domain_index+10+destination_index, 
221                step=20,
222                title=f"Destination '{destination_config["name"]}' "
223                      f"in '{name}' domain")
224     
225    if (config["with_graphes"]):
226        domain.plot(
227            id=100*domain_index, 
228            title=f"Geometry of the '{name}' domain")
229        domain.plot_wall_dist(
230            id=100*domain_index+1, 
231            step=20,
232            title=f"Distance from walls for the '{name}' domain")
233       
234    domains[name] = domain
235
236
237"""
238    To create the sensors to measure the pedestrian flows
239"""
240
241all_sensors = {domain_name: [] for domain_name in domains}
242
243for sensor_config in config["sensors"]:
244    sensor = dict(sensor_config)
245    sensor["id"] = []
246    sensor["times"] = []
247    sensor["xy"] = []
248    sensor["dir"] = []
249    LOGGER.info(
250        f"Initializing the '{sensor["name"]}' sensor "
251        f"located in the '{sensor["domain"]}' domain")
252    all_sensors[sensor["domain"]].append(sensor)
253
254
255"""
256    Initialize pedestrian populations in all configured domains
257"""
258
259# Current time
260t = 0.0
261counter = 0
262dt = config["dt"]
263Tf = config["Tf"]
264
265all_people = {}
266
267for people_index, people_config in enumerate(config["people_init"]):
268    domain_name = people_config["domain"]
269    domain = domains[domain_name]
270    LOGGER.info(
271        f"Initializing people configuration in the '{domain_name}' domain"
272    )
273
274    people = people_initialization(
275        domain,
276        people_config["groups"],
277        dt,
278        dmin_people=config["dmin_people"],
279        dmin_walls=config["dmin_walls"],
280        seed=config["seed"],
281        itermax=10,
282        projection_method="uzawa",
283        flow_duration=60,
284        log=True,
285    )
286    if people is None:
287        raise RuntimeError(
288            f"No population was created for domain {domain_name!r}."
289        )
290    else:
291        all_people[domain_name] = people
292
293
294"""
295    Visualization: either using Matplotlib or pyvista 
296"""
297
298all_visu = {} # one per domain...
299
300for domain_index, (name, domain) in enumerate(domains.items()):
301
302    LOGGER.info(
303        f"Creating the view vindow for the '{domain_name}' domain "
304        f"(render={config["plot_module"]})"
305    )
306    
307    if not config["with_graphes"]:
308    
309        all_visu[name] = None
310    
311    else:
312
313        people = all_people[name]
314        colors = people["xyrv"][:, 2]
315        if colors.size>0:
316            cmin, cmax = colors.min(), colors.max()
317        else:
318            cmin, cmax = 0, 1
319                
320        II, JJ, Vd = domain.people_desired_velocity(
321            people["xyrv"], 
322            people["destinations"]
323        )
324        people["Vd"] = Vd
325                
326        contacts = compute_contacts(domain, people["xyrv"], config["dmax"])      
327
328        if config["plot_module"] == "matplotlib":
329            
330            from cromosim.visu import matplotlibVisualization
331            
332            visu = matplotlibVisualization(
333                domain,
334                figure_id=100 * (1+domain_index) + 20,
335                cmap="viridis",
336                velocity_reference=config["plot_velocity_reference"],
337                velocity_arrow_length=config["plot_velocity_arrow_length"]
338            )
339            visu.update(
340                people=people,
341                contacts=contacts,
342                sensors=all_sensors[name],
343                time=t,
344                plot_people=config["plot_people"],
345                plot_contacts=config["plot_contacts"],
346                plot_paths=config["plot_paths"],
347                plot_velocities=config["plot_velocities"],
348                plot_desired_velocities=config["plot_desired_velocities"],
349                plot_sensors=config["plot_sensors"],
350                color_scalar_values=colors,
351                color_scalar_name="r",
352                color_minmax=(cmin, cmax),
353                color_cmap="viridis",
354                show_scalar_bar=config["plot_people"],
355                savefig=True,
356                filename=os.path.join(
357                    config["prefix"],
358                    f"{name}_fig_{counter:06d}.png",
359                ),
360            )
361            plt.pause(0.001)
362            visu.show()
363            all_visu[name] = visu
364
365        elif config["plot_module"] == "pyvista":
366        
367            from cromosim.visu import pyvistaVisualization
368        
369            target_width = 1920
370            target_height = round(
371                target_width * domain.height / domain.width
372            )    
373            # scale = max(1, round(1920 / domain.width))
374            movie_path = os.path.join(config["prefix"], f"{name}.mp4")
375            visu = pyvistaVisualization(
376                domain,
377                image_path=domain._Domain__image_filename,
378                pixel_size=domain.pixel_size,
379                movie_path=movie_path,
380                fps=4, # max(1,int(round(1/(config["plot_period"]*dt)))),
381                # window_size=(scale * domain.width, scale * domain.height),
382                window_size=(target_width, target_height),
383                movie_macro_block_size=16,
384                velocity_reference=config["plot_velocity_reference"],
385                velocity_arrow_length=config["plot_velocity_arrow_length"]
386            )
387            visu.update(
388                people=people,
389                contacts=contacts,
390                sensors=all_sensors[name],
391                time=t,
392                plot_people=config["plot_people"],
393                plot_contacts=config["plot_contacts"],
394                plot_velocities=config["plot_velocities"],
395                plot_desired_velocities=config["plot_desired_velocities"],
396                plot_paths=config["plot_paths"],
397                plot_sensors=config["plot_sensors"],
398                color_scalar_values=colors,
399                color_scalar_name="r",
400                color_minmax=(cmin, cmax),
401                color_cmap="viridis",
402                show_scalar_bar=config["plot_people"],
403                write_frame=True,
404            )
405            visu.show()
406            all_visu[name] = visu
407                
408
409# It may be helpful to pause here to check the calculation 
410# domains and the positions of the individuals before 
411# continuing the simulation...
412# sys.exit()
413
414"""
415    Main loop
416"""
417
418cc = 0
419draw = True
420
421while (t < Tf):
422    
423    LOGGER.info( f"-------------------------------------\n" )
424    LOGGER.info( f"===> Time = {t}" )
425
426    # If there are inflows of people arriving in the area, we calculate 
427    # the number of people to be added.
428    
429    for people_index, people_config in enumerate(config["people_init"]):
430        domain_name = people_config["domain"]
431        domain = domains[domain_name]
432        people_to_add = add_people_flows(
433            domain,
434            people_config["groups"],
435            dt,
436            dmin_people=config["dmin_people"],
437            dmin_walls=config["dmin_walls"],
438            rng=all_people[domain_name]["rng"],
439            itermax=10,
440            projection_method="uzawa",
441            flow_duration=config["dt"],
442            log=True,
443        )
444        if people_to_add is not None: 
445            LOGGER.info(
446                f"People flow in the '{domain_name}' domain: "
447                f"number of people added = {people_to_add["xyrv"].shape[0]}"
448            )
449            append_people(all_people[domain_name], people_to_add)
450                            
451 
452     # Compute people desired velocity
453     
454    for domain_index, (name, domain) in enumerate(domains.items()):
455        LOGGER.info(
456            f"Compute people desired velocity in the '{name}' domain"
457            )
458        people = all_people[name]
459        II, JJ, Vd = domain.people_desired_velocity(
460            people["xyrv"],
461            people["destinations"])
462        people["Vd"][:,0] = people["beta"]*Vd[:,0]
463        people["Vd"][:,1] = people["beta"]*Vd[:,1]
464        people["I"] = II
465        people["J"] = JJ
466
467    # Look at if there are people in the transit boxes
468    
469    LOGGER.info(
470                f"Find the people who need to be duplicated when "
471                f"they are in an area where two domains overlap "
472                f"(such as when changing floors)"
473            )
474    virtual_people = find_duplicate_people(all_people, domains)
475 
476    # Projection
477    
478    for domain_index, (name, domain) in enumerate(domains.items()):
479        LOGGER.info(
480                    f"Start projection step for the '{name}' domain "
481                    f"to prevent individuals from overlapping..."
482                )
483        people = all_people[name]
484
485        if virtual_people[name]["xyrv"] is not None:
486            LOGGER.info(
487                        f"Projection step for the '{name}' domain: "
488                        f"{virtual_people[name]["xyrv"].shape[0]} "
489                        f"virtual guys found"
490                    )
491            xyrv = np.concatenate(
492                (people["xyrv"], virtual_people[name]["xyrv"]))
493            Vd = np.concatenate(
494                (people["Vd"], virtual_people[name]["Vd"]))
495            Mass = np.concatenate(
496                (people["mass"], virtual_people[name]["mass"]))
497        else:
498            xyrv = people["xyrv"]
499            Vd = people["Vd"]
500            Mass = people["mass"]
501
502        if (xyrv.shape[0] > 0):
503
504            if (np.unique(xyrv, axis=0).shape[0] != xyrv.shape[0]):
505                raise ValueError( 
506                                 "There are two identical lines in "
507                                 "the array xyrv used to determine the "
508                                 "contacts between individuals and this "
509                                 "is not normal!"
510                            )
511                         
512            contacts = compute_contacts(domain, xyrv, config["dmax"])
513            LOGGER.info(
514                f"Projection step for the '{name}' domain: "
515                f"{contacts.shape[0]} contacts computed"
516            )
517
518            start = time.perf_counter()
519            info, B, U, L, P, frustrations = projection(
520                dt, xyrv,
521                contacts, Vd, Mass,
522                dmin_people=config["dmin_people"],
523                dmin_walls=config["dmin_walls"],
524                method=config["projection_method"],
525                log=True)
526            cputime_projection = time.perf_counter() - start
527            LOGGER.info(
528                        f"Projection step for the '{name}' domain: "
529                        f"info={info} (iterations for the projection algorithm) "
530                        f"cputime={np.round(cputime_projection,4)}"
531                    )
532              
533            LOGGER.info(
534                    f"Update the positions of the individuals using "
535                    f"the new permissible velocities obtained"
536                )
537                                    
538            nn = people["xyrv"].shape[0]
539            all_people[name]["U"] = U[:nn, :]
540            
541            virtual_people[name]["U"] = U[nn:, :]
542               
543            # We retain the last frustration values and compute beta 
544                 
545            all_people[name]["frustration_last_pos"] = (all_people[name]["frustration_last_pos"] + 1) % \
546                config["frustration_memory_size"]
547            
548            all_people[name]["frustration_memory"      ][ np.arange(nn), all_people[name]["frustration_last_pos"] ] = frustrations[:nn]
549            all_people[name]["frustration_memory_times"][ np.arange(nn), all_people[name]["frustration_last_pos"] ] = t
550            
551            # The correcting factor, beta, is activated only 
552            # for individual who have neighbors in their cone
553            # of vision and at a distance less than a disc diameter.
554            
555            all_people[name]["beta"] = np.ones((nn,))
556            
557            active = np.zeros((nn,))
558            for ic, cont in enumerate(contacts):
559                i = int(cont[0])
560                j = int(cont[1])
561                dij = cont[2]
562                eij_x = cont[3]
563                eij_y = cont[4]
564                if (j > -1):  # contact person/person
565                    # theta_ij is the angle at which i (assumed to look 
566                    # ahead along their desired velocity) sees j
567                    # Angular dependence
568                    norm_Vdi = np.sqrt(Vd[i, 0]**2+Vd[i, 1]**2)
569                    if (norm_Vdi > 0) and (i<nn): # not virtual people...
570                        theta_ij = np.arccos( max(-1.0, min(1.0, ( (Vd[i,0]*eij_x+Vd[i,1]*eij_y)/norm_Vdi ))))
571                        if theta_ij<=all_people[name]["viewing_angle"][i]/180.0*np.pi:
572                            active[i] = 1
573                                                    
574                    norm_Vdj = np.sqrt(Vd[j, 0]**2+Vd[j, 1]**2)
575                    if (norm_Vdj > 0) and (j<nn): # not virtual people...
576                        theta_ji = np.arccos( max(-1.0, min(1.0, ( -(Vd[j,0]*eij_x+Vd[j,1]*eij_y)/norm_Vdj ))))
577                        if theta_ji<=all_people[name]["viewing_angle"][j]/180.0*np.pi:
578                            active[j] = 1
579
580            ind_frust = np.where( (all_people[name]["with_frustration"]==1) * \
581                (active==1) )[0]
582                                                                
583            all_people[name]["beta"][ind_frust] = 1.0 - \
584                np.mean(all_people[name]["frustration_memory"], axis=1)[ind_frust]
585                                    
586            # Update people positions
587            
588            all_people[name], all_sensors[name] = move_people(
589                t, dt,
590                all_people[name],
591                all_sensors[name])
592
593    # Update people destinations
594    
595    # the distance at which one changes destinations:
596    dist_to_change_dest = 2*dt
597    all_people = people_update_destination(
598        all_people, domains, 
599        dist_to_change_dest)
600
601  
602    t += dt
603    cc += 1
604    counter += 1
605  
606    for domain_index, (name, domain) in enumerate(domains.items()):
607        LOGGER.info(
608            f"--- Number of persons in the '{name}' "
609            f"domain: {all_people[name]["xyrv"].shape[0]}"
610        )    
611
612    if (cc >= config["plot_period"]):
613        draw = True
614        cc = 0
615    else:
616        draw = False
617    
618    # Ntot = sum([people["xyrv"].shape[0] for (name,people) in all_people.items()])
619    # if (Ntot==0):
620    #     LOGGER.info(
621    #         f"No one left in any domain! END..."
622    #     )
623    #     break    
624    
625    """
626    Visualization: either using Matplotlib or pyvista 
627    """
628
629    if draw and config["with_graphes"]:
630
631        for domain_index, (name, domain) in enumerate(domains.items()):
632
633            people = all_people[name]
634            
635            colors = people["xyrv"][:, 2]
636            if colors.size>0:
637                cmin, cmax = colors.min(), colors.max()
638            else:
639                cmin, cmax = 0, 1
640            
641            # Since people have moved, we recalculate Vd and 
642            # the contacts to ensure a correct display
643            II, JJ, Vd = domain.people_desired_velocity(
644                people["xyrv"], 
645                people["destinations"]
646            )
647            people["Vd"][:,0] = people["beta"]*Vd[:,0]
648            people["Vd"][:,1] = people["beta"]*Vd[:,1]        
649            
650            contacts = compute_contacts(domain, people["xyrv"], config["dmax"]) 
651            
652            if config["plot_module"] == "matplotlib":
653            
654                all_visu[name].update(
655                    people=people,
656                    contacts=contacts,
657                    virtual_people=virtual_people[name],
658                    sensors=all_sensors[name],
659                    time=t,
660                    plot_people=config["plot_people"],
661                    plot_contacts=config["plot_contacts"],
662                    plot_paths=config["plot_paths"],
663                    plot_velocities=config["plot_velocities"],
664                    plot_desired_velocities=config["plot_desired_velocities"],
665                    plot_sensors=config["plot_sensors"],
666                    color_scalar_values=colors,
667                    color_scalar_name="r",
668                    color_minmax=(cmin, cmax),
669                    color_cmap="viridis",
670                    show_scalar_bar=config["plot_people"],
671                    savefig=True,
672                    filename=os.path.join(
673                        config["prefix"],
674                        f"{name}_fig_{counter:06d}.png",
675                    ),
676                )
677                plt.pause(0.001)
678                all_visu[name].show()
679
680            elif config["plot_module"] == "pyvista":
681                       
682                all_visu[name].update(
683                    people=people,
684                    contacts=contacts,
685                    virtual_people=virtual_people[name],
686                    sensors=all_sensors[name],
687                    time=t,
688                    plot_people=config["plot_people"],
689                    plot_contacts=config["plot_contacts"],
690                    plot_velocities=config["plot_velocities"],
691                    plot_desired_velocities=config["plot_desired_velocities"],
692                    plot_paths=config["plot_paths"],
693                    plot_sensors=config["plot_sensors"],
694                    color_scalar_values=colors,
695                    color_scalar_name="r",
696                    color_minmax=(cmin, cmax),
697                    color_cmap="viridis",
698                    show_scalar_bar=config["plot_people"],
699                    write_frame=True,
700                )
701                all_visu[name].process_events()
702                all_visu[name].plotter.update()
703                all_visu[name].show()      
704
705
706# End...
707
708LOGGER.info(
709    f"============ It's over! ============\n"
710)    
711
712
713if config["with_graphes"]:
714    
715    for domain_index, (name, domain) in enumerate(domains.items()):
716        
717        if config["plot_module"] == "pyvista":
718            LOGGER.info(
719                f"Close the pyvista visualization window for the '{name}' domain"
720            )    
721            all_visu[name].close()
722        
723        if config["plot_sensors"]:
724            sensors = all_sensors[domain_name]
725            LOGGER.info(
726                f"Plot sensor data for the '{name}' domain "
727                f"({len(sensors)} sensors)"
728            )    
729            plot_sensors(
730                100*domain_index+40, 
731                sensors, 
732                t, 
733                savefig=True,
734                filename=os.path.join(
735                    config["prefix"],
736                    f"sensors_{domain_name}.png"
737                    )
738                )
739            plt.pause(0.001)
740
741plt.show()