diff --git a/.env b/.env index 1c09a84e..320b2d34 100644 --- a/.env +++ b/.env @@ -1,101 +1,67 @@ -# Basic simulation parameters - -# Number of agents and resc. patches -N = 3 -N_RESOURCES = 10 - -# Simulation time -T = 1000 - -# Framerate (only matter if visualization is turned on) -INIT_FRAMERATE = 25 - -# Visualization ON/OFF -WITH_VISUALIZATION = 1 - -# Resolution (px) of agents' visual projection fields -VISUAL_FIELD_RESOLUTION = 1200 - -# Enviornment size (width and height in px) -ENV_WIDTH = 600 -ENV_HEIGHT = 600 - -# Agent and Resource sizes -RADIUS_AGENT = 10 -RADIUS_RESOURCE = 45 -MIN_RESOURCE_PER_PATCH = 400 -MAX_RESOURCE_PER_PATCH = 450 -REGENERATE_PATCHES = 1 -PATCH_BORDER_OVERLAP = 1 - -# Exploitation -# max consumption of an agent per timestep -AGENT_CONSUMPTION = 1 -# resource patch quality -MIN_RESOURCE_QUALITY = 0.05 -MAX_RESOURCE_QUALITY = 0.3 -TELEPORT_TO_MIDDLE = 0 -GHOST_WHILE_EXPLOIT = 1 -AGENT_AGENT_COLLISION = 1 -PATCHWISE_SOCIAL_EXCLUSION =1 - -# Vision -# fov is symmetric and defined with percent of pi. e.g -# if 0.6 then fov is (-0.6*pi, 0.6*pi). 1 is full 360 degree vision -AGENT_FOV = 0.45 -VISION_RANGE = 200 -VISUAL_EXCLUSION = 1 - -# Visualize agents' visual field on screen by default or when return is pressed -SHOW_VISUAL_FIELDS = 0 -SHOW_VISUAL_FIELDS_RETURN = 1 -SHOW_VISION_RANGE = 1 - -# Time needed to acquire info from env and initial probability for pooling process -POOLING_TIME = 0 -POOLING_PROBABILITY = 0.05 - -# Monitoring Related parameters (True only tested on Ubuntu machines with -# initialized grafana and ifdb instances as in readme) -# log data in influxdb -USE_IFDB_LOGGING = 0 -# logs data in RAM -USE_RAM_LOGGING = 1 -# using zarr compression -USE_ZARR_FORMAT = 1 -# saves csv files (if ifdb logging is used, json files otherwise) -SAVE_CSV_FILES = 0 - -# Decision process parameters -DEC_TW = 0.5 -DEC_EPSW = 3 -DEC_GW = 0.085 -DEC_BW = 0 -DEC_WMAX = 1 - -DEC_TU = 0.5 -DEC_EPSU = 3 -DEC_GU = 0.085 -DEC_BU = 0 -DEC_UMAX = 1 - -DEC_SWU = 0.25 -DEC_SUW = 0.01 - -DEC_TAU = 10 -DEC_FN = 2 -DEC_FR = 1 - -# Movement parameters -# Exploration -MOV_EXP_VEL_MIN = 1 -MOV_EXP_VEL_MAX = 1 -MOV_EXP_TH_MIN = -0.3 -MOV_EXP_TH_MAX = 0.3 - -# Relocation -MOV_REL_DES_VEL = 1 -MOV_REL_TH_MAX = 0.5 - -# Exploitation -CONS_STOP_RATIO = 0.08 \ No newline at end of file +N=10 +N_RESOURCES=1 +T=5000 +INIT_FRAMERATE=25 +WITH_VISUALIZATION=0 +VISUAL_FIELD_RESOLUTION=320 +ENV_WIDTH=900 +ENV_HEIGHT=900 +RADIUS_AGENT=10 +RADIUS_RESOURCE=20 +MIN_RESOURCE_PER_PATCH=100 +MAX_RESOURCE_PER_PATCH=-1 +REGENERATE_PATCHES=1 +PATCH_BORDER_OVERLAP=1 +MAX_SPEED=1.5 +AGENT_CONSUMPTION=1 +MIN_RESOURCE_QUALITY=0.25 +MAX_RESOURCE_QUALITY=-1 +TELEPORT_TO_MIDDLE=0 +GHOST_WHILE_EXPLOIT=1 +AGENT_AGENT_COLLISION=0 +PATCHWISE_SOCIAL_EXCLUSION=1 +AGENT_FOV=0.5 +VISION_RANGE=2000 +VISUAL_EXCLUSION=0 +SHOW_VISUAL_FIELDS=0 +SHOW_VISUAL_FIELDS_RETURN=0 +SHOW_VISION_RANGE=0 +POOLING_TIME=0 +POOLING_PROBABILITY=0.05 +USE_IFDB_LOGGING=1 +USE_RAM_LOGGING=1 +USE_ZARR_FORMAT=1 +SAVE_CSV_FILES=1 +DEC_TW=0.5 +DEC_EPSW=0 +DEC_GW=0.085 +DEC_BW=0 +DEC_WMAX=1 +DEC_TU=0.5 +DEC_EPSU=1 +DEC_GU=0.085 +DEC_BU=0 +DEC_UMAX=1 +DEC_SWU=0 +DEC_SUW=0 +DEC_TAU=10 +DEC_FN=0.5 +DEC_FR=0.5 +MOV_EXP_VEL_MIN=3 +MOV_EXP_VEL_MAX=3 +MOV_EXP_TH_MIN=-0.5 +MOV_EXP_TH_MAX=0.5 +MOV_REL_DES_VEL=3 +MOV_REL_TH_MAX=1.8 +MEMORY_DEPTH=1 +CONS_STOP_RATIO=0.175 +APP_VERSION=VisualFlocking +VF_GAM=0.1 +VF_V0=1 +VF_ALP0=0 +VF_ALP1=0.0014 +VF_ALP2=0 +VF_BET0=0 +VF_BET1=0.0014 +VF_BET2=0 +SAVE_ROOT_DIR=abm/data/simulation_data/VFinfinite/batch_0 diff --git a/abm/agent/supcalc.py b/abm/agent/supcalc.py index c20a6214..ff82adb5 100644 --- a/abm/agent/supcalc.py +++ b/abm/agent/supcalc.py @@ -2,89 +2,7 @@ calc.py : Supplementary methods and calculations necessary for agents """ import numpy as np -from scipy import integrate -from abm.contrib import vswrm, movement_params - - -# Functions needed for VSWRM functionality -def VSWRM_flocking_state_variables(vel_now, Phi, V_now, t_now=None, V_prev=None, t_prev=None): - """Calculating state variables of a given agent according to the main algorithm as in - https://advances.sciencemag.org/content/6/6/eaay0792. - Args: - vel_now: current speed of the agent - V_now: current binary visual projection field array - Phi: linspace numpy array of visual field axis - t_now: current time - V_prev: previous binary visual projection field array - t_prev: previous time - Returns: - dvel: temporal change in agent velocity - dpsi: temporal change in agent heading angle - - """ - # # Deriving over t - # if V_prev is not None and t_prev is not None and t_now is not None: - # dt = t_now - t_prev - # logger.debug('Movement calculation called with NONE as time-related parameters.') - # joined_V = np.vstack((V_prev, t_prev)) - # dt_V = dt_V_of(dt, joined_V) - # else: - # dt_V = np.zeros(len(Phi)) - - dt_V = np.zeros(len(Phi)) - - # Deriving over Phi - dPhi_V = dPhi_V_of(Phi, V_now) - - # Calculating series expansion of functional G - G_vel = (-V_now + vswrm.ALP2 * dt_V) - - # Spikey parts shall be handled separately because of numerical integration - G_vel_spike = np.square(dPhi_V) - - G_psi = (-V_now + vswrm.BET2 * dt_V) - - # Spikey parts shall be handled separately because of numerical integration - G_psi_spike = np.square(dPhi_V) - - # Calculating change in velocity and heading direction - dPhi = Phi[-1] - Phi[-2] - FOV_rescaling_cos = 1 - FOV_rescaling_sin = 1 - - # print(f"alp0 : {vswrm.ALP0 * integrate.trapz(np.cos(FOV_rescaling_cos * Phi) * G_vel, Phi)}", ) - # print(f'alp1 : {vswrm.ALP0 * vswrm.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) * dPhi}') - - dvel = vswrm.GAM * (vswrm.V0 - vel_now) + \ - vswrm.ALP0 * integrate.trapz(np.cos(FOV_rescaling_cos * Phi) * G_vel, Phi) + \ - vswrm.ALP0 * vswrm.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) * dPhi - dpsi = vswrm.BET0 * integrate.trapz(np.sin(Phi) * G_psi, Phi) + \ - vswrm.BET0 * vswrm.BET1 * np.sum(np.sin(FOV_rescaling_sin * Phi) * G_psi_spike) * dPhi - - return dvel, dpsi - - -def dPhi_V_of(Phi, V): - """Calculating derivative of VPF according to Phi visual angle array at a given timepoint t - Args: - Phi: linspace numpy array of visual field axis - V: binary visual projection field array - Returns: - dPhi_V: derivative array of V w.r.t Phi - """ - # circular padding for edge cases - padV = np.pad(V, (1, 1), 'wrap') - dPhi_V_raw = np.diff(padV) - - # we want to include non-zero value if it is on the edge - if dPhi_V_raw[0] > 0 and dPhi_V_raw[-1] > 0: - dPhi_V_raw = dPhi_V_raw[0:-1] - - else: - dPhi_V_raw = dPhi_V_raw[1:, ...] - - dPhi_V = dPhi_V_raw / (Phi[-1] - Phi[-2]) - return dPhi_V +from abm.contrib import movement_params def find_nearest(array, value): @@ -128,6 +46,15 @@ def random_walk(desired_vel=None): movement_params.exp_theta_max) return dvel, dtheta +def distance_torus(p0, p1, dimensions): + """Calculating distance between 2 2D points p0 and p1 as nparrays in an arena + with dimensions as in dimensions with infinite boundary conditions. + po and p1 can be a set of 2D coordinates e.g. np.array([[x0, y0],[x1, y1],[x2, y2]]) + """ + delta = np.abs(p0 - p1) + delta = np.where(delta > 0.5 * dimensions, delta - dimensions, delta) + return np.sqrt((delta ** 2).sum(axis=1)) + def distance_coords(x1, y1, x2, y2, vectorized=False): """Distance between 2 points in 2D space calculated from point coordinates. diff --git a/abm/app_visual_flocking.py b/abm/app_visual_flocking.py new file mode 100644 index 00000000..ed16dd12 --- /dev/null +++ b/abm/app_visual_flocking.py @@ -0,0 +1,108 @@ +import os +import shutil +from pathlib import Path + +from contextlib import ExitStack +from dotenv import dotenv_values + +from abm.app import save_isims_env +from abm.projects.visual_flocking.vf_simulation.vf_isims import VFPlaygroundSimulation +from abm.projects.visual_flocking.vf_contrib.vf_playgroundtool import setup_visflock_playground +from abm.projects.visual_flocking.vf_simulation.vf_sims import VFSimulation + + +def setup_environment(): + EXP_NAME = os.getenv("EXPERIMENT_NAME", "visual_flocking") + EXP_NAME_COPY = f"{EXP_NAME}_copy" + os.path.dirname(os.path.realpath(__file__)) + root_abm_dir = Path(__file__).parent.parent + env_file_dir = root_abm_dir / "abm" / "projects" / "visual_flocking" # Path(__file__).parent + env_path = env_file_dir / f"{EXP_NAME}.env" + env_path_copy = env_file_dir / f"{EXP_NAME_COPY}.env" + # make a duplicate of the env file to be used by the playground + shutil.copyfile(env_path, env_path_copy) + envconf = dotenv_values(env_path) + return env_file_dir, EXP_NAME_COPY, envconf + + +def start_playground(): + """starting simulation with interactive interface""" + env_file_dir, EXP_NAME_COPY, envconf = setup_environment() + # changing env file according to playground default parameters before + # running any component of the SW + pgt = setup_visflock_playground() + save_isims_env(env_file_dir, EXP_NAME_COPY, pgt, envconf) + # Start interactive simulation + sim = VFPlaygroundSimulation() + sim.start() + + +def start(parallel=False, headless=False, agent_behave_param_list=None): + """starting simulation without interactive interface""" + # Define root abm directory from which env file is read out + root_abm_dir = os.path.dirname(os.path.dirname(os.path.realpath(__file__))) + + # Finding env file + EXP_NAME = os.getenv("EXPERIMENT_NAME", "") + env_path = os.path.join(root_abm_dir, f"{EXP_NAME}.env") + if os.path.isfile(env_path): + print(f"Read env vars from {env_path}") + envconf = dotenv_values(env_path) + app_version = envconf.get("APP_VERSION", "Base") + if app_version != "VisualFlocking": + raise Exception(".env file was not created for project visual " + "flocking or no APP_VERSION parameter found!") + else: + raise Exception(f"Could not find .env file under path {env_path}") + + # Initializing virtual display + window_pad = 30 + vscreen_width = int(float(envconf["ENV_WIDTH"])) + 2 * window_pad + 10 + vscreen_height = int(float(envconf["ENV_HEIGHT"])) + 2 * window_pad + 10 + + # Running headless simulation on virtual display + if headless: + # required to start pygame in headless mode + os.environ['SDL_VIDEODRIVER'] = 'dummy' + from xvfbwrapper import Xvfb + + with ExitStack() if not headless else Xvfb(width=vscreen_width, height=vscreen_height) as xvfb: + sim = VFSimulation(N=int(float(envconf["N"])), + T=int(float(envconf["T"])), + v_field_res=int(envconf["VISUAL_FIELD_RESOLUTION"]), + agent_fov=float(envconf['AGENT_FOV']), + framerate=int(float(envconf["INIT_FRAMERATE"])), + with_visualization=bool(int(float(envconf["WITH_VISUALIZATION"]))), + width=int(float(envconf["ENV_WIDTH"])), + height=int(float(envconf["ENV_HEIGHT"])), + show_vis_field=bool(int(float(envconf["SHOW_VISUAL_FIELDS"]))), + show_vis_field_return=bool(int(envconf['SHOW_VISUAL_FIELDS_RETURN'])), + pooling_time=int(float(envconf["POOLING_TIME"])), + pooling_prob=float(envconf["POOLING_PROBABILITY"]), + agent_radius=int(float(envconf["RADIUS_AGENT"])), + N_resc=int(float(envconf["N_RESOURCES"])), + allow_border_patch_overlap=bool(int(float(envconf["PATCH_BORDER_OVERLAP"]))), + min_resc_perpatch=int(float(envconf["MIN_RESOURCE_PER_PATCH"])), + max_resc_perpatch=int(float(envconf["MAX_RESOURCE_PER_PATCH"])), + min_resc_quality=float(envconf["MIN_RESOURCE_QUALITY"]), + max_resc_quality=float(envconf["MAX_RESOURCE_QUALITY"]), + patch_radius=int(float(envconf["RADIUS_RESOURCE"])), + regenerate_patches=bool(int(float(envconf["REGENERATE_PATCHES"]))), + agent_consumption=int(float(envconf["AGENT_CONSUMPTION"])), + ghost_mode=bool(int(float(envconf["GHOST_WHILE_EXPLOIT"]))), + patchwise_exclusion=bool(int(float(envconf["PATCHWISE_SOCIAL_EXCLUSION"]))), + teleport_exploit=bool(int(float(envconf["TELEPORT_TO_MIDDLE"]))), + vision_range=int(float(envconf["VISION_RANGE"])), + visual_exclusion=bool(int(float(envconf["VISUAL_EXCLUSION"]))), + show_vision_range=bool(int(float(envconf["SHOW_VISION_RANGE"]))), + use_ifdb_logging=bool(int(float(envconf["USE_IFDB_LOGGING"]))), + use_ram_logging=bool(int(float(envconf["USE_RAM_LOGGING"]))), + save_csv_files=bool(int(float(envconf["SAVE_CSV_FILES"]))), + use_zarr=bool(int(float(envconf["USE_ZARR_FORMAT"]))), + parallel=parallel, + window_pad=window_pad, + agent_behave_param_list=agent_behave_param_list, + collide_agents=bool(int(float(envconf["AGENT_AGENT_COLLISION"]))) + ) + sim.write_batch_size = 100 + sim.start() diff --git a/abm/contrib/ifdb_params.py b/abm/contrib/ifdb_params.py index 928de886..18833a61 100644 --- a/abm/contrib/ifdb_params.py +++ b/abm/contrib/ifdb_params.py @@ -22,7 +22,7 @@ if WRITE_EACH_POINT is not None: write_batch_size = 1 else: - T = float(int(envconf["T"])) + T = float(int(envconf.get("T", 1000))) if T <= 1000: write_batch_size = T else: diff --git a/abm/contrib/vswrm.py b/abm/contrib/vswrm.py deleted file mode 100644 index c38aed66..00000000 --- a/abm/contrib/vswrm.py +++ /dev/null @@ -1,11 +0,0 @@ -"""parameters for visdualswarm algorithm""" - -GAM = 0.1 -V0 = 1 -# todo: why does it work fine with negative sign instead of positive one? -ALP0 = 0 -ALP1 = 0.09 -ALP2 = 0 -BET0 = 0.01 -BET1 = 0.1 -BET2 = 0 \ No newline at end of file diff --git a/abm/data/metaprotocol/experiments/VFExp4c.py b/abm/data/metaprotocol/experiments/VFExp4c.py new file mode 100644 index 00000000..905bcb49 --- /dev/null +++ b/abm/data/metaprotocol/experiments/VFExp4c.py @@ -0,0 +1,114 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 03.04.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastien and Romanczuk paper, + and checking the effect of field of view + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 900 +arena_h = 900 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 600), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.09), # 1/(2BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.09), # 1/(2BL) + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VFExp4b.py b/abm/data/metaprotocol/experiments/archive/VFExp4b.py new file mode 100644 index 00000000..d081fb94 --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VFExp4b.py @@ -0,0 +1,114 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 22.03.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastein and Romanczuk paper, + and checking the effect of field of view + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 500 +arena_h = 500 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 600), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.09), # 1/(2BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.09), # 1/(2BL) + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VFExp5b.py b/abm/data/metaprotocol/experiments/archive/VFExp5b.py new file mode 100644 index 00000000..aa9acd68 --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VFExp5b.py @@ -0,0 +1,114 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 22.03.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastein and Romanczuk paper, + and checking the effect of field of view with walls + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 500 +arena_h = 500 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 600), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.09), # 1/(2BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.09), # 1/(2BL) + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "walls"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp1.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp1.py new file mode 100644 index 00000000..2bc5173d --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp1.py @@ -0,0 +1,113 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: DD.MM.YYYY +Parameters: + Testing new subversion of software stack with collective signalling on computainbg cluster + +Project Maintainers (CoopSignalling): mezdahun & vagechrikov +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 900 +arena_h = 900 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("VISUAL_FIELD_RESOLUTION", 320), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.0014), + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.0014), + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.25, 0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp2.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp2.py new file mode 100644 index 00000000..8c1e390c --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp2.py @@ -0,0 +1,113 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 27.02.2023 +Parameters: + Comparing different conditions regarding FOV, alpha0 vs. beta0 with Wall reflections + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 900 +arena_h = 900 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("VISUAL_FIELD_RESOLUTION", 320), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.0014), + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.0014), + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.25, 0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "walls"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp3.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp3.py new file mode 100644 index 00000000..5ca5be57 --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp3.py @@ -0,0 +1,113 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 10.03.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastein and Romanczuk paper + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 1000 +arena_h = 1000 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 2000), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.036), # 1/(2.5BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.036), # 1/(2.5BL) + Constant("VF_BET2", 0), + Constant("AGENT_FOV", 1), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp4.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp4.py new file mode 100644 index 00000000..2d51248a --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp4.py @@ -0,0 +1,114 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 12.03.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastein and Romanczuk paper, + and checking the effect of field of view + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 1000 +arena_h = 1000 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 320), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.036), # 1/(2.5BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.036), # 1/(2.5BL) + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.25, 0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp5.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp5.py new file mode 100644 index 00000000..e5f62e7e --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/VFExp5.py @@ -0,0 +1,114 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: 12.03.2023 +Parameters: + Validating diff calculation error by reproducing main results from Bastein and Romanczuk paper, + and checking the effect of field of view with walls + +Project Maintainers (VisualFlockig): mezdahun +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 1000 +arena_h = 1000 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 0), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + # Constant("INIT_FRAMERATE", 100), + Constant("VISUAL_FIELD_RESOLUTION", 320), + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.036), # 1/(2.5BL) + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.036), # 1/(2.5BL) + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.25, 0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "walls"), + Constant("T", 25000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=True) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/archive/VisualFlocking/template/VisualFlockingTemplate.py b/abm/data/metaprotocol/experiments/archive/VisualFlocking/template/VisualFlockingTemplate.py new file mode 100644 index 00000000..8f18e58d --- /dev/null +++ b/abm/data/metaprotocol/experiments/archive/VisualFlocking/template/VisualFlockingTemplate.py @@ -0,0 +1,113 @@ +from abm.metarunner.metarunner import Tunable, Constant, MetaProtocol, TunedPairRestrain +import numpy as np +import os + +# On HPC the filename will be used as experiment name +# otherwise we will need to pass it for the script, such as +# EXPERIMENT_NAME=template python ExperimentCoopSigTemplate.py +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") +if EXP_NAME == "": + raise Exception("No experiment name has been passed") + +description_text = f""" +Experiment file using the MetaRunner interfacing language to define a set of criteria for batch simulations. +The filename shall not include underscore characters when running on HPC. This description will be included +as a txt file in the generated experiment data folder. + +Title: Experiment : {EXP_NAME} +Date: DD.MM.YYYY +Parameters: + Testing new subversion of software stack with collective signalling on computainbg cluster + +Project Maintainers (CoopSignalling): mezdahun & vagechrikov +""" + +# Defining fixed criteria for all automized simulations/experiments +arena_w = 900 +arena_h = 900 +fixed_criteria = [ + Constant("ENV_WIDTH", arena_w), + Constant("ENV_HEIGHT", arena_h), + Constant("USE_IFDB_LOGGING", 1), + Constant("USE_RAM_LOGGING", 1), # as we have plenty of resources we don't have to deal with IFDB on HPC + Constant("USE_ZARR_FORMAT", 1), + Constant("SAVE_CSV_FILES", 1), + Constant("WITH_VISUALIZATION", 1), + Constant("SHOW_VISUAL_FIELDS", 0), + Constant("SHOW_VISUAL_FIELDS_RETURN", 0), + Constant("SHOW_VISION_RANGE", 0), + Constant("TELEPORT_TO_MIDDLE", 0), + Constant("PATCHWISE_SOCIAL_EXCLUSION", 1), + Constant("POOLING_TIME", 0), + Constant("VISUAL_FIELD_RESOLUTION", 320), + Constant("AGENT_CONSUMPTION", 1), + Constant("MAX_RESOURCE_QUALITY", -1), # so that the minimum value will be used as definite + Constant("MAX_RESOURCE_PER_PATCH", -1), # so that the minimum value will be used as definite + Constant("MOV_EXP_VEL_MIN", 3), + Constant("MOV_EXP_VEL_MAX", 3), + Constant("MOV_REL_DES_VEL", 3), + Constant("MOV_EXP_TH_MIN", -0.5), + Constant("MOV_EXP_TH_MAX", 0.5), + Constant("MOV_REL_TH_MAX", 1.8), + Constant("CONS_STOP_RATIO", 0.175), + Constant("REGENERATE_PATCHES", 1), + Constant("DEC_FN", 0.5), + Constant("DEC_FR", 0.5), + Constant("DEC_TAU", 10), + Constant("DEC_BW", 0), + Constant("DEC_WMAX", 1), + Constant("DEC_BU", 0), + Constant("DEC_UMAX", 1), + Constant("DEC_GW", 0.085), + Constant("DEC_GU", 0.085), + Constant("DEC_TW", 0.5), + Constant("DEC_TU", 0.5), + Constant("PATCH_BORDER_OVERLAP", 1), + Constant("DEC_EPSW", 0), + Constant("DEC_EPSU", 1), + Constant("AGENT_AGENT_COLLISION", 0), + Constant("MIN_RESOURCE_PER_PATCH", 100), + Constant("N_RESOURCES", 1), + Constant("VISUAL_EXCLUSION", 0), # no visual occlusion + Constant("VISION_RANGE", 2000), # unlimited visual range + Constant("GHOST_WHILE_EXPLOIT", 1), + Constant("MIN_RESOURCE_QUALITY", 0.25), + Constant("RADIUS_RESOURCE", 20), + Constant("DEC_SWU", 0), # no cross-inhibition + Constant("DEC_SUW", 0), # no cross-inhibition +] + +# Defining decision param +arena_size = arena_w * arena_h +criteria_exp = [ + Constant("APP_VERSION", "VisualFlocking"), # Enabling project specific run in headless mode + Constant("INIT_FRAMERATE", 100), + Constant("RADIUS_AGENT", 5.5), + Constant("VF_GAM", 0.1), + Constant("VF_V0", 1), + Constant("VF_ALP1", 0.0014), + Constant("VF_ALP2", 0), + Constant("VF_BET1", 0.0014), + Constant("VF_BET2", 0), + Tunable("AGENT_FOV", values_override=[0.25, 0.5, 0.75, 1]), # unlimited FOV + Constant("BOUNDARY", "infinite"), + Constant("T", 2000), + Constant("N", 10), # Can not be tuned, must be fixed for Replay tool + Tunable("VF_ALP0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]), + Tunable("VF_BET0", values_override=[0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 3]) +] + +# Creating metaprotocol and add defined criteria +mp = MetaProtocol(experiment_name=EXP_NAME, num_batches=1, parallel=True, + description=description_text, headless=False) +for crit in fixed_criteria: + mp.add_criterion(crit) +for crit in criteria_exp: + mp.add_criterion(crit) + +# Generating temporary env files with criterion combinations. Comment this out if you want to continue simulating due +# to interruption +mp.generate_temp_env_files() + +# Running the simulations with project specific application +mp.run_protocols(project="VisualFlocking") diff --git a/abm/data/metaprotocol/experiments/scripts/VisualFlocking/figExp0prelim_create.py b/abm/data/metaprotocol/experiments/scripts/VisualFlocking/figExp0prelim_create.py new file mode 100644 index 00000000..5a7b0bfc --- /dev/null +++ b/abm/data/metaprotocol/experiments/scripts/VisualFlocking/figExp0prelim_create.py @@ -0,0 +1,178 @@ +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import numpy as np +import os +import json +import glob + +# under this path the individual summary statistics are saved +# with _N post tag. +exp_name = "VFprelimsummaryBugfix" +data_path = f"/media/david/DMezeySCIoI/ABMData/VisFlock/{exp_name}" + +titles = ["Polarization", "Inter-individual Distance", "Agent Collisions"] +file_names = ["polarization", "meaniid", "aacoll"] +conditions = ["infinite", "walls"] + +FOVs = [0.25, 0.5, 0.75, 1] +alphas = [0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5] #[0.0, 0.01, 0.1, 0.25, 0.5, 0.75, 1.0, 1.5, 3.0] +betas = [0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.75, 1, 2, 5] #[0.0, 0.01, 0.1, 0.25, 0.5, 0.75, 1.0, 1.5, 3.0] +fov_dim = 1 +batch_dim = 0 + +min_data_val = 0 +max_data_val = 1 + +# figure shape +for fovi, fov in enumerate(FOVs): + fig_shape = [len(titles), len(conditions)] + fig, ax = plt.subplots(fig_shape[1], fig_shape[0], + constrained_layout=True, figsize=(fig_shape[0] * 3, fig_shape[1] * 3), + sharey=True, sharex=True) + plt.suptitle(f"FOV={fov}") + gs1 = gridspec.GridSpec(fig_shape[0], fig_shape[1]) + gs1.update(wspace=0, hspace=0) + for wi in range(fig_shape[0]): # metrics + metric_min = None + metric_max = None + + for hi in range(fig_shape[1]): # conditions (wall vs infinite) + raw_data = np.load(os.path.join(data_path, file_names[wi] + f"_{conditions[hi]}.npy")) + print(raw_data.shape) + if file_names[wi] == "polarization": + # mean over all agents and runs which is first and last dimesnions + mean_data = np.mean(raw_data[:, fovi, ...], axis=0) + elif file_names[wi] == "meaniid": + mean_data = np.mean(raw_data[fovi, ...], axis=-1) + else: + mean_data = np.mean(np.mean(raw_data[fovi, ...], axis=-1), axis=0) + print(f"mean shape: ", mean_data.shape) + + if metric_max is None: + metric_max = np.max(mean_data) + else: + metric_max = max(np.max(mean_data), metric_max) + + if metric_min is None: + metric_min = np.min(mean_data) + else: + metric_min = min(np.min(mean_data), metric_min) + + + hi = 0 + for hi in range(fig_shape[1]): # conditions (wall vs infinite) + + plt.axes(ax[hi, wi]) + if hi == 0: + plt.title(titles[wi]) + elif hi == fig_shape[1]-1 and wi == 0: + plt.xticks([i for i in range(len(betas))], betas) + plt.xlabel(f"Beta0") + if wi == 0: + plt.yticks([i for i in range(len(alphas))], alphas) + plt.ylabel(f"{conditions[hi]}\nAlpha0") + + + print(os.path.join(data_path, file_names[wi]+f"_{conditions[hi]}.npy")) + raw_data = np.load(os.path.join(data_path, file_names[wi]+f"_{conditions[hi]}.npy")) + print(raw_data.shape) + if file_names[wi] == "polarization": + # mean over all agents and runs which is first and last dimesnions + mean_data = np.mean(raw_data[:, fovi, ...], axis=0) + elif file_names[wi] == "meaniid": + mean_data = np.mean(raw_data[fovi, ...], axis=-1) + else: + mean_data = np.mean(np.mean(raw_data[fovi, ...], axis=-1), axis=0) + if mean_data.ndim == 2: + im = plt.imshow(mean_data, vmin=metric_min, vmax=metric_max) + fig.colorbar(im) + else: + print("resulting data is of non image shape") + + + + + # if ni < len(Ns): + # print(hi, wi, ni) + # plt.axes(ax[wi, hi]) + # collapsed_data = np.load(os.path.join(data_path, f"coll_eff_N{Ns[ni]}.npy")) + # # column-wise normalization + # for coli in range(collapsed_data.shape[1]): + # print(f"Normalizing column {coli}") + # minval = np.min(collapsed_data[:, coli]) + # maxval = np.max(collapsed_data[:, coli]) + # collapsed_data[:, coli] = (collapsed_data[:, coli] - minval) / (maxval - minval) + # agent_dim = 4 + # im = plt.imshow(collapsed_data, vmin=min_data_val, vmax=max_data_val) + # + # # creating x-axis and labels + # # reading original labels with long names + # labels = np.loadtxt(os.path.join(data_path, f"coll_eff_N{Ns[ni]}.txt"), dtype=str) + # + # # simplifying and concatenating labels + # conc_x_labels = [] + # for li in range(0, len(labels), 2): + # conc_x_labels.append(labels[li + 1].replace("N_RESOURCES=", "").replace(".0,", " x ") + + # labels[li].replace("MIN_RESOURCE_PER_PATCH=", "").replace(".0", "U")) + # + # # creating y-axis labels + # tuned_env_pattern = os.path.join(data_path, "tuned_env*.json") + # print("Patterns: ", tuned_env_pattern) + # json_files = glob.glob(tuned_env_pattern) + # for json_path in json_files: + # with open(json_path, "r") as f: + # envvars = json.load(f) + # y_labels = envvars["DEC_EPSW"] + # + # # Manually definig the sparified y ticks and their rotations + # sparse_y_indices = [0, 2, 4, 6, 8] # int((len(y_labels)-1)/2), len(y_labels)-1] + # y_tick_rotations = [45, 33.75, 22.5, 11.25, 0, -11.25, -22.5, -33.75, -45] + # sparese_y_tick_rotations = [45, 0, 0, 0, -45] + # + # # The first plot is detailed shows all ticks and labels + # if hi == 0 and wi == 0: + # plt.ylabel("Social Excitability ($\epsilon_w$)") + # plt.yticks([i for i in range(len(y_labels))], y_labels, ha='right', rotation_mode='anchor') + # # for ticki, tick in enumerate(ax[wi, hi].get_yticklabels()): + # # print(ticki, tick) + # # tick.set_rotation(y_tick_rotations[ticki]) + # + # # The others are sparsified for better readability + # elif hi == 0: + # plt.yticks([i for i in sparse_y_indices], + # [y_labels[i] for i in sparse_y_indices], ha='right', rotation_mode='anchor') + # for ticki, tick in enumerate(ax[wi, hi].get_yticklabels()): + # tick.set_rotation(sparese_y_tick_rotations[ticki]) + # + # # The ones in second or other columns have same y axis + # else: + # plt.yticks([], []) + # + # if hi == 0 and wi == fig_shape[0] - 1: + # plt.xlabel("Environment") + # plt.xticks([i for i in range(0, len(conc_x_labels), 2)], + # [conc_x_labels[i] for i in range(0, len(conc_x_labels), 2)], + # ha='right', rotation_mode='anchor') + # elif wi == fig_shape[0]-1: + # plt.xticks([i for i in range(0, len(conc_x_labels), 2)], + # [conc_x_labels[i] for i in range(0, len(conc_x_labels), 2)], + # ha='right', rotation_mode='anchor') + # for tick in ax[wi, hi].get_xticklabels(): + # tick.set_rotation(45) + # else: + # plt.axes(ax[wi, hi]) + # plt.yticks([], []) + # plt.xticks([], []) + +# # Adding the colorbar +# # [left, bottom, width, height] +# cbaxes = fig.add_axes([0.2, 0.805, 0.6, 0.03]) +# # position for the colorbar +# cb = plt.colorbar(im, cax=cbaxes, orientation='horizontal') +# cb.ax.set_xticks([min_data_val, (min_data_val+max_data_val)/2, max_data_val]) +# cb.ax.xaxis.set_ticks_position("top") +# cb.ax.xaxis.set_label_position('top') +# cb.ax.set_title('Relative Search Efficiency') +# # fig.colorbar(im, ax=ax.ravel().tolist(), orientation = 'horizontal') +# plt.subplots_adjust(hspace=0, wspace=0, top=0.8, bottom=0.2, left=0.2, right=0.8) +plt.show() \ No newline at end of file diff --git a/abm/loader/data_loader.py b/abm/loader/data_loader.py index 3c97d461..16dfd077 100644 --- a/abm/loader/data_loader.py +++ b/abm/loader/data_loader.py @@ -116,6 +116,7 @@ def __init__(self, data_folder_path, only_env=False, only_agent=False, undersamp self.load_files() self.preprocess_data() + def agent_json_to_csv_format(self): """transforming a read in json format dictionary to the standard data structure we use when read in a csv file""" new_dict = {} @@ -248,29 +249,32 @@ def load_files(self): print("agent_data loaded") if not self.only_agent: - if self.resource_csv_path is not None: - self.resource_data = dh.load_csv_file(self.resource_csv_path, undersample=self.undersample) - print("OLD") - elif self.resource_json_path is not None: - with open(self.resource_json_path, "r") as f: - self.resource_data = json.load(f) - self.resource_data = self.resource_json_to_csv_format() - else: - if self.zarr_compressed_runs: - self.resource_data = {} - self.resource_data['posx'] = zarr.open(os.path.join(self.data_folder_path, f"res_posx{self.zarr_extension}"), - mode='r') - self.resource_data['posy'] = zarr.open(os.path.join(self.data_folder_path, f"res_posy{self.zarr_extension}"), - mode='r') - self.resource_data['radius'] = zarr.open( - os.path.join(self.data_folder_path, f"res_rad{self.zarr_extension}"), - mode='r') - if self.project_version=="Base": - self.resource_data['resc_left'] = zarr.open( - os.path.join(self.data_folder_path, f"res_left{self.zarr_extension}"), + try: + if self.resource_csv_path is not None: + self.resource_data = dh.load_csv_file(self.resource_csv_path, undersample=self.undersample) + print("OLD") + elif self.resource_json_path is not None: + with open(self.resource_json_path, "r") as f: + self.resource_data = json.load(f) + self.resource_data = self.resource_json_to_csv_format() + else: + if self.zarr_compressed_runs: + self.resource_data = {} + self.resource_data['posx'] = zarr.open(os.path.join(self.data_folder_path, f"res_posx{self.zarr_extension}"), + mode='r') + self.resource_data['posy'] = zarr.open(os.path.join(self.data_folder_path, f"res_posy{self.zarr_extension}"), + mode='r') + self.resource_data['radius'] = zarr.open( + os.path.join(self.data_folder_path, f"res_rad{self.zarr_extension}"), mode='r') - self.resource_data['quality'] = zarr.open(os.path.join(self.data_folder_path, f"res_qual{self.zarr_extension}"), - mode='r') + if self.project_version=="Base": + self.resource_data['resc_left'] = zarr.open( + os.path.join(self.data_folder_path, f"res_left{self.zarr_extension}"), + mode='r') + self.resource_data['quality'] = zarr.open(os.path.join(self.data_folder_path, f"res_qual{self.zarr_extension}"), + mode='r') + except: + pass print("resource data loaded") @@ -302,7 +306,7 @@ def preprocess_data(self): self.env_data["SHOW_VISUAL_FIELDS"] = bool(int(self.env_data["SHOW_VISUAL_FIELDS"])), self.env_data["POOLING_TIME"] = int(self.env_data["POOLING_TIME"]), self.env_data["POOLING_PROBABILITY"] = float(self.env_data["POOLING_PROBABILITY"]), - self.env_data["RADIUS_AGENT"] = int(self.env_data["RADIUS_AGENT"]), + self.env_data["RADIUS_AGENT"] = int(float(self.env_data["RADIUS_AGENT"])), self.env_data["N_RESOURCES"] = int(self.env_data["N_RESOURCES"]), self.env_data["MIN_RESOURCE_PER_PATCH"] = int(self.env_data["MIN_RESOURCE_PER_PATCH"]), self.env_data["MAX_RESOURCE_PER_PATCH"] = int(self.env_data["MAX_RESOURCE_PER_PATCH"]), @@ -361,9 +365,12 @@ class ExperimentLoader: def __init__(self, experiment_path, enforce_summary=False, undersample=1, with_plotting=False, collapse_plot=None, t_start=None, t_end=None): # experiment data after summary + self.clustering_data = None + self.clustering_ids = None self.project_version = None self.zarr_extension = ".zarr" self.mean_iid = None + self.mean_nn_dist = None self.iid_matrix = None self.undersample = int(undersample) self.chunksize = None # chunk size in zarr array @@ -790,7 +797,10 @@ def read_all_data(self, only_res=False, project_version="Base"): r_rescleft_array[ind] += data else: - num_res_in_run = res_data['posx'].shape[0] + try: + num_res_in_run = res_data['posx'].shape[0] + except: + num_res_in_run = 0 for pid in range(num_res_in_run): ind = (i,) + tuple(index) + (pid,) r_posx_array[ind] = res_data['posx'][pid, self.t_start:self.t_end:self.undersample] @@ -897,10 +907,17 @@ def reload_summarized_data(self): self.agent_summary['mode'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_mode{extension}"), mode='r') # self.experiment.agent_summary['mode'] if self.project_version=="Base": - self.agent_summary['u'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_u{extension}"), - mode='r') - self.agent_summary['w'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_w{extension}"), - mode='r') + # u and w are not crucial for replay, so we can skip them if they are not available + try: + self.agent_summary['u'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_u{extension}"), + mode='r') + except zarr.errors.PathNotFoundError: + print("No agent_u found!") + try: + self.agent_summary['w'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_w{extension}"), + mode='r') + except zarr.errors.PathNotFoundError: + print("No agent_w found!") self.agent_summary['explpatch'] = zarr.open(os.path.join(self.experiment_path, "summary", f"agent_explpatch{extension}"), mode='r') self.agent_summary['collresource'] = zarr.open( @@ -937,8 +954,12 @@ def reload_summarized_data(self): self.res_summary['resc_left'] = zarr.open( os.path.join(self.experiment_path, "summary", "res_rescleft.zarr"), mode='r') - self.res_summary['quality'] = zarr.open(os.path.join(self.experiment_path, "summary", "res_qual.zarr"), - mode='r') + try: + self.res_summary['quality'] = zarr.open(os.path.join(self.experiment_path, "summary", "res_qual.zarr"), + mode='r') + except zarr.errors.PathNotFoundError: + print("No res_qual found!") + else: # found npz summary for resources @@ -966,6 +987,310 @@ def reload_summarized_data(self): print("___END_README___\n") print("Experiment loaded") + def calculate_pairwise_pol_matrix_supervect(self, t, undersample=1, batchi=0): + """ + Calculating NxN matrix of pairwise polarizations between agents in a given condition. + The condition idx is a tuple of the varying parameter indices. + t is the time index. + """ + # Defining dimesnions in the orientation array + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + + # reshape the orientation array to have the first dimensions batch_dim and num_var_params + print("Loading Orientation data!") + ori_data = self.agent_summary["orientation"][batchi, ..., ::undersample] + ori_shape = ori_data.shape + num_timesteps = ori_shape[-1] + new_shape = ori_shape[0:agent_dim] + (num_agents, num_timesteps) + print(ori_shape, new_shape) + ori_reshaped = np.reshape(ori_data, new_shape) + + # calculate pairwise polarizations + ag_uni = np.stack((np.cos(ori_reshaped), np.sin(ori_reshaped)), axis=-1) + ag_uni = ag_uni.reshape((-1, num_agents, num_timesteps, 2)) + normed_sum = np.linalg.norm(ag_uni[:, :, None, :] + ag_uni[:, None, :, :], axis=-1) / 2 + pol_matrix = normed_sum.reshape((new_shape[:-1] + (num_agents, num_agents, num_timesteps))) + + # reshape pol_matrix back to the desired shape + pol_matrix = np.moveaxis(pol_matrix, -1, agent_dim) + return pol_matrix + + def calculate_clustering(self): + """Using hierarhical clustering according to inter-individual distance and order scores to get number of + subgroups""" + # Checking if in the previous run a cluster_dict has already been saved and reloading it if exists + clustering_data_path = os.path.join(self.experiment_path, "summary", "clustering_data.npy") + largest_clustering_data_path = os.path.join(self.experiment_path, "summary", "largest_clustering_data.npy") + clustering_id_path = os.path.join(self.experiment_path, "summary", "clustering_id.npy") + if self.iid_matrix is None: + self.calculate_interindividual_distance() + + if os.path.isfile(clustering_data_path): + print("Cluster dict reloaded!") + # loading cluster data numpy array + self.clustering_data = np.load(clustering_data_path, mmap_mode="r+") + self.largest_clustering_data = np.load(largest_clustering_data_path, mmap_mode="r+") + self.clustering_ids = np.load(clustering_id_path, mmap_mode="r+") + return + + from fastcluster import linkage + from scipy.cluster.hierarchy import dendrogram + clustering_dict = {} + num_agents = int(self.env.get("N")) + # calculating the maximum distance on the torus as the half diameter of the arena + max_dist = (self.env.get("ENV_WIDTH") ** 2 + self.env.get("ENV_HEIGHT") ** 2) ** 0.5 / 2 + # The shape will depend on the IID matrix which is costly. + # If we use undersample there, we also use undersample here + num_timesteps_orig = self.env.get("T") + num_timesteps = self.iid_matrix.shape[-1] + undersample = int(num_timesteps_orig // num_timesteps) + clustering_dict['num_subgroups'] = [] + print(f"Num varying params: {len(self.varying_params)}") + num_varying_params = len(self.varying_params) + clustering_data = np.zeros(tuple(list(self.iid_matrix.shape[:num_varying_params+1])+[num_timesteps])) + largest_clustering_data = np.zeros(tuple(list(self.iid_matrix.shape[:num_varying_params+1])+[num_timesteps])) + clustering_ids = np.zeros(tuple(list(self.iid_matrix.shape[:num_varying_params+1])+[num_agents, num_timesteps])) + # calculating the number of parameter combinations along iidm.shape[:num_varying_params+1] + num_combinations = np.prod(self.iid_matrix.shape[:num_varying_params+1]) + curr_param_comb = 0 + print("Calculating polarizations first") + t_slice = slice(0, num_timesteps_orig, undersample) + + # calculating the normalized iid distance matrix with punishing large distances more + niidm = self.iid_matrix.copy() + niidm /= max_dist + niidm = (niidm - np.mean(niidm)) / np.std(niidm) + # punishing large distances more + niidm = niidm ** 2 + + for batchi in range(self.agent_summary["orientation"].shape[0]): + print("Batchi: ", batchi) + # if batchi < 5: + # pass + # else: + # break + condition_idx = (batchi, slice(None), slice(None), slice(None)) + pol_m_large = self.calculate_pairwise_pol_matrix_vectorized(condition_idx, t_slice) + for idx_base_ in np.ndindex(*self.iid_matrix.shape[1:num_varying_params+1]): + print("Progress: ", curr_param_comb, "/", num_combinations) + idx_base = tuple([batchi] + list(idx_base_)) + for t in range(num_timesteps): + idx = tuple(list(idx_base) + [slice(None), slice(None), t]) + idx_ = tuple(list(idx_base_) + [slice(None), slice(None), t]) + + # getting the distance matrix for current condition and timestep + normiidm = niidm[idx] + + # calculating the distance matrix for polarization + pm = pol_m_large[idx_] + # standardizing pm and normiidm so that their mean is 0 and their std is 1 + pm = (pm - np.mean(pm)) / np.std(pm) + dist_pm = 1 - pm.astype('float') + # squaring distances to punish large distances more + dist_pm = dist_pm ** 2 + + # calculating the final distance matrix as a weighted average of the two + dist = (dist_pm + normiidm) / 2 + + if idx_base == (0, 1, 5, 2) and 326 < t < 310: + with_plotting = True + print(dist_pm) + print(normiidm) + else: + with_plotting = False + + # clustering + linkage_matrix = linkage(dist, "ward") + ret = dendrogram(linkage_matrix, color_threshold=15, labels=[i for i in range(num_agents)], + show_leaf_counts=True, no_plot=not with_plotting) + colors = [color for _, color in sorted(zip(ret['leaves'], ret['leaves_color_list']))] + group_ids = np.array([int(a.split("C")[1]) for a in colors]) + for i, gid in enumerate(group_ids): + # when gid is zero we make it the maximum element + # because it means it is an independent leaf + if gid == 0: + group_ids[i] = np.max(group_ids) + 1 + + group_ids -= 1 + clustering_data[tuple(list(idx_base) + [t])] = len(list(set(colors))) + largest_clustering_data[tuple(list(idx_base) + [t])] = self.calculate_largest_subcluster_size(group_ids) + clustering_ids[tuple(list(idx_base) + [slice(None), t])] = group_ids + if with_plotting: + plt.show() + input(group_ids) + curr_param_comb += 1 + + self.clustering_data = clustering_data + self.largest_clustering_data = largest_clustering_data + self.clustering_ids = clustering_ids + print(f"Saving clustering data...") + # save clustering_data as numpy array + np.save(clustering_data_path, clustering_data) + np.save(largest_clustering_data_path, largest_clustering_data) + np.save(clustering_id_path, clustering_ids) + + return clustering_dict + + def calculate_largest_subcluster_size(self, cluster_id_list): + """Calculating the size (i.e. number of agents) of largest subcluster of the group""" + cluster_sizes = [] + for cluster_id in set(cluster_id_list): + cluster_sizes.append(np.sum(cluster_id_list == cluster_id)) + return np.max(cluster_sizes) + + def plot_largest_subclusters(self): + """Method to plot size of largest subclusters irrespectively of how many parameters have been tuned during the + experiments.""" + cbar = None + self.calculate_clustering() + self.mean_largest_clusters = np.mean(self.largest_clustering_data, axis=0) + min_data = np.min(self.mean_largest_clusters) + max_data = np.max(self.mean_largest_clusters) + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if num_var_params == 1: + # fig, ax = plt.subplots(1, 1) + # plt.title("Number of Subclusters") + # plt.plot(self.mean_clusters) + # for run_i in range(self.efficiency.shape[0]): + # plt.plot(np.mean(self.efficiency, axis=agent_dim)[run_i, ...], marker=".", linestyle='None') + # ax.set_xticks(range(len(self.varying_params[list(self.varying_params.keys())[0]]))) + # ax.set_xticklabels(self.varying_params[list(self.varying_params.keys())[0]]) + # plt.xlabel(list(self.varying_params.keys())[0]) + pass + + elif num_var_params == 2: + fig, ax = plt.subplots(1, 1) + keys = sorted(list(self.varying_params.keys())) + im = ax.imshow(self.mean_largest_clusters) + + ax.set_yticks(range(len(self.varying_params[keys[0]]))) + ax.set_yticklabels(self.varying_params[keys[0]]) + ax.set_ylabel(keys[0]) + + ax.set_xticks(range(len(self.varying_params[keys[1]]))) + ax.set_xticklabels(self.varying_params[keys[1]]) + ax.set_xlabel(keys[1]) + + elif num_var_params == 3 or num_var_params == 4: + if len(self.mean_largest_clusters.shape) == 4: + # reducing the number of variables to 3 by connecting 2 of the dimensions + self.new_mean_largest_clusters = np.zeros((self.mean_largest_clusters.shape[0:3])) + print(self.new_mean_largest_clusters.shape) + for j in range(self.mean_largest_clusters.shape[0]): + for i in range(self.mean_largest_clusters.shape[1]): + self.new_mean_largest_clusters[j, i, :] = self.mean_largest_clusters[j, i, :, i] + self.mean_largest_clusters = self.new_mean_largest_clusters + if self.collapse_plot is None: + num_plots = self.mean_largest_clusters.shape[0] + fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + for i in range(num_plots): + img = ax[i].imshow(self.mean_largest_clusters[i, :, :], vmin=min_data, vmax=max_data) + ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") + + if i == 0: + ax[i].set_yticks(range(len(self.varying_params[keys[1]]))) + ax[i].set_yticklabels(self.varying_params[keys[1]]) + ax[i].set_ylabel(keys[1]) + + ax[i].set_xticks(range(len(self.varying_params[keys[2]]))) + ax[i].set_xticklabels(self.varying_params[keys[2]]) + ax[i].set_xlabel(keys[2]) + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + description_text = "" + self.add_plot_interaction(description_text, fig, ax, show=True) + return fig, ax, cbar + + + def plot_clustering(self): + """Method to plot clustering irrespectively of how many parameters have been tuned during the + experiments.""" + cbar = None + self.calculate_clustering() + self.mean_clusters = np.mean(self.clustering_data, axis=0) + min_data = np.min(self.mean_clusters) + max_data = np.max(self.mean_clusters) + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if num_var_params == 1: + # fig, ax = plt.subplots(1, 1) + # plt.title("Number of Subclusters") + # plt.plot(self.mean_clusters) + # for run_i in range(self.efficiency.shape[0]): + # plt.plot(np.mean(self.efficiency, axis=agent_dim)[run_i, ...], marker=".", linestyle='None') + # ax.set_xticks(range(len(self.varying_params[list(self.varying_params.keys())[0]]))) + # ax.set_xticklabels(self.varying_params[list(self.varying_params.keys())[0]]) + # plt.xlabel(list(self.varying_params.keys())[0]) + pass + + elif num_var_params == 2: + fig, ax = plt.subplots(1, 1) + keys = sorted(list(self.varying_params.keys())) + im = ax.imshow(self.mean_clusters) + + ax.set_yticks(range(len(self.varying_params[keys[0]]))) + ax.set_yticklabels(self.varying_params[keys[0]]) + ax.set_ylabel(keys[0]) + + ax.set_xticks(range(len(self.varying_params[keys[1]]))) + ax.set_xticklabels(self.varying_params[keys[1]]) + ax.set_xlabel(keys[1]) + + elif num_var_params == 3 or num_var_params == 4: + if len(self.mean_clusters.shape) == 4: + # reducing the number of variables to 3 by connecting 2 of the dimensions + self.new_mean_clusters = np.zeros((self.mean_clusters.shape[0:3])) + print(self.new_mean_clusters.shape) + for j in range(self.mean_clusters.shape[0]): + for i in range(self.mean_clusters.shape[1]): + self.new_mean_clusters[j, i, :] = self.mean_clusters[j, i, :, i] + self.mean_clusters = self.new_mean_clusters + if self.collapse_plot is None: + num_plots = self.mean_clusters.shape[0] + fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + for i in range(num_plots): + img = ax[i].imshow(self.mean_clusters[i, :, :], vmin=min_data, vmax=max_data) + ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") + + if i == 0: + ax[i].set_yticks(range(len(self.varying_params[keys[1]]))) + ax[i].set_yticklabels(self.varying_params[keys[1]]) + ax[i].set_ylabel(keys[1]) + + ax[i].set_xticks(range(len(self.varying_params[keys[2]]))) + ax[i].set_xticklabels(self.varying_params[keys[2]]) + ax[i].set_xlabel(keys[2]) + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + description_text = "" + self.add_plot_interaction(description_text, fig, ax, show=True) + return fig, ax, cbar + + def calculate_search_efficiency(self, t_start_plot=0, t_end_plot=-1, used_batches=None): """Method to calculate search efficiency throughout the experiments as the sum of collected resorces normalized with the travelled distance. The timestep in which the efficiency is calculated. This might mismatch from @@ -1039,11 +1364,12 @@ def calculate_interindividual_distance_slice(self, posx, posy): iid[i, j] = distance[i] return iid - def calculate_interindividual_distance(self, undersample=1, avg_over_time=False): + def calculate_interindividual_distance(self, undersample=1, avg_over_time=False, periodic_boundary=False): """Method to calculate inter-individual distance array from posx and posy arrays of agents. The final array has the same dimension as any of the input arrays, i.e.: (num_batches, *[dims of varying params], num_agents, time) and defines the mean (across group members) inter-individual distance for a given agent i in timestep t. + If periodic boundary is true, we calculate distances on a torus. """ summary_path = os.path.join(self.experiment_path, "summary") iidpath = os.path.join(summary_path, "iid.npy") @@ -1052,6 +1378,12 @@ def calculate_interindividual_distance(self, undersample=1, avg_over_time=False) print("Found saved I.I.D array in summary, reloading it...") self.iid_matrix = np.load(iidpath) else: + if self.env.get("BOUNDARY") == "infinite": + print( + "Dataset was generated with periodic boundary conditions. Do you want to calculate I.I.D. according" + "to torus environment? (Y/N)") + periodic_boundary = input().lower() == "y" + agposx = self.agent_summary['posx'] agposy = self.agent_summary['posy'] @@ -1078,7 +1410,7 @@ def calculate_interindividual_distance(self, undersample=1, avg_over_time=False) iid = np.zeros(new_shape) for batchi in range(num_batches): - print(f"Calculating iid for batch {batchi}") + print(f"Calculating iid for batch {batchi}, torus={periodic_boundary}") for agi in range(num_agents): for agj in range(num_agents): if agj > agi: @@ -1086,7 +1418,16 @@ def calculate_interindividual_distance(self, undersample=1, avg_over_time=False) y1s = agposy[batchi, ..., agi, ::undersample] x2s = agposx[batchi, ..., agj, ::undersample] y2s = agposy[batchi, ..., agj, ::undersample] - distance_matrix = supcalc.distance_coords(x1s, y1s, x2s, y2s, vectorized=True) + if not periodic_boundary: + distance_matrix = supcalc.distance_coords(x1s, y1s, x2s, y2s, vectorized=True) + else: + p1 = np.array([[x1s[ii], y1s[ii]] for ii in range(len(x1s))]) + p2 = np.array([[x2s[ii], y2s[ii]] for ii in range(len(x2s))]) + dim = np.array([int(self.env["ENV_WIDTH"]), int(self.env["ENV_HEIGHT"])]) + dima = np.zeros_like(p1) + dima[:, 0, ...] = dim[0] + dima[:, 1, ...] = dim[1] + distance_matrix = supcalc.distance_torus(p1, p2, dima) if not avg_over_time: iid[batchi, ..., agi, agj, :] = distance_matrix else: @@ -1101,6 +1442,7 @@ def calculate_interindividual_distance(self, undersample=1, avg_over_time=False) # for the mean we restrict the iid matrix as upper triangular in the agent dimensions so that we # don't calculate IIDs in mean twice (as IID is symmetric, i.e. the distance between agent i and j # is the same as between j and i) + num_agents = self.agent_summary['posx'].shape[-2] restr_m = self.iid_matrix[..., np.triu_indices(num_agents, k=1)[0], np.triu_indices(num_agents, k=1)[1], :] # Then we take the mean along the flattened dimension in which we defined the previous restriction, and we also # take the mean along all the repeated batches (0dim) @@ -1111,10 +1453,452 @@ def calculate_interindividual_distance(self, undersample=1, avg_over_time=False) else: self.mean_iid = np.mean(np.mean(restr_m, axis=-2), axis=0) - # Saving calculated arrays for future use - print("Saving I.I.D and mean I.I.D arrays into summary!") - np.save(iidpath, self.iid_matrix) - np.save(meaniid_path, self.mean_iid) + # Saving calculated arrays for future use + print("Saving I.I.D and mean I.I.D arrays into summary!") + np.save(iidpath, self.iid_matrix) + np.save(meaniid_path, self.mean_iid) + + def calculate_mean_NN_dist(self, undersample=1, avg_over_time=False): + """Calculating the mean nearest neighbor metric over time. We take the average over + agents of the min neighbor distance per agent.""" + summary_path = os.path.join(self.experiment_path, "summary") + iidpath = os.path.join(summary_path, "iid.npy") + meanNNd_path = os.path.join(summary_path, "meanNNd.npy") + if os.path.isfile(meanNNd_path): + print("Found saved mean NND array in summary, reloading it...") + + self.mean_nn_dist = np.load(meanNNd_path) + if self.iid_matrix is None: + if os.path.isfile(iidpath): + print("Found saved IID array in summary, reloading it...") + self.iid_matrix = np.load(iidpath) + else: + print("Didn't find saved IID array in summary, calculating them...") + self.calculate_interindividual_distance(undersample=undersample, avg_over_time=False) + iid = self.iid_matrix.copy() + num_agents = iid.shape[-2] + for i in range(num_agents): + iid[..., i, i, :] = np.nan + nearest_per_agents = np.nanmin(iid, axis=-2) + mean_nearest = np.mean(nearest_per_agents, axis=-2) + self.mean_nn_dist = np.mean(mean_nearest, axis=0) + print("Saving mean NNd array under ", meanNNd_path) + np.save(meanNNd_path, self.mean_nn_dist) + + + + def plot_mean_polarization(self, t_start=0, t_end=-1, from_script=False, used_batches=None): + """Method to plot polarization irrespectively of how many parameters have been tuned during the + experiments.""" + cbar = None + self.calculate_polarization() + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if num_var_params == 1: + fig, ax = plt.subplots(1, 1) + plt.title("Polarization") + plt.plot(self.mean_pol) + plt.plot(self.mean_pol + self.pol_std) + plt.plot(self.mean_pol - self.pol_std) + for run_i in range(self.efficiency.shape[0]): + plt.plot(np.mean(self.efficiency, axis=agent_dim)[run_i, ...], marker=".", linestyle='None') + ax.set_xticks(range(len(self.varying_params[list(self.varying_params.keys())[0]]))) + ax.set_xticklabels(self.varying_params[list(self.varying_params.keys())[0]]) + plt.xlabel(list(self.varying_params.keys())[0]) + + elif num_var_params == 2: + fig, ax = plt.subplots(1, 1) + keys = sorted(list(self.varying_params.keys())) + im = ax.imshow(self.mean_pol) + + ax.set_yticks(range(len(self.varying_params[keys[0]]))) + ax.set_yticklabels(self.varying_params[keys[0]]) + ax.set_ylabel(keys[0]) + + ax.set_xticks(range(len(self.varying_params[keys[1]]))) + ax.set_xticklabels(self.varying_params[keys[1]]) + ax.set_xlabel(keys[1]) + + elif num_var_params == 3 or num_var_params == 4: + if len(self.mean_pol.shape) == 4: + # reducing the number of variables to 3 by connecting 2 of the dimensions + self.new_mean_pol = np.zeros((self.mean_pol.shape[0:3])) + print(self.new_mean_pol.shape) + for j in range(self.mean_pol.shape[0]): + for i in range(self.mean_pol.shape[1]): + self.new_mean_pol[j, i, :] = self.mean_pol[j, i, :, i] + self.mean_pol = self.new_mean_pol + if self.collapse_plot is None: + num_plots = self.mean_pol.shape[0] + fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + for i in range(num_plots): + img = ax[i].imshow(self.mean_pol[i, :, :], vmin=0, vmax=1) + ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") + + if i == 0: + ax[i].set_yticks(range(len(self.varying_params[keys[1]]))) + ax[i].set_yticklabels(self.varying_params[keys[1]]) + ax[i].set_ylabel(keys[1]) + + ax[i].set_xticks(range(len(self.varying_params[keys[2]]))) + ax[i].set_xticklabels(self.varying_params[keys[2]]) + ax[i].set_xlabel(keys[2]) + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + else: + fig, ax = plt.subplots(1, 1, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + + collapsed_data, labels = self.collapse_mean_data(self.mean_pol, save_name="coll_pol.npy") + coll_std, _ = self.collapse_mean_data(self.pol_std, save_name="coll_polstd.npy") + + # # column-wise normalization + # for coli in range(collapsed_data.shape[1]): + # print(f"Normalizing column {coli}") + # minval = np.min(collapsed_data[:, coli]) + # maxval = np.max(collapsed_data[:, coli]) + # collapsed_data[:, coli] = (collapsed_data[:, coli] - minval) / (maxval - minval) + + img = ax.imshow(collapsed_data) + ax.set_yticks(range(len(self.varying_params[keys[self.collapse_fixedvar_ind]]))) + ax.set_yticklabels(self.varying_params[keys[self.collapse_fixedvar_ind]]) + ax.set_ylabel(keys[self.collapse_fixedvar_ind]) + + ax.set_xticks(range(len(labels))) + ax.set_xticklabels(labels, rotation=45) + ax.set_xlabel("Combined Parameters") + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + fig.set_tight_layout(True) + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + description_text = "" + self.add_plot_interaction(description_text, fig, ax, show=True, from_script=from_script) + return fig, ax, cbar + + def plot_mean_collision_time(self, t_start=0, t_end=-1, from_script=False, used_batches=None, undersample=1): + cbar = None + self.calculate_collision_time(undersample=undersample) + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if num_var_params == 1: + fig, ax = plt.subplots(1, 1) + plt.title("Collision Time") + plt.plot(self.mean_aacoll) + plt.plot(self.mean_aacoll + self.aacoll_std) + plt.plot(self.mean_aacoll - self.aacoll_std) + for run_i in range(self.efficiency.shape[0]): + plt.plot(np.mean(self.efficiency, axis=agent_dim)[run_i, ...], marker=".", linestyle='None') + ax.set_xticks(range(len(self.varying_params[list(self.varying_params.keys())[0]]))) + ax.set_xticklabels(self.varying_params[list(self.varying_params.keys())[0]]) + plt.xlabel(list(self.varying_params.keys())[0]) + + elif num_var_params == 2: + fig, ax = plt.subplots(1, 1) + keys = sorted(list(self.varying_params.keys())) + im = ax.imshow(self.mean_aacoll) + + ax.set_yticks(range(len(self.varying_params[keys[0]]))) + ax.set_yticklabels(self.varying_params[keys[0]]) + ax.set_ylabel(keys[0]) + + ax.set_xticks(range(len(self.varying_params[keys[1]]))) + ax.set_xticklabels(self.varying_params[keys[1]]) + ax.set_xlabel(keys[1]) + + elif num_var_params == 3 or num_var_params == 4: + if len(self.mean_aacoll.shape) == 4: + # reducing the number of variables to 3 by connecting 2 of the dimensions + self.new_mean_pol = np.zeros((self.mean_aacoll.shape[0:3])) + print(self.new_mean_pol.shape) + for j in range(self.mean_aacoll.shape[0]): + for i in range(self.mean_aacoll.shape[1]): + self.new_mean_pol[j, i, :] = self.mean_aacoll[j, i, :, i] + self.mean_aacoll = self.new_mean_pol + if self.collapse_plot is None: + num_plots = self.mean_aacoll.shape[0] + fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + for i in range(num_plots): + img = ax[i].imshow(self.mean_aacoll[i, :, :], vmin=0, vmax=np.max(self.mean_aacoll)) + ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") + + if i == 0: + ax[i].set_yticks(range(len(self.varying_params[keys[1]]))) + ax[i].set_yticklabels(self.varying_params[keys[1]]) + ax[i].set_ylabel(keys[1]) + + ax[i].set_xticks(range(len(self.varying_params[keys[2]]))) + ax[i].set_xticklabels(self.varying_params[keys[2]]) + ax[i].set_xlabel(keys[2]) + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + else: + fig, ax = plt.subplots(1, 1, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + + collapsed_data, labels = self.collapse_mean_data(self.mean_aacoll, save_name="coll_aacoll.npy") + coll_std, _ = self.collapse_mean_data(self.aacoll_std, save_name="coll_polstd.npy") + + # # column-wise normalization + # for coli in range(collapsed_data.shape[1]): + # print(f"Normalizing column {coli}") + # minval = np.min(collapsed_data[:, coli]) + # maxval = np.max(collapsed_data[:, coli]) + # collapsed_data[:, coli] = (collapsed_data[:, coli] - minval) / (maxval - minval) + + img = ax.imshow(collapsed_data) + ax.set_yticks(range(len(self.varying_params[keys[self.collapse_fixedvar_ind]]))) + ax.set_yticklabels(self.varying_params[keys[self.collapse_fixedvar_ind]]) + ax.set_ylabel(keys[self.collapse_fixedvar_ind]) + + ax.set_xticks(range(len(labels))) + ax.set_xticklabels(labels, rotation=45) + ax.set_xlabel("Combined Parameters") + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + fig.set_tight_layout(True) + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + description_text = "" + self.add_plot_interaction(description_text, fig, ax, show=True, from_script=from_script) + return fig, ax, cbar + + def calculate_pairwise_pol_matrix_vectorized(self, condition_idx, t): + """t is also slice""" + # Defining dimesnions in the orientation array + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + + # Get the orientations of all agents at time t + agent_ori = self.agent_summary["orientation"][condition_idx + (slice(None), t)] + # print(agent_ori.shape) + # input() + + # Calculate the univectors of all agent pairs + agent_uni = np.stack([np.cos(agent_ori), np.sin(agent_ori)], axis=-1) + # print(agent_uni.shape) + # input() + if isinstance(t, slice): + # print("t is slice") + agi_uni = np.expand_dims(agent_uni, axis=-3) + # print(agi_uni.shape) + agj_uni = np.expand_dims(agent_uni, axis=-4) + # print(agj_uni.shape) + # input() + else: + # print("t is not slice") + agi_uni = np.expand_dims(agent_uni, axis=-2) + # print(agi_uni.shape) + agj_uni = np.expand_dims(agent_uni, axis=-3) + # print(agj_uni.shape) + # input() + + # Calculate the pairwise polarizations using matrix multiplication + pol_matrix = np.linalg.norm(agi_uni + agj_uni, axis=-1) / 2 + + # print(pol_matrix.shape) + # input() + # print(pol_matrix[..., 0]) + + return pol_matrix + + def calculate_pairwise_pol_matrix(self, condition_idx, t): + """ + Calculating NxN matrix of pairwise polarizations between agents in a given condition. + The condition idx is a tuple of the varying parameter indices. + t is the time index. + """ + # Defining dimesnions in the orientation array + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + num_agents = self.agent_summary["orientation"].shape[agent_dim] + + # calculate pairwise polarizations + pol_matrix = np.zeros((num_agents, num_agents)) + for i in range(num_agents): + for j in range(num_agents): + # getting orientation of 2 agents i and j + agi_ori = self.agent_summary["orientation"][condition_idx + (i, t)] + agj_ori = self.agent_summary["orientation"][condition_idx + (j, t)] + + # calculating univectors with orientations of agent i and j + agi_uni = np.array([np.cos(agi_ori), np.sin(agi_ori)]) + agj_uni = np.array([np.cos(agj_ori), np.sin(agj_ori)]) + + # calculating the absolute sum of the univectors normed by the number of agents + # contributing to the sum (=2, pairwise) + normed_sum = np.linalg.norm(agi_uni + agj_uni) / 2 + + pol_matrix[i, j] = normed_sum + + return pol_matrix + + def calculate_polarization(self, undersample=1, filtered_by_wallcoll=0, filtering_window=50): + """Calculating polarization of agents in the environment used to + quantify e.g. flocking models""" + summary_path = os.path.join(self.experiment_path, "summary") + if filtered_by_wallcoll: + self.find_wall_collisions() + polpath = os.path.join(summary_path, f"polarization_wallcoll.npy") + else: + polpath = os.path.join(summary_path, f"polarization_us{undersample}.npy") + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if os.path.isfile(polpath): + print("Found saved polarization array in summary, reloading it...") + self.pol_matrix = np.load(polpath) + self.mean_pol = np.nanmean(np.nanmean(self.pol_matrix, axis=batch_dim), axis=-1) + self.pol_std = np.nanmean(np.nanstd(self.pol_matrix, axis=batch_dim), axis=-1) + + else: + num_agents = self.agent_summary["orientation"].shape[agent_dim] + num_timesteps = self.agent_summary["orientation"].shape[time_dim]/undersample + ori_shape = list(self.agent_summary["orientation"].shape) + new_shape = ori_shape[0:num_var_params+1] + [int(num_timesteps/undersample)] + + self.pol_matrix = np.zeros(new_shape) + unitvec_shape = ori_shape[1:-2] + [2] + [num_agents, int(num_timesteps/undersample)] + + for runi in range(self.num_batches): + print(f"Calculating polarization for batch {runi}") + if filtered_by_wallcoll: + + # wrefs = list(tuple(self.wrefs[runi][:])) + # for ti, t in enumerate(wrefs_old[-1]): + # print(ti) + # wrefs[-1] = np.append(wrefs[1], np.array([tk for tk in range(t, t+100)])) + # for wi in range(len(wrefs)-1): + # wrefs[wi] = np.append(wrefs[wi], np.array([wrefs[wi][t] for tk in range(100)])) + + orif = self.agent_summary["orientation"][runi, ...] + print(f"Extending collision filtering with time window {filtering_window}") + wrefs = self.wrefs[runi][:] + for wi in range(filtering_window): + print(wi) + new_times = wrefs[-1] + 1 + new_times[new_times>self.env["T"]-1]=self.env["T"]-1 + wrefs[-1] = new_times + new_wrefs = tuple(wrefs) + orif[new_wrefs] = np.nan + + unitvecs = np.zeros(unitvec_shape) + for robi in range(num_agents): + if filtered_by_wallcoll: + ori = orif[..., robi, :] + else: + ori = self.agent_summary["orientation"][runi, ..., robi, ::undersample] + unitvecs[..., 0, robi, :] = np.array([np.cos(ang) for ang in ori]) + unitvecs[..., 1, robi, :] = np.array([np.sin(ang) for ang in ori]) + + unitsum = np.nansum(unitvecs, axis=-2) # summing for all robots + unitsum_norm = np.linalg.norm(unitsum, axis=-2) / num_agents # getting norm in x and y + self.pol_matrix[runi, ...] = unitsum_norm # np.nanmean(unitsum_norm, axis=-1) + + # for runi in range(self.num_batches): + # pol_matrix[runi, ...] = np.mean(np.array( + # [np.linalg.norm([unitsum[runi, 0, t], unitsum[runi, 1, t]]) / num_agents for t in + # range(num_timesteps)]), axis=-1) + + self.mean_pol = np.nanmean(np.nanmean(self.pol_matrix, axis=batch_dim), axis=-1) + self.pol_std = np.nanmean(np.nanstd(self.pol_matrix, axis=batch_dim), axis=-1) + print("Saving calculated polarization time matrix!") + np.save(polpath, self.pol_matrix) + + return self.pol_matrix, self.mean_pol + + def calculate_collision_time(self, undersample=1): + summary_path = os.path.join(self.experiment_path, "summary") + iidpath = os.path.join(summary_path, "iid.npy") + aacollpath = os.path.join(summary_path, "aacoll.npy") + batch_dim = 0 + + if os.path.isfile(aacollpath): + print("Found saved collision time array in summary, reloading it...") + self.aacoll_matrix = np.load(aacollpath) + else: + if self.iid_matrix is None and not os.path.isfile(iidpath): + self.calculate_interindividual_distance(undersample=undersample) + elif os.path.isfile(iidpath): + self.iid_matrix = np.load(iidpath) + aacoll = np.zeros(list(self.agent_summary['orientation'].shape)[0:-2]) + num_timesteps = self.iid_matrix.shape[-1] + if self.iid_matrix.shape[-1] < 1000: + raise Exception(f"The IID matrix has been collapsed on the time axis to len {self.iid_matrix.shape[-1]} when was calculated, can not" + " calculate collision matrix. Please delete iid.npy from the summary folder and" + "try again!") + else: + iid_individual_coll = np.any(self.iid_matrix < (2*self.env["RADIUS_AGENT"]), axis=-2) # containing info about which agent has been collided + iid_sum_coll = np.any(np.logical_and(self.iid_matrix > 0, self.iid_matrix < (2*self.env["RADIUS_AGENT"])), axis=-2) # is there collision or not in every timestep + aacoll = np.count_nonzero(iid_sum_coll, axis=-1) / num_timesteps + self.aacoll_matrix = aacoll + + self.mean_aacoll = np.mean(np.mean(self.aacoll_matrix, axis=batch_dim), axis=-1) + self.aacoll_std = np.std(np.mean(self.aacoll_matrix, axis=batch_dim), axis=-1) + print("Saving calculated mean aacoll time matrix!") + np.save(aacollpath, self.aacoll_matrix) + + return self.aacoll_matrix, self.mean_aacoll + + def find_wall_collisions(self, undersample=1): + """Finding the timepoints of the data where ANY of the agents have been reflected from ANY of the walls""" + summary_path = os.path.join(self.experiment_path, "summary") + wrefpaths = [os.path.join(summary_path, f"wallref_b{bi}.zarr") for bi in range(self.num_batches)] + self.wrefs = {} + if os.path.isdir(wrefpaths[0]): + print("Found saved wall reflection array in summary, reloading it...") + for bi in range(self.num_batches): + self.wrefs[bi] = zarr.open(wrefpaths[bi], mode='r', dtype='int') + else: + boundaries_x = [0, int(float(self.env.get("ENV_WIDTH")))] + boundaries_y = [0, int(float(self.env.get("ENV_HEIGHT")))] + for bi in range(self.num_batches): + print(f"Calculating wall reflection times to batch {bi}") + agposx = self.agent_summary['posx'][bi, ...] + agposy = self.agent_summary['posy'][bi, ...] + + wrefs = np.array(np.where(np.logical_or(np.logical_or( + agposx < boundaries_x[0] - self.env.get("RADIUS_AGENT"), + agposx > boundaries_x[1] - self.env.get("RADIUS_AGENT")), + np.logical_or( + agposy < boundaries_y[0] - self.env.get("RADIUS_AGENT"), + agposy > boundaries_y[1] - self.env.get("RADIUS_AGENT"))))) + + wrefszarr = zarr.open(wrefpaths[bi], mode='w', + shape=wrefs.shape, dtype='int') + wrefszarr[...] = wrefs[...] + self.wrefs[bi] = wrefszarr + + + def calculate_relocation_time(self, undersample=1): """Calculating relocation time matrix over the given data""" @@ -1189,13 +1973,109 @@ def collapse_mean_data(self, mean_data, save_name=None): return collapsed_data, labels + def plot_mean_NN_dist(self, from_script=False, undersample=1): + """Method to plot mean inter-individual distance irrespectively of how many parameters have been tuned during the + experiments.""" + cbar = None + if self.mean_nn_dist is None: + # self.calculate_interindividual_distance(undersample=undersample) + self.calculate_mean_NN_dist(undersample=undersample) + + batch_dim = 0 + num_var_params = len(list(self.varying_params.keys())) + agent_dim = batch_dim + num_var_params + 1 + time_dim = agent_dim + 1 + + if num_var_params == 1: + fig, ax = plt.subplots(1, 1) + plt.title("Inter-individual distance (mean)") + plt.plot(self.mean_nn_dist) + num_agents = self.iid_matrix.shape[-2] + restr_m = self.iid_matrix[..., np.triu_indices(num_agents, k=1)[0], np.triu_indices(num_agents, k=1)[1]] + for run_i in range(self.iid_matrix.shape[0]): + plt.plot(restr_m[run_i, :, :], marker=".", linestyle='None') + ax.set_xticks(range(len(self.varying_params[list(self.varying_params.keys())[0]]))) + ax.set_xticklabels(self.varying_params[list(self.varying_params.keys())[0]]) + plt.xlabel(list(self.varying_params.keys())[0]) + + elif num_var_params == 2: + fig, ax = plt.subplots(1, 1) + plt.title("Inter-individual distance (mean)") + keys = sorted(list(self.varying_params.keys())) + im = ax.imshow(self.mean_nn_dist) + + ax.set_yticks(range(len(self.varying_params[keys[0]]))) + ax.set_yticklabels(self.varying_params[keys[0]]) + ax.set_ylabel(keys[0]) + + ax.set_xticks(range(len(self.varying_params[keys[1]]))) + ax.set_xticklabels(self.varying_params[keys[1]]) + ax.set_xlabel(keys[1]) + + elif num_var_params == 3 or num_var_params == 4: + if len(self.mean_nn_dist.shape) == 4: + # reducing the number of variables to 3 by connecting 2 of the dimensions + self.new_mean_nn_dist = np.zeros((self.mean_nn_dist.shape[0:3])) + print(self.new_mean_nn_dist.shape) + for j in range(self.mean_nn_dist.shape[0]): + for i in range(self.mean_nn_dist.shape[1]): + self.new_mean_nn_dist[j, i, :] = self.mean_nn_dist[j, i, :, i] + self.mean_nn_dist = self.new_mean_nn_dist + + if self.collapse_plot is None: + num_plots = self.mean_nn_dist.shape[0] + fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + for i in range(num_plots): + img = ax[i].imshow(self.mean_nn_dist[i, :, :]) + ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") + + if i == 0: + ax[i].set_yticks(range(len(self.varying_params[keys[1]]))) + ax[i].set_yticklabels(self.varying_params[keys[1]]) + ax[i].set_ylabel(keys[1]) + + ax[i].set_xticks(range(len(self.varying_params[keys[2]]))) + ax[i].set_xticklabels(self.varying_params[keys[2]]) + ax[i].set_xlabel(keys[2]) + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + else: + fig, ax = plt.subplots(1, 1, sharex=True, sharey=True) + keys = sorted(list(self.varying_params.keys())) + + collapsed_data, labels = self.collapse_mean_data(self.mean_nn_dist, save_name="coll_iid.npy") + + img = ax.imshow(collapsed_data) + ax.set_yticks(range(len(self.varying_params[keys[self.collapse_fixedvar_ind]]))) + ax.set_yticklabels(self.varying_params[keys[self.collapse_fixedvar_ind]]) + ax.set_ylabel(keys[self.collapse_fixedvar_ind]) + + ax.set_xticks(range(len(labels))) + ax.set_xticklabels(labels, rotation=45) + ax.set_xlabel("Combined Parameters") + + fig.subplots_adjust(right=0.8) + cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) + cbar = fig.colorbar(img, cax=cbar_ax) + + fig.set_tight_layout(True) + + num_agents = self.agent_summary["posx"].shape[agent_dim] + description_text = f"Showing the mean (over {self.num_batches} batches and {num_agents} agents)\n" \ + f"of inter-individual distance between agents.\n" + self.add_plot_interaction(description_text, fig, ax, show=True, from_script=from_script) + return fig, ax, cbar + def plot_mean_iid(self, from_script=False, undersample=1): """Method to plot mean inter-individual distance irrespectively of how many parameters have been tuned during the experiments.""" cbar = None if self.iid_matrix is None: - self.calculate_interindividual_distance(avg_over_time=True, undersample=undersample) + self.calculate_interindividual_distance(undersample=undersample) batch_dim = 0 num_var_params = len(list(self.varying_params.keys())) @@ -1243,7 +2123,7 @@ def plot_mean_iid(self, from_script=False, undersample=1): fig, ax = plt.subplots(1, num_plots, sharex=True, sharey=True) keys = sorted(list(self.varying_params.keys())) for i in range(num_plots): - img = ax[i].imshow(self.mean_iid[i, :, :]) + img = ax[i].imshow(self.mean_iid[i, :, :], vmin=np.min(self.mean_iid), vmax=np.max(self.mean_iid)) ax[i].set_title(f"{keys[0]}={self.varying_params[keys[0]][i]}") if i == 0: @@ -1279,7 +2159,7 @@ def plot_mean_iid(self, from_script=False, undersample=1): fig.set_tight_layout(True) - num_agents = self.agent_summary["collresource"].shape[agent_dim] + num_agents = self.agent_summary["posx"].shape[agent_dim] description_text = f"Showing the mean (over {self.num_batches} batches and {num_agents} agents)\n" \ f"of inter-individual distance between agents.\n" self.add_plot_interaction(description_text, fig, ax, show=True, from_script=from_script) @@ -1378,7 +2258,7 @@ def plot_search_efficiency(self, t_start=0, t_end=-1, from_script=False, used_ba cbar = fig.colorbar(img, cax=cbar_ax) fig.set_tight_layout(True) - num_agents = self.agent_summary["collresource"].shape[agent_dim] + num_agents = self.agent_summary["posx"].shape[agent_dim] description_text = f"Showing the mean (over {self.num_batches} batches and {num_agents} agents)\n" \ f"of total collected resource units normalized with the mean total\n" \ f"distance travelled by agents over the experiments.\n" @@ -1475,7 +2355,7 @@ def plot_mean_relocation_time(self): fig.set_tight_layout(True) - num_agents = self.agent_summary["collresource"].shape[agent_dim] + num_agents = self.agent_summary["posx"].shape[agent_dim] description_text = f"Showing the mean (over {self.num_batches} batches and {num_agents} agents)\n" \ f"of relative relocation time, i.e. ratio of time spent in relocation\n" self.add_plot_interaction(description_text, fig, ax, show=True) @@ -1514,7 +2394,7 @@ def plot_mean_travelled_distances(self): ax.set_yticklabels(self.varying_params[list(self.varying_params.keys())[1]]) plt.ylabel(list(self.varying_params.keys())[1]) - num_agents = self.agent_summary["collresource"].shape[agent_dim] + num_agents = self.agent_summary["posx"].shape[agent_dim] description_text = f"Showing the mean (over {self.num_batches} batches and {num_agents} agents)\n" \ f"travelled distance\n" self.add_plot_interaction(description_text, fig, ax, show=True) diff --git a/abm/metarunner/metarunner.py b/abm/metarunner/metarunner.py index 657e67ac..590a4ba7 100644 --- a/abm/metarunner/metarunner.py +++ b/abm/metarunner/metarunner.py @@ -229,6 +229,10 @@ def run_protocol(self, env_path, project="Base"): elif project == "CoopSignaling": from abm import app_collective_signaling app_collective_signaling.start(parallel=self.parallel_run, headless=self.headless) + elif project == "VisualFlocking": + from abm import app_visual_flocking + app_visual_flocking.start(parallel=self.parallel_run, headless=self.headless) + os.remove(default_env_path) shutil.copyfile(backup_default_env, default_env_path) sleep(2) diff --git a/abm/monitoring/ifdb.py b/abm/monitoring/ifdb.py index 9f45cfe6..e280da7d 100644 --- a/abm/monitoring/ifdb.py +++ b/abm/monitoring/ifdb.py @@ -436,34 +436,36 @@ def save_ifdb_as_csv(exp_hash="", use_ram=False, as_zar=True, save_extracted_vfi import zarr, json print("Saving resource data as compressed zarr arrays...") num_res = len(resources_dict) - t_len = len(resources_dict[list(resources_dict.keys())[0]]['pos_x']) - posxzarr = zarr.open(os.path.join(save_dir, "res_posx.zarr"), mode='w', shape=(num_res, t_len), - chunks=(num_res, t_len), dtype='float') - posyzarr = zarr.open(os.path.join(save_dir, "res_posy.zarr"), mode='w', shape=(num_res, t_len), - chunks=(num_res, t_len), dtype='float') - resrad = zarr.open(os.path.join(save_dir, "res_rad.zarr"), mode='w', shape=(num_res, t_len), - chunks=(num_res, t_len), dtype='float') - if project_version == "Base": - rleftzarr = zarr.open(os.path.join(save_dir, "res_left.zarr"), mode='w', shape=(num_res, t_len), - chunks=(num_res, t_len), dtype='float') - qualzarr = zarr.open(os.path.join(save_dir, "res_qual.zarr"), mode='w', shape=(num_res, t_len), + if num_res > 0: + t_len = len(resources_dict[list(resources_dict.keys())[0]]['pos_x']) + posxzarr = zarr.open(os.path.join(save_dir, "res_posx.zarr"), mode='w', shape=(num_res, t_len), chunks=(num_res, t_len), dtype='float') - elif project_version == "CooperativeSignaling": - # Here we can save project specific resource data in individual zarr arrays - pass - for res_id, res_dict in resources_dict.items(): - posxzarr[res_id - 1, :] = resources_dict[res_id]['pos_x'] - posyzarr[res_id - 1, :] = resources_dict[res_id]['pos_y'] - resrad[res_id - 1, :] = [resources_dict[res_id]['radius'] for i in range(t_len)] + posyzarr = zarr.open(os.path.join(save_dir, "res_posy.zarr"), mode='w', shape=(num_res, t_len), + chunks=(num_res, t_len), dtype='float') + resrad = zarr.open(os.path.join(save_dir, "res_rad.zarr"), mode='w', shape=(num_res, t_len), + chunks=(num_res, t_len), dtype='float') if project_version == "Base": - rleftzarr[res_id - 1, :] = resources_dict[res_id]['resc_left'] - qualzarr[res_id - 1, :] = resources_dict[res_id]['quality'] + rleftzarr = zarr.open(os.path.join(save_dir, "res_left.zarr"), mode='w', shape=(num_res, t_len), + chunks=(num_res, t_len), dtype='float') + qualzarr = zarr.open(os.path.join(save_dir, "res_qual.zarr"), mode='w', shape=(num_res, t_len), + chunks=(num_res, t_len), dtype='float') elif project_version == "CooperativeSignaling": # Here we can save project specific resource data in individual zarr arrays pass + for res_id, res_dict in resources_dict.items(): + posxzarr[res_id - 1, :] = resources_dict[res_id]['pos_x'] + posyzarr[res_id - 1, :] = resources_dict[res_id]['pos_y'] + resrad[res_id - 1, :] = [resources_dict[res_id]['radius'] for i in range(t_len)] + if project_version == "Base": + rleftzarr[res_id - 1, :] = resources_dict[res_id]['resc_left'] + qualzarr[res_id - 1, :] = resources_dict[res_id]['quality'] + elif project_version == "CooperativeSignaling": + # Here we can save project specific resource data in individual zarr arrays + pass print(f"Saving agent data as compressed zarr arrays in {save_dir}...") num_ag = len(agents_dict) + t_len = len(agents_dict[list(agents_dict.keys())[0]]['posx']) v_field_len = int(float(ifdbp.envconf.get("VISUAL_FIELD_RESOLUTION"))) aposxzarr = zarr.open(os.path.join(save_dir, "ag_posx.zarr"), mode='w', shape=(num_ag, t_len), chunks=(num_ag, t_len), dtype='float') diff --git a/abm/projects/visual_flocking/__init__.py b/abm/projects/visual_flocking/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/abm/projects/visual_flocking/vf_agent/__init__.py b/abm/projects/visual_flocking/vf_agent/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/abm/projects/visual_flocking/vf_agent/vf_agent.py b/abm/projects/visual_flocking/vf_agent/vf_agent.py new file mode 100644 index 00000000..f74b80ee --- /dev/null +++ b/abm/projects/visual_flocking/vf_agent/vf_agent.py @@ -0,0 +1,338 @@ +import matplotlib +import numpy as np +import pygame + +from abm.projects.visual_flocking.vf_agent import vf_supcalc +from abm.agent.agent import Agent +from abm.projects.visual_flocking.vf_contrib import vf_params +from abm.contrib import colors +import importlib + + +class VFAgent(Agent): + def __init__(self, **kwargs): + super().__init__(**kwargs) + importlib.reload(vf_supcalc) + importlib.reload(vf_params) + + # creating agent status + # flocking or emergency + self.agent_state = "flocking" + + # introducing agent-specific movement params + # if they are None, they will be read from vf_contrib.vf_params + self.ALP0 = None + self.BET0 = None + self.V0 = None + + # lines to follow + self.lines = [] + self.sensor_distance = 20 # sensor distance to follow lines + self.sensor_size = 9 + self.line_map = np.zeros((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + + # getting additional information from supplementary calculation functions + self.verbose_supcalc = False + # boundary conditions + # infinite or walls + self.boundary_cond = vf_params.BOUNDARY + self.limit_movement = vf_params.LIMIT_MOVEMENT + self.max_vel = vf_params.MAX_VEL + self.max_th = vf_params.MAX_TH + + # preparing phi values for algorithm according to FOV + self.PHI = np.arange(-np.pi, np.pi, (2*np.pi)/self.v_field_res) + + # social information: visual field projections + self.soc_v_field_proj = np.zeros(self.v_field_res) + + # saving collision times + self.coll_times = [] + + def update(self, agents): + """ + main update method of the agent. This method is called in every timestep + to calculate the new state/position + of the agent and visualize it in the environment + :param agents: a list of all obstacle/agents coordinates as (X, Y) in + the environment. These are not necessarily socially relevant, i.e. all + agents. + """ + # update agent information + self.update_social_info(agents) + + # update agent's state + self.update_state() + + # perform agent's action i.e. exploration, taxis, relocation or flocking + self.perform_action() + + # boundary conditions if applicable + if self.boundary_cond == "walls": + self.reflect_from_walls() + elif self.boundary_cond == "infinite": + self.teleport_infinite_arena() + + # updating agent visualization + self.draw_update() + + + def reflect_from_walls(self): + """reflecting agent from environment boundaries according to a desired x, y coordinate. If this is over any + boundaries of the environment, the agents position and orientation will be changed such that the agent is + reflected from these boundaries.""" + if len(self.lines) == 0: + super().reflect_from_walls() + else: + # Boundary conditions according to center of agent (simple) + x = self.position[0] + self.radius + y = self.position[1] + self.radius + + # Reflection from left wall + if x < self.boundaries_x[0]: + self.position[0] = self.boundaries_x[0] - self.radius + + if np.pi / 2 <= self.orientation < np.pi: + self.orientation -= np.pi + elif np.pi <= self.orientation <= 3 * np.pi / 2: + self.orientation += np.pi + self.prove_orientation() # bounding orientation into 0 and 2pi + + # Reflection from right wall + if x > self.boundaries_x[1]: + + self.position[0] = self.boundaries_x[1] - self.radius - 1 + + if 3 * np.pi / 2 <= self.orientation < 2 * np.pi: + self.orientation -= np.pi + elif 0 <= self.orientation <= np.pi / 2: + self.orientation += np.pi + self.prove_orientation() # bounding orientation into 0 and 2pi + + # Reflection from upper wall + if y < self.boundaries_y[0]: + self.position[1] = self.boundaries_y[0] - self.radius + + if np.pi / 2 <= self.orientation <= np.pi: + self.orientation += np.pi + elif 0 <= self.orientation < np.pi / 2: + self.orientation -= np.pi + self.prove_orientation() # bounding orientation into 0 and 2pi + + # Reflection from lower wall + if y > self.boundaries_y[1]: + self.position[1] = self.boundaries_y[1] - self.radius - 1 + if 3 * np.pi / 2 <= self.orientation <= 2 * np.pi: + self.orientation += np.pi + elif np.pi <= self.orientation < 3 * np.pi / 2: + self.orientation -= np.pi + self.prove_orientation() # bounding orientation into 0 and 2pi + + + def draw_update(self): + """ + updating the outlook of the agent according to position and orientation + """ + # update position + self.rect.x = self.position[0] + self.rect.y = self.position[1] + + # change agent color according to mode + self.change_color() + + # update surface according to new orientation + # creating visualization surface for agent as a filled circle + self.image = pygame.Surface([self.radius * 2, self.radius * 2]) + self.image.fill(colors.BACKGROUND) + self.image.set_colorkey(colors.BACKGROUND) + if self.is_moved_with_cursor: + pygame.gfxdraw.filled_circle( + self.image, + self.radius, + self.radius, + self.radius, + self.selected_color + ) + pygame.gfxdraw.aacircle(self.image, + self.radius, + self.radius, + self.radius, + colors.BACKGROUND) + else: + pygame.gfxdraw.filled_circle( + self.image, + self.radius, + self.radius, + self.radius-1, + self.color[0:-1] + ) + pygame.gfxdraw.aacircle(self.image, + self.radius, + self.radius, + self.radius-1, + colors.BLACK) + + # pygame.draw.circle( + # self.image, self.color, (self.radius, self.radius), self.radius + # ) + + # showing agent orientation with a line towards agent orientation + new_white = (255, 255, 254) + pygame.draw.line(self.image, new_white, (self.radius, self.radius), + ((1 + np.cos(self.orientation)) * self.radius, (1 - np.sin(self.orientation)) * self.radius), + 3) + self.image = pygame.transform.smoothscale(self.image, [self.radius * 2, self.radius * 2]) + self.mask = pygame.mask.from_surface(self.image) + + + def teleport_infinite_arena(self): + """In case the boundary conditions are infinite (as in now reflection from walls is requested) the + agents are teleported on a torus when reaching walls.""" + + # Boundary conditions according to center of agent (simple) + x = self.position[0] + self.radius + y = self.position[1] + self.radius + + if x < self.boundaries_x[0]: + self.position[0] = self.boundaries_x[1] - self.radius + elif x > self.boundaries_x[1]: + self.position[0] = self.boundaries_x[0] + self.radius + + if y < self.boundaries_y[0]: + self.position[1] = self.boundaries_y[1] - self.radius + elif y > self.boundaries_y[1]: + self.position[1] = self.boundaries_y[0] + self.radius + + def update_social_info(self, agents): + # calculate socially relevant projection field (e.g. according to + # signalling agents) + self.soc_v_field = self.calc_soc_v_proj(agents) + try: + dphi2 = vf_supcalc.follow_lines_local(self.position, self.radius, self.orientation, self.line_map, self.velocity, + sensor_radius=self.sensor_size, sensor_distance=self.sensor_distance) + except: + print("Error during line following!") + + def calc_soc_v_proj(self, agents): + """ + :param agents: all agents + """ + agents_of_interest = [ag for ag in agents if ag.id != self.id] + visual_field = vf_supcalc.projection_field( + fov=self.FOV, + v_field_resolution=self.v_field_res, + position=np.array(self.position), + radius=self.radius, + orientation=self.orientation, + object_positions=[np.array(ag.position) for ag in agents_of_interest], + object_sizes=[ag.radius for ag in agents_of_interest], + boundary_cond=self.boundary_cond, + arena_width=self.WIDTH, + arena_height=self.HEIGHT, + ag_id=self.id) + # sum of all agents projections at each point in visual field + svfield = visual_field.sum(axis=0) + # binarize v field + svfield[svfield > 0] = 1 + return svfield + + def update_state(self): + # update agent state based on the decision-making process + self.agent_state = "flocking" + + def update_linemap(self): + """Updating background line map to follow""" + subsurface = pygame.Surface((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + subsurface.fill(colors.BACKGROUND) + subsurface.set_colorkey(colors.WHITE) + subsurface.set_alpha(255) + for line in self.lines: + for pi in range(1, len(line)): + point1 = line[pi-1] + point2 = line[pi] + color = colors.BLACK + pygame.draw.line(subsurface, color, point1, point2, 1) + from scipy.ndimage.filters import gaussian_filter + line_map = pygame.surfarray.array3d(subsurface) + line_map = gaussian_filter(line_map, sigma=1) + line_map[line_map<253] = 0 + line_map = gaussian_filter(line_map, sigma=1) + self.line_map = 1 - line_map.swapaxes(0, 1)[:, :, 0]/255 + + + def perform_action(self): + # we freeze agents when we move them + if not self.is_moved_with_cursor: + # perform the agent's action according to the current state + if self.agent_state == "flocking": + if len(self.PHI) == len(self.soc_v_field): + if not self.verbose_supcalc: + dv, dphi = vf_supcalc.VSWRM_flocking_state_variables(self.velocity, self.PHI, np.flip(self.soc_v_field), + vf_params, verbose=self.verbose_supcalc, + ALP0=self.ALP0, BET0=self.BET0, V0=self.V0) + if len(self.lines) != 0: + dphi2 = vf_supcalc.follow_lines_local(self.position, self.radius, self.orientation, self.line_map, self.velocity, + sensor_radius=self.sensor_size, sensor_distance=self.sensor_distance) + dphi = dphi2 + else: + self.dv, self.dphi, self.ablob, self.aedge, self.bblob, self.bedge = vf_supcalc.VSWRM_flocking_state_variables(self.velocity, self.PHI, np.flip(self.soc_v_field), + vf_params, verbose=self.verbose_supcalc, + ALP0=self.ALP0, BET0=self.BET0, V0=self.V0) + dv, dphi = self.dv, self.dphi + else: + dv, dphi = 0, 0 + print("Warning: temporary skip calculation due to mismatch in phi and vfield resolutions!") + self.update_agent_position(dphi, dv) + elif self.agent_state == "emergency": + pass + + def update_agent_position(self, theta, vel): + # updating agent's state variables according to calculated vel and + # theta + # maximum turning angle per timestep + if self.limit_movement: + theta = self.prove_turning(theta, theta_lim=self.max_th) + self.orientation += theta + self.prove_orientation() # bounding orientation into 0 and 2pi + + self.velocity += vel + if self.limit_movement: + self.prove_velocity(velocity_limit=self.max_vel) # possibly bounding velocity of agent + + # new agent's position + new_pos = ( + self.position[0] + self.velocity * np.cos(self.orientation), + self.position[1] - self.velocity * np.sin(self.orientation) + ) + # update the agent's position with constraints (reflection from the + # walls) or with the new position + self.position = list(new_pos) + + def prove_turning(self, theta, theta_lim=0.2): + t_sign = np.sign(theta) + if t_sign == 0: + t_sign = +1 + if np.abs(theta) > theta_lim: + theta = theta_lim * t_sign + return theta + + + def prove_velocity(self, velocity_limit=1): + """Restricting the absolute velocity of the agent""" + vel_sign = np.sign(self.velocity) + if vel_sign == 0: + vel_sign = +1 + if np.abs(self.velocity) > velocity_limit: + self.velocity = velocity_limit * vel_sign + + + def change_color(self): + """ + Changing color of agent according to the behavioral mode the agent is + currently in. + """ + cmap = matplotlib.cm.get_cmap('jet') + rgba = np.array(cmap(self.orientation / (2 * np.pi))) + # rescaling color for pygame + rgba[0:3] *= 255 + self.color = rgba diff --git a/abm/projects/visual_flocking/vf_agent/vf_supcalc.py b/abm/projects/visual_flocking/vf_agent/vf_supcalc.py new file mode 100644 index 00000000..54d7aa52 --- /dev/null +++ b/abm/projects/visual_flocking/vf_agent/vf_supcalc.py @@ -0,0 +1,329 @@ +from math import atan2 + +import numpy as np + +from abm.agent.supcalc import angle_between, find_nearest + +from scipy import integrate + +def distance_infinite(p1, p2, L=500, dim=2): + """ Returns the distance vector of two position vectors x,y + by tanking periodic boundary conditions into account. + Input parameters: L - system size, dim - no. of dimension + """ + distvec = p2 - p1 + distvec_periodic = np.copy(distvec) + distvec_periodic[distvec < -0.5*L] += L + distvec_periodic[distvec > 0.5*L] -= L + return distvec_periodic + +def projection_field(fov, v_field_resolution, position, radius, + orientation, object_positions, object_sizes=None, + boundary_cond="walls", arena_width=None, arena_height=None, vision_range=None, ag_id=0): + """ + Calculating visual projection field for the agent given the visible + obstacles in the environment + obstacle sprites to generate projection field + :param fov: tuple of number with borders of fov such as (-np.pi, np.pi) + :param v_field_resolution: visual field resolution in pixels + :param position: np.xarray of agent's position + :param radius: radius of the agent + :param orientation: orientation angle between 0 and 2pi + :param max_proj_size: maximum projection size to include in the visual proj. field + :param boundary_cond: boundary condition either infinite or walls. If walls projection field + is calculated according to euclidean coordinates in 2D space, otherwise on a torus. + :return: projection field np.xarray with shape (n objects, field resolution) + """ + # initializing visual field and relative angles + v_field = np.zeros((len(object_positions), v_field_resolution)) + phis = np.linspace(-np.pi, np.pi, v_field_resolution) + + # center point + agents_center = position + radius + + # point on agent's edge circle according to it's orientation + agent_edge = position + np.array( + [1 + np.cos(orientation), 1 - np.sin(orientation)]) * radius + + # vector between center and edge according to orientation + v1 = agent_edge - agents_center + + # 1. Calculating closed angle between object and agent according to the + # position of the object. + # 2. Calculating visual projection size according to visual angle on + # the agent's retina according to distance between agent and object + for i, obj_position in enumerate(object_positions): + # continue if the object and the agent position completely coincide + if obj_position[0] == position[0] and obj_position[1] == position[1]: + continue + + # center of obstacle (as it must be another agent) + if object_sizes is None: + object_center = obj_position + radius + else: + object_center = obj_position + object_sizes[i] + + # vector between agent center and object center + v2 = object_center - agents_center + + # in case torus, positions might change + if boundary_cond == "infinite": + if np.abs(v2[0]) > arena_width/2: + if agents_center[0] < object_center[0]: + object_center[0] -= arena_width + elif agents_center[0] > object_center[0]: + object_center[0] += arena_width + if np.abs(v2[1]) > arena_height/2: + if agents_center[1] < object_center[1]: + object_center[1] -= arena_height + elif agents_center[1] > object_center[1]: + object_center[1] += arena_height + + # recalculating v2 after teleporting on torus + v2 = object_center - agents_center + + # calculating closed angle between v1 and v2 + closed_angle = calculate_closed_angle(v1, v2) + + distance = np.linalg.norm(object_center - agents_center) + + # limiting vision range if requested + if vision_range is not None: + if distance > vision_range: + continue + + # calculating the visual angle from focal agent to target + if object_sizes is None: + vis_angle = 2 * np.arctan(radius / (1 * distance)) + else: + vis_angle = 2 * np.arctan(object_sizes[i] / (1 * distance)) + + # finding where in the retina the projection belongs to + phi_target = find_nearest(phis, closed_angle) + + # transforming fov boundaries to pixel boundaries + fov_px = [find_nearest(phis, fov[0]), find_nearest(phis, fov[1])] + + # # if target is visible we save its projection into the VPF + # # source data + # if fov[0] <= closed_angle <= fov[1]: + # the projection size is proportional to the visual angle. + # If the projection is maximal (i.e. taking each pixel of the + # retina) the angle is 2pi from this we just calculate the + # projection size using a single proportion + proj_size = (vis_angle / (2 * np.pi)) * v_field_resolution + + proj_start = int(phi_target - np.floor(proj_size / 2)) + proj_end = int(phi_target + np.floor(proj_size / 2)) + + if fov_px[0] < proj_start < fov_px[1] or fov_px[0] < proj_end < fov_px[1]: + + # circular boundaries to the VPF as there is 360 degree vision + if proj_start < 0: + v_field[i, v_field_resolution + proj_start:v_field_resolution] = 1 + proj_start = 0 + if proj_end >= v_field_resolution: + v_field[i, 0:proj_end - (v_field_resolution - 1)] = 1 + proj_end = v_field_resolution + + v_field[i, proj_start:proj_end] = 1 + + + # post_processing and limiting FOV + # flip field data along second dimension + v_field_post = np.flip(v_field, axis=1) + + # v_field_post[:, phis < fov[0]] = 0 + # v_field_post[:, phis > fov[1]] = 0 + return v_field_post + + + +def calculate_closed_angle(v1, v2): + """ + Calculating closed angle between two vectors v1 and v2. Rotated with the orientation of the agent as it is relative. + :param v1: vector 1; np.xarray + :param v2: vector 2; np.xarray + :return: closed angle between v1 and v2 + """ + closed_angle = angle_between(v1, v2) + closed_angle = (closed_angle % (2 * np.pi)) + # at this point closed angle between 0 and 2pi, but we need it between -pi and pi + # we also need to take our orientation convention into consideration to + # recalculate theta=0 is pointing to the right + if 0 <= closed_angle <= np.pi: + closed_angle = -closed_angle + else: + closed_angle = 2 * np.pi - closed_angle + return closed_angle + +# Functions needed for VSWRM functionality +def VSWRM_flocking_state_variables(vel_now, Phi, V_now, vf_params, t_now=None, V_prev=None, t_prev=None, verbose=False, + ALP0=None, BET0=None, V0=None): + """Calculating state variables of a given agent according to the main algorithm as in + https://advances.sciencemag.org/content/6/6/eaay0792. + Args: + vel_now: current speed of the agent + V_now: current binary visual projection field array + Phi: linspace numpy array of visual field axis + vf_params: parameters of the flocking algorithm as in main article, e.g.: V0, GAM, ALP0, etc. + t_now: current time + V_prev: previous binary visual projection field array + t_prev: previous time + ALP0: overwriting alpha0 parameter for heterogenity + BET: overwriting bet0 for heterogeneity + V0: overwriting self-propelled speed for heterogenity + Returns: + dvel: temporal change in agent velocity + dpsi: temporal change in agent heading angle + + """ + # # Deriving over t is omitted in simplest case + # if V_prev is not None and t_prev is not None and t_now is not None: + # dt = t_now - t_prev + # logger.debug('Movement calculation called with NONE as time-related parameters.') + # joined_V = np.vstack((V_prev, t_prev)) + # dt_V = dt_V_of(dt, joined_V) + # else: + # dt_V = np.zeros(len(Phi)) + + # Overwriting default homogeneous variables if ALP0, BET0 or V0 is provided + if ALP0 is None: + ALP0 = vf_params.ALP0 + if BET0 is None: + BET0 = vf_params.BET0 + if V0 is None: + V0 = vf_params.V0 + + # Using only zeros for temporal derivative in simplest case to keep formulation general + dt_V = np.zeros(len(Phi)) + + # Deriving over Phi + dPhi_V = dPhi_V_of(Phi, V_now) + + # Calculating series expansion of functional G + G_vel = (-V_now + vf_params.ALP2 * dt_V) + + # Spikey parts shall be handled separately because of numerical integration + G_vel_spike = np.square(dPhi_V) + + G_psi = (-V_now + vf_params.BET2 * dt_V) + + # Spikey parts shall be handled separately because of numerical integration + G_psi_spike = np.square(dPhi_V) + + # Calculating change in velocity and heading direction + dPhi = Phi[-1] - Phi[-2] + FOV_rescaling_cos = 1 + FOV_rescaling_sin = 1 + + # print(f"alp0 : {vswrm.ALP0 * integrate.trapz(np.cos(FOV_rescaling_cos * Phi) * G_vel, Phi)}", ) + # print(f'alp1 : {vswrm.ALP0 * vswrm.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) * dPhi}') + + # dvel = vf_params.GAM * (vf_params.V0 - vel_now) + \ + # vf_params.ALP0 * integrate.trapz(np.cos(FOV_rescaling_cos * Phi) * G_vel, Phi) + \ + # vf_params.ALP0 * vf_params.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) * dPhi + # dpsi = vf_params.BET0 * integrate.trapz(np.sin(Phi) * G_psi, Phi) + \ + # vf_params.BET0 * vf_params.BET1 * np.sum(np.sin(FOV_rescaling_sin * Phi) * G_psi_spike) * dPhi + + + + if not verbose: + # without reacling edge information + dvel = vf_params.GAM * (V0 - vel_now) + \ + ALP0 * integrate.trapz(np.cos(Phi) * G_vel, Phi) + \ + ALP0 * vf_params.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) + + dpsi = BET0 * integrate.trapz(np.sin(Phi) * G_psi, Phi) + \ + BET0 * vf_params.BET1 * np.sum(np.sin(Phi) * G_psi_spike) + + return dvel, dpsi + else: + alpha_blob = ALP0 * integrate.trapz(np.cos(Phi) * G_vel, Phi) + alpha_edge = ALP0 * vf_params.ALP1 * np.sum(np.cos(Phi) * G_vel_spike) + + beta_blob = BET0 * integrate.trapz(np.sin(Phi) * G_psi, Phi) + beta_edge = BET0 * vf_params.BET1 * np.sum(np.sin(Phi) * G_psi_spike) + # without reacling edge information + dvel = vf_params.GAM * (V0 - vel_now) + \ + alpha_blob + \ + alpha_edge + + dpsi = beta_blob + \ + beta_edge + return dvel, dpsi, alpha_blob, alpha_edge, beta_blob, beta_edge + + +def dPhi_V_of(Phi, V): + """Calculating derivative of VPF according to Phi visual angle array at a given timepoint t + Args: + Phi: linspace numpy array of visual field axis + V: binary visual projection field array + Returns: + dPhi_V: derivative array of V w.r.t Phi + """ + # circular padding for edge cases + padV = np.pad(V, (1, 1), 'wrap') + dPhi_V_raw = np.diff(padV) + + # we want to include non-zero value if it is on the edge + if dPhi_V_raw[0] > 0 and dPhi_V_raw[-1] > 0: + dPhi_V_raw = dPhi_V_raw[0:-1] + + else: + dPhi_V_raw = dPhi_V_raw[1:, ...] + + dPhi_V = dPhi_V_raw #/ (Phi[-1] - Phi[-2]) + return dPhi_V + +def distance_coords(x1, y1, x2, y2, vectorized=False): + """Distance between 2 points in 2D space calculated from point coordinates. + if vectorized is True, we use multidimensional (i.e. vectorized) form of distance + calculation that preserved original dimensions of coordinate arrays in the dimensions of the output and the output + will contain pairwise distance measures according to coordinate matrices.""" + c1 = np.array([x1, y1]) + c2 = np.array([x2, y2]) + if not vectorized: + distance = np.linalg.norm(c2 - c1) + else: + distance = np.linalg.norm(c2 - c1, axis=0) + return distance + + +def follow_lines_local(agposition, agradius, agorientation, linemap, agvel, sensor_radius=10, sensor_distance=5): + """Following line with 2 sensors""" + sensor1_pos = [agposition[1] + agradius - sensor_distance + ( + 1 + np.sin(agorientation + (3*np.pi / 4))) * sensor_distance, + agposition[0] + agradius - sensor_distance + ( + 1 - np.cos(agorientation + (3*np.pi / 4))) * sensor_distance] + sensor2_pos = [agposition[1] + agradius - sensor_distance + ( + 1 + np.sin(agorientation - (3*np.pi / 4))) * sensor_distance, + agposition[0] + agradius - sensor_distance + ( + 1 - np.cos(agorientation - (3*np.pi / 4))) * sensor_distance] + # superline = [] + # for line in lines: + # superline.extend(line) + + # line = superline + # points_s1_range = np.array([point for point in line if distance_coords(sensor1_pos[1], sensor1_pos[0], point[1], point[0]) s2: + return ori_change + elif s1 < s2: + return ori_change + else: + if s1 != 0: + return 0.01 + else: + return 0 \ No newline at end of file diff --git a/abm/projects/visual_flocking/vf_contrib/__init__.py b/abm/projects/visual_flocking/vf_contrib/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/abm/projects/visual_flocking/vf_contrib/vf_params.py b/abm/projects/visual_flocking/vf_contrib/vf_params.py new file mode 100644 index 00000000..eab4e623 --- /dev/null +++ b/abm/projects/visual_flocking/vf_contrib/vf_params.py @@ -0,0 +1,23 @@ +# loading env variables from dotenv file +from dotenv import dotenv_values +import os + +EXP_NAME = os.getenv("EXPERIMENT_NAME", "") + +root_abm_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname( + os.path.dirname(os.path.realpath(__file__)))))) +env_path = os.path.join(root_abm_dir, f"{EXP_NAME}.env") +envconf = dotenv_values(env_path) + +GAM = float(envconf.get("VF_GAMMA", 0.1)) +V0 = float(envconf.get("VF_V0", 1)) +ALP0 = float(envconf.get("VF_ALP0", 1)) +ALP1 = float(envconf.get("VF_ALP1", 0.09)) +ALP2 = float(envconf.get("VF_ALP2", 0)) +BET0 = float(envconf.get("VF_BET0", 1)) +BET1 = float(envconf.get("VF_BET1", 0.09)) +BET2 = float(envconf.get("VF_BET2", 0)) +BOUNDARY = envconf.get("BOUNDARY", "walls") +LIMIT_MOVEMENT = bool(float(envconf.get("VF_LIMIT_MOVEMENT", "0"))) +MAX_VEL = float(envconf.get("VF_MAX_VEL", "3")) +MAX_TH = float(envconf.get("VF_MAX_TH", "0.1")) diff --git a/abm/projects/visual_flocking/vf_contrib/vf_playgroundtool.py b/abm/projects/visual_flocking/vf_contrib/vf_playgroundtool.py new file mode 100644 index 00000000..09b75e4a --- /dev/null +++ b/abm/projects/visual_flocking/vf_contrib/vf_playgroundtool.py @@ -0,0 +1,90 @@ +from abm.contrib import playgroundtool as pgt + + +def setup_visflock_playground(): + playground_tool = pgt + # default_params + # update default parameters + playground_tool.default_params["framerate"] = 60 + playground_tool.default_params["N_resc"] = 0 + playground_tool.default_params["patch_radius"] = 0 + playground_tool.default_params["width"] = 500 + playground_tool.default_params["height"] = 500 + # # new default parameters + playground_tool.default_params["collide_agents"] = False + playground_tool.default_params["agent_radius"] = 8 + playground_tool.default_params["N"] = 10 + playground_tool.default_params["v_field_res"] = 800 + #playground_tool.default_params["alpha_0"] = 4 + # playground_tool.default_params["phototaxis_theta_step"] = 0.2 + # playground_tool.default_params["detection_range"] = 120 + # playground_tool.default_params["resource_meter_multiplier"] = 1 + # playground_tool.default_params["des_velocity_res"] = 1.5 + # playground_tool.default_params["res_theta_abs"] = 0.2 + # playground_tool.default_params["signalling_cost"] = 0.2 + # playground_tool.default_params["probability_of_starting_signaling"] = 0.5 + # playground_tool.default_params["agent_signaling_rand_event_update"] = 10 + + # # update def_params_to_env_vars + # playground_tool.def_params_to_env_vars[ + # "phototaxis_theta_step"] = "PHOTOTAX_THETA_FAC" + # playground_tool.def_params_to_env_vars[ + # "detection_range"] = "DETECTION_RANGE" + # playground_tool.def_params_to_env_vars[ + # "resource_meter_multiplier"] = "METER_TO_RES_MULTI" + # playground_tool.def_params_to_env_vars[ + # "signalling_cost"] = "SIGNALLING_COST" + # playground_tool.def_params_to_env_vars[ + # "probability_of_starting_signaling"] = "SIGNALLING_PROB" + # playground_tool.def_params_to_env_vars["des_velocity_res"] = "RES_VEL" + # playground_tool.def_params_to_env_vars["res_theta_abs"] = "RES_THETA" + # playground_tool.def_params_to_env_vars["agent_signaling_rand_event_update"] = "SIGNAL_PROB_UPDATE_FREQ" + + # # update help_messages + playground_tool.help_messages["V.Resol."] = ''' + Resolution of Visual Projection Field [px]: + + Desired resolution of visual projection field, i.e. for how many pixel + is the retina of the agnt is divided. Ultimately this affects the visual range + of agents and their movement precision in very large distances and small resolutions. + ''' + playground_tool.help_messages["R ag."] = ''' + Agent radius [px]: + + Radius of the agents in pixels. + ''' + + playground_tool.help_messages["Alpha0"] = ''' + Alpha 0 [AU]: + + Sensitivity of agents according to the main algorithm as in + https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7002123/ + ''' + + playground_tool.help_messages["Beta0"] = ''' + Beta 0 [AU]: + + Maneauverability of agents according to the main algorithm as in + https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7002123/ + ''' + + playground_tool.help_messages["Alpha1Beta1"] = ''' + Alpha 1 = Beta 1 [AU]: + + Higher level parameters of the flocking algorithm controlling front-back, and left-right + equilibrium distances respectively. + ''' + + playground_tool.help_messages["V0"] = ''' + V0 [px/ts]: + + Preferred/self-propelled speed of the agents in pixels/timestep + ''' + + playground_tool.help_messages["V0"] = ''' + Path length [ts]: + + Path length to save and visualize + ''' + + return playground_tool diff --git a/abm/projects/visual_flocking/vf_simulation/__init__.py b/abm/projects/visual_flocking/vf_simulation/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/abm/projects/visual_flocking/vf_simulation/vf_isims.py b/abm/projects/visual_flocking/vf_simulation/vf_isims.py new file mode 100644 index 00000000..9a2c5ec2 --- /dev/null +++ b/abm/projects/visual_flocking/vf_simulation/vf_isims.py @@ -0,0 +1,444 @@ +import os + +import numpy as np +import pygame +import pygame_widgets +from pygame_widgets.button import Button +from pygame_widgets.slider import Slider +from pygame_widgets.textbox import TextBox + +from abm.contrib import colors +from abm.monitoring.ifdb import pad_to_n_digits +from abm.projects.visual_flocking.vf_simulation.vf_sims import VFSimulation +from abm.simulation.isims import PlaygroundSimulation +from abm.projects.visual_flocking.vf_contrib import vf_params + + + +class VFPlaygroundSimulation(PlaygroundSimulation, VFSimulation): + """ + SEE: https://docs.python.org/3/tutorial/classes.html#multiple-inheritance + """ + + def __init__(self): + super().__init__() + # Hiding unused buttons and sliders + self.fix_SUM_res_button.hide() + self.visual_exclusion_button.hide() + self.IFDB_button.hide() + + + # Replacing NRES slider to Visual Resolution Slider + self.NRES_help.hide() + self.NRES_textbox.hide() + self.NRES_slider.hide() + x = self.NRES_slider.getX() + y = self.NRES_slider.getY() + w = self.NRES_slider.getWidth() + h = self.NRES_slider.getHeight() + self.VISRES_slider = Slider(self.screen, x, y, w, h, min=20, max=3000, step=10, initial=self.v_field_res) + x = self.NRES_textbox.getX() + y = self.NRES_textbox.getY() + w = self.NRES_textbox.getWidth() + h = self.NRES_textbox.getHeight() + self.VISRES_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.NRES_help.getX() + y = self.NRES_help.getY() + w = self.NRES_help.getWidth() + h = self.NRES_help.getHeight() + self.VISRES_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.VISRES_help.onClick = lambda: self.show_help('V.Resol.', self.VISRES_help) + self.VISRES_help.onRelease = lambda: self.unshow_help(self.VISRES_help) + + self.help_buttons.append(self.VISRES_help) + self.sliders.append(self.VISRES_slider) + self.slider_texts.append(self.VISRES_textbox) + + + # Replacing Res radius slider with agent radius slider + self.RESradius_slider.hide() + self.RESradius_textbox.hide() + self.RES_help.hide() + x = self.RESradius_slider.getX() + y = self.RESradius_slider.getY() + w = self.RESradius_slider.getWidth() + h = self.RESradius_slider.getHeight() + self.ARAD_slider = Slider(self.screen, x, y, w, h, min=1, max=30, step=1, initial=self.agent_radii) + x = self.RESradius_textbox.getX() + y = self.RESradius_textbox.getY() + w = self.RESradius_textbox.getWidth() + h = self.RESradius_textbox.getHeight() + self.ARAD_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.RES_help.getX() + y = self.RES_help.getY() + w = self.RES_help.getWidth() + h = self.RES_help.getHeight() + self.ARAD_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.ARAD_help.onClick = lambda: self.show_help('R ag.', self.VISRES_help) + self.ARAD_help.onRelease = lambda: self.unshow_help(self.VISRES_help) + + self.help_buttons.append(self.ARAD_help) + self.sliders.append(self.ARAD_slider) + self.slider_texts.append(self.ARAD_textbox) + + + # Replacing Ew slider with Alpha0 + self.Epsw_slider.hide() + self.Epsw_textbox.hide() + self.Epsw_help.hide() + x = self.Epsw_slider.getX() + y = self.Epsw_slider.getY() + w = self.Epsw_slider.getWidth() + h = self.Epsw_slider.getHeight() + self.ALP0_slider = Slider(self.screen, x, y, w, h, min=0, max=5, step=0.01, initial=vf_params.ALP0) + x = self.Epsw_textbox.getX() + y = self.Epsw_textbox.getY() + w = self.Epsw_textbox.getWidth() + h = self.Epsw_textbox.getHeight() + self.ALP0_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.Epsw_help.getX() + y = self.Epsw_help.getY() + w = self.Epsw_help.getWidth() + h = self.Epsw_help.getHeight() + self.ALP0_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.ALP0_help.onClick = lambda: self.show_help('Alpha0', self.ALP0_help) + self.ALP0_help.onRelease = lambda: self.unshow_help(self.ALP0_help) + + self.help_buttons.append(self.ALP0_help) + self.sliders.append(self.ALP0_slider) + self.slider_texts.append(self.ALP0_textbox) + + # Replacing Eu slider with Beta0 + self.Epsu_slider.hide() + self.Epsu_textbox.hide() + self.Epsu_help.hide() + x = self.Epsu_slider.getX() + y = self.Epsu_slider.getY() + w = self.Epsu_slider.getWidth() + h = self.Epsu_slider.getHeight() + self.BET0_slider = Slider(self.screen, x, y, w, h, min=0, max=5, step=0.01, initial=vf_params.BET0) + x = self.Epsu_textbox.getX() + y = self.Epsu_textbox.getY() + w = self.Epsu_textbox.getWidth() + h = self.Epsu_textbox.getHeight() + self.BET0_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.Epsu_help.getX() + y = self.Epsu_help.getY() + w = self.Epsu_help.getWidth() + h = self.Epsu_help.getHeight() + self.BET0_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.BET0_help.onClick = lambda: self.show_help('Beta0', self.BET0_help) + self.BET0_help.onRelease = lambda: self.unshow_help(self.BET0_help) + + self.help_buttons.append(self.BET0_help) + self.sliders.append(self.BET0_slider) + self.slider_texts.append(self.BET0_textbox) + + # Replacing Suw slider with A1B1 + self.SWU_slider.hide() + self.SWU_textbox.hide() + self.SWU_help.hide() + x = self.SWU_slider.getX() + y = self.SWU_slider.getY() + w = self.SWU_slider.getWidth() + h = self.SWU_slider.getHeight() + self.ALP1BET1_slider = Slider(self.screen, x, y, w, h, min=0, max=0.1, step=0.001, initial=vf_params.ALP1) + x = self.SWU_textbox.getX() + y = self.SWU_textbox.getY() + w = self.SWU_textbox.getWidth() + h = self.SWU_textbox.getHeight() + self.ALP1BET1_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.SWU_help.getX() + y = self.SWU_help.getY() + w = self.SWU_help.getWidth() + h = self.SWU_help.getHeight() + self.ALP1BET1_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.ALP1BET1_help.onClick = lambda: self.show_help('Alpha1Beta1', self.ALP1BET1_help) + self.ALP1BET1_help.onRelease = lambda: self.unshow_help(self.ALP1BET1_help) + + self.help_buttons.append(self.ALP1BET1_help) + self.sliders.append(self.ALP1BET1_slider) + self.slider_texts.append(self.ALP1BET1_textbox) + + # Replacing Suw slider with V0 + self.SUW_slider.hide() + self.SUW_textbox.hide() + self.SUW_help.hide() + x = self.SUW_slider.getX() + y = self.SUW_slider.getY() + w = self.SUW_slider.getWidth() + h = self.SUW_slider.getHeight() + self.V0_slider = Slider(self.screen, x, y, w, h, min=0, max=5, step=0.01, initial=vf_params.V0) + x = self.SUW_textbox.getX() + y = self.SUW_textbox.getY() + w = self.SUW_textbox.getWidth() + h = self.SUW_textbox.getHeight() + self.V0_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.SUW_help.getX() + y = self.SUW_help.getY() + w = self.SUW_help.getWidth() + h = self.SUW_help.getHeight() + self.V0_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.V0_help.onClick = lambda: self.show_help('V0', self.V0_help) + self.V0_help.onRelease = lambda: self.unshow_help(self.V0_help) + + self.help_buttons.append(self.V0_help) + self.sliders.append(self.V0_slider) + self.slider_texts.append(self.V0_textbox) + + # Replacing SUMR slider with Path Length + self.SUMR_slider.hide() + self.SUMR_textbox.hide() + self.SUMR_help.hide() + x = self.SUMR_slider.getX() + y = self.SUMR_slider.getY() + w = self.SUMR_slider.getWidth() + h = self.SUMR_slider.getHeight() + self.memory_length = 0 + self.show_path_history = False + self.PLEN_slider = Slider(self.screen, x, y, w, h, min=0, max=100, step=1, initial=self.memory_length) + x = self.SUMR_textbox.getX() + y = self.SUMR_textbox.getY() + w = self.SUMR_textbox.getWidth() + h = self.SUMR_textbox.getHeight() + self.PLEN_textbox = TextBox(self.screen, x, y, w, h, fontSize=self.textbox_height - 2, borderThickness=1) + x = self.SUMR_help.getX() + y = self.SUMR_help.getY() + w = self.SUMR_help.getWidth() + h = self.SUMR_help.getHeight() + self.PLEN_help = Button(self.screen, x, y, w, h, text='?', fontSize=self.help_height - 2, + inactiveColour=colors.GREY, borderThickness=1, ) + self.PLEN_help.onClick = lambda: self.show_help('PLEN', self.PLEN_help) + self.PLEN_help.onRelease = lambda: self.unshow_help(self.PLEN_help) + + self.help_buttons.append(self.PLEN_help) + self.sliders.append(self.PLEN_slider) + self.slider_texts.append(self.PLEN_textbox) + + + def draw_frame(self, stats, stats_pos): + """ + Overwritten method of sims drawframe adding possibility to update + pygame widgets + """ + super().draw_frame(stats, stats_pos) + self.framerate_textbox.setText(f"Framerate: {self.framerate}") + self.N_textbox.setText(f"N: {self.N}") + self.VISRES_textbox.setText(f"V.Resol.: {self.v_field_res}") + self.FOV_textbox.setText(f"FOV: {int(self.fov_ratio * 100)}%") + self.ARAD_textbox.setText(f"R ag.: {self.agent_radii:.2f}") + self.ALP0_textbox.setText(f"A0: {vf_params.ALP0:.2f}") + self.BET0_textbox.setText(f"B0: {vf_params.BET0:.2f}") + self.ALP1BET1_textbox.setText(f"A1B1: {vf_params.ALP1:.2f}") + self.V0_textbox.setText(f"V0: {vf_params.V0:.2f}") + self.PLEN_textbox.setText(f"path len.: {self.memory_length:.2f}") + if self.SUM_res == 0: + self.update_SUMR() + self.SUMR_textbox.setText(f"SUM R: {self.SUM_res:.2f}") + for sl in self.sliders: + sl.draw() + for slt in self.slider_texts: + slt.draw() + for hb in self.help_buttons: + hb.draw() + for fb in self.function_buttons: + fb.draw() + if self.is_help_shown: + self.draw_help_message() + self.draw_global_stats() + if self.is_recording: + self.draw_record_circle() + if self.save_video: + # Showing the help message before the screen freezes + self.help_message = "\n\n\n Saving video, please wait..." + self.draw_help_message() + pygame.display.flip() + # Save video (freezes update process for a while) + self.saved_images_to_video() + self.save_video = False + # self.keep_agents_in_center() + #self.screen = pygame.transform.smoothscale(self.screen, (self.screen.get_width(), self.screen.get_height())) + + # def keep_agents_in_center(self): + # """Translating all agents with the center of mass, so they keep in the center of the arena""" + # posx_coords = np.array([agent.position[0] for agent in self.agents]) + # posy_coords = np.array([agent.position[1] for agent in self.agents]) + # COM = (np.mean(posx_coords), np.mean(posy_coords)) + # arena_center_x = int(self.vis_area_end_width/2) + # arena_center_y = int(self.vis_area_end_height / 2) + # dx = COM[0] - arena_center_x + # dy = COM[1] - arena_center_y + # for ai, agent in enumerate(list(self.agents)): + # agent.position[0] = agent.position[0] - dx + # agent.position[1] = agent.position[1] - dy + # if self.pos_memory is not None: + # self.pos_memory[ai, 0, -1] += dx + # self.pos_memory[ai, 1, -1] += dy + + + + def interact_with_event(self, events): + """Carry out functionality according to user's interaction""" + super().interact_with_event(events) + pygame_widgets.update(events) + + for event in events: + if event.type == pygame.KEYUP and event.key == pygame.K_l: + print("L UP, breaking line") + self.agent_for_line_id = list(set(self.agent_for_line_id)) + for agid, ag in enumerate(list(self.agents)): + if agid in self.agent_for_line_id: + print(f"Updating linemap for agent {agid}") + ag.lines.append(self.lines[-1]) + ag.update_linemap() + self.agent_for_line_id.remove(agid) + self.lines.append([]) + + self.framerate = self.framerate_slider.getValue() + self.N = self.N_slider.getValue() + self.N_resc = self.NRES_slider.getValue() + self.fov_ratio = self.FOV_slider.getValue() + if self.N != len(self.agents): + self.act_on_N_mismatch() + if self.fov_ratio != self.agent_fov[1] / np.pi: + self.update_agent_fovs() + if self.v_field_res != self.VISRES_slider.getValue(): + self.v_field_res = self.VISRES_slider.getValue() + self.update_agent_vfiled_res() + if self.agent_radii != self.ARAD_slider.getValue(): + self.agent_radii = self.ARAD_slider.getValue() + self.update_agent_radii() + if vf_params.ALP0 != self.ALP0_slider.getValue(): + vf_params.ALP0 = self.ALP0_slider.getValue() + if vf_params.BET0 != self.BET0_slider.getValue(): + vf_params.BET0 = self.BET0_slider.getValue() + if vf_params.ALP1 != self.ALP1BET1_slider.getValue(): + vf_params.ALP1 = self.ALP1BET1_slider.getValue() + vf_params.BET1 = self.ALP1BET1_slider.getValue() + if vf_params.V0 != self.V0_slider.getValue(): + vf_params.V0 = self.V0_slider.getValue() + if self.memory_length != self.PLEN_slider.getValue(): + self.memory_length = self.PLEN_slider.getValue() + if self.memory_length == 0: + self.show_path_history = False + else: + self.show_path_history = True + self.ori_memory = None + self.pos_memory = None + + if self.is_recording: + filename = f"{pad_to_n_digits(self.t, n=6)}.jpeg" + path = os.path.join(self.image_save_path, filename) + pygame.image.save(self.screen, path) + + + def update_agent_radii(self): + """Updating agent radius according to slider value""" + for agent in self.agents: + agent.radius = self.agent_radii + + def update_agent_vfiled_res(self): + """Updating visual field resolution of agents""" + for agent in self.agents: + agent.v_field_res = self.v_field_res + # preparing phi values for algorithm according to FOV + agent.PHI = np.arange(-np.pi, np.pi, (2 * np.pi) / self.v_field_res) + + # social information: visual field projections + agent.soc_v_field_proj = np.zeros(self.v_field_res) + + def act_on_NRES_mismatch(self): + """ + method is called if the requested amount of patches is not the same + as what the playground already has + """ + pass + + def act_on_N_mismatch(self): + """ + method is called if the requested amount of agents is not the same as + what the playground already has + """ + if self.N > len(self.agents): + diff = self.N - len(self.agents) + for i in range(diff): + ag_id = len(self.agents) + orient = np.random.uniform(0, 2 * np.pi) + dist = np.random.rand() * ( + self.HEIGHT / 2 - 2 * self.window_pad - + self.agent_radii) + x = np.cos(orient) * dist + self.WIDTH / 2 + y = np.sin(orient) * dist + self.HEIGHT / 2 + self.add_new_agent(ag_id, x, y, orient, with_proove=False) + self.update_agent_decision_params() + self.update_agent_fovs() + else: + while self.N < len(self.agents): + for i, ag in enumerate(self.agents): + if i == len(self.agents) - 1: + ag.kill() + if self.show_all_stats: + for ag in self.agents: + ag.show_stats = True + self.stats, self.stats_pos = self.create_vis_field_graph() + + def draw_visual_fields(self): + """ + Visualizing the range of vision for agents as opaque circles around the + agents + """ + r = 100 + + for agent in self.agents: + visual_field = agent.soc_v_field_proj + phis = np.linspace(agent.orientation - np.pi, agent.orientation + np.pi, + visual_field.shape[0]) + fov = agent.FOV + + # Center and radius of pie chart + cx = agent.position[0] + agent.radius + cy = agent.position[1] + agent.radius + + p_proj = [(cx + int(agent.radius * np.cos(phis[0])), + cy + int(agent.radius * - np.sin(phis[0])))] + for i, ang in enumerate(phis): + height = visual_field[i] + x = cx + int(r * np.cos(ang) * height) + y = cy + int(r * - np.sin(ang) * height) + p_proj.append((x, y)) + p_proj.append((cx + int(agent.radius * np.cos(phis[-1])), + cy + int(agent.radius * - np.sin(phis[-1])))) + + image = pygame.Surface( + [self.vis_area_end_width, self.vis_area_end_height]) + image.fill(colors.BACKGROUND) + image.set_colorkey(colors.BACKGROUND) + image.set_alpha(10) + + if 0 < fov[1] < np.pi: + p_fov = [(cx, cy)] + # Get points on arc + angles = [agent.orientation + fov[0], + agent.orientation + fov[1]] + step_size = (angles[1] - angles[0]) / 50 + angles_array = np.arange(angles[0], angles[1] + step_size, + step_size) + for n in angles_array: + x = cx + int(r * np.cos(n)) + y = cy + int(r * - np.sin(n)) + p_fov.append((x, y)) + p_fov.append((cx, cy)) + pygame.draw.polygon(image, colors.GREEN, p_fov) + elif fov[1] == np.pi: + cx, cy, r = agent.position[0] + agent.radius, agent.position[ + 1] + agent.radius, 100 + pygame.draw.circle(image, colors.GREEN, (cx, cy), r) + + pygame.draw.polygon(image, colors.BLACK, p_proj) + self.screen.blit(image, (0, 0)) diff --git a/abm/projects/visual_flocking/vf_simulation/vf_sims.py b/abm/projects/visual_flocking/vf_simulation/vf_sims.py new file mode 100644 index 00000000..c0bafbd3 --- /dev/null +++ b/abm/projects/visual_flocking/vf_simulation/vf_sims.py @@ -0,0 +1,397 @@ +from datetime import datetime + +import numpy as np +import pygame + +from abm.agent import supcalc +from abm.contrib import colors +from abm.monitoring import ifdb, env_saver +from abm.projects.visual_flocking.vf_agent.vf_agent import VFAgent +from abm.projects.visual_flocking.vf_contrib import vf_params +from abm.simulation.sims import Simulation +from matplotlib import cm as colmaps + + +class VFSimulation(Simulation): + def __init__(self, **kwargs): + """ + Inherited from Simulation class + :param phototaxis_theta_step: rotational speed scaling factor during + phototaxis + :param detection_range: detection range of resource patches (in pixels) + :param resource_meter_multiplier: scaling factor of how much resource is + extraxted for a detected resource unit + :param signalling_cost: cost of signalling in resource units + :param probability_of_starting_signaling: probability of starting signaling + :param des_velocity_res: desired velocity of resource patch in pixel per + timestep + :param res_theta_abs: change in orientation will be pulled from uniform + -res_theta_abs to res_theta_abs + :param agent_signaling_rand_event_update: updating agent's + random number for signalling probability in every N simulation time step + """ + super().__init__(**kwargs) + self.show_all_stats = False + # line following prototype + self.agent_for_line_id = [] + self.lines = [[]] # a list of guide lines that agents will follow where each line is a list of (x, y) points + self.line_map = np.zeros((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + + # making only the used part of retina to the given resolution + print(f"Original retina resolution: {self.v_field_res}") + self.v_field_res *= (1/self.fov_ratio) + self.v_field_res = int(self.v_field_res) + print(f"Due to limited FOV: {self.fov_ratio} increasing overall resolution to {self.v_field_res}") + + self.show_path_history = False + self.memory_length = 30 + self.ori_memory = None + self.pos_memory = None + + def update_lines_to_follow(self): + """Updating background line map to follow""" + subsurface = pygame.Surface((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + subsurface.fill(colors.BACKGROUND) + subsurface.set_colorkey(colors.WHITE) + subsurface.set_alpha(255) + for line in self.lines: + for pi in range(1, len(line)): + point1 = line[pi-1] + point2 = line[pi] + color = colors.BLACK + pygame.draw.line(subsurface, color, point1, point2, 4) + + for ag in self.agents: + if len(ag.lines)!=0: + sensor1_pos = [ag.position[1] + ag.radius - ag.sensor_distance + (1 + np.sin(ag.orientation + (3*np.pi / 4))) * ag.sensor_distance, + ag.position[0] + ag.radius - ag.sensor_distance + (1 - np.cos(ag.orientation + (3*np.pi / 4))) * ag.sensor_distance] + sensor2_pos = [ag.position[1] + ag.radius - ag.sensor_distance + (1 + np.sin(ag.orientation - (3*np.pi / 4))) * ag.sensor_distance, + ag.position[0] + ag.radius - ag.sensor_distance + (1 - np.cos(ag.orientation - (3*np.pi / 4))) * ag.sensor_distance] + + pygame.draw.circle( + subsurface, colors.GREEN, sensor2_pos, ag.sensor_size + ) + pygame.draw.circle( + subsurface, colors.GREEN, sensor1_pos, ag.sensor_size + ) + line_map = pygame.surfarray.array3d(subsurface) + self.line_map = line_map.swapaxes(0, 1)[:, :, 0]/255 + + + def handle_mouse_wheel_event(self, event): + """Handling event if mouse wheel is moved""" + if event.type == pygame.MOUSEWHEEL: + print(event.y) + keys = pygame.key.get_pressed() + if keys[pygame.K_a]: + for agid, ag in enumerate(list(self.agents)): + if ag.rect.collidepoint(pygame.mouse.get_pos()): + if ag.ALP0 is None: + ag.ALP0 = vf_params.ALP0 + else: + ag.ALP0 += event.y * 0.05 + elif keys[pygame.K_b]: + for agid, ag in enumerate(list(self.agents)): + if ag.rect.collidepoint(pygame.mouse.get_pos()): + if ag.BET0 is None: + ag.BET0 = vf_params.BET0 + else: + ag.BET0 += event.y * 0.05 + elif keys[pygame.K_r]: + for agid, ag in enumerate(list(self.agents)): + if ag.rect.collidepoint(pygame.mouse.get_pos()): + ag.radius += int(event.y) + elif keys[pygame.K_v]: + for agid, ag in enumerate(list(self.agents)): + if ag.rect.collidepoint(pygame.mouse.get_pos()): + if ag.V0 is None: + ag.V0 = vf_params.V0 + else: + ag.V0 += event.y * 0.05 + else: + if event.y == -1: + event.y = 0 + for ag in self.agents: + ag.move_with_mouse(pygame.mouse.get_pos(), event.y, 1 - event.y) + + + def handle_cursor_event(self, event): + """Handling event if cursor buttons are clicked overwriting base to add line functionalities""" + if pygame.mouse.get_pressed()[0]: + keys = pygame.key.get_pressed() + if keys[pygame.K_l]: + for agid, ag in enumerate(list(self.agents)): + if ag.rect.collidepoint(pygame.mouse.get_pos()): + self.agent_for_line_id.append(agid) + print(f"Assigned agent {agid} for current line, will finalize when release l button!") + # assigning the last line to the given agent + ag.color = colors.RED + # ag.lines.append(self.lines[-1]) + # print(ag.lines) + + # tart adding points to it + if (pygame.mouse.get_pos()[1], pygame.mouse.get_pos()[0]) not in self.lines[-1]: + self.lines[-1].append((pygame.mouse.get_pos()[1], pygame.mouse.get_pos()[0])) + print(len(self.lines), len(self.lines[-1])) + else: + try: + for ag in self.agents: + ag.move_with_mouse(event.pos, 0, 0) + for res in self.rescources: + res.update_clicked_status(event.pos) + except AttributeError: + for ag in self.agents: + ag.move_with_mouse(pygame.mouse.get_pos(), 0, 0) + else: + for ag in self.agents: + ag.is_moved_with_cursor = False + ag.draw_update() + for res in self.rescources: + res.is_clicked = False + res.draw_update() + + def draw_lines_to_follow(self): + """Draw lines that agents follow""" + self.update_lines_to_follow() + subsurface = pygame.surfarray.make_surface(self.line_map*255) + subsurface.set_colorkey(colors.WHITE) + subsurface.set_alpha(255) + self.screen.blit(subsurface, (0, 0)) + + + def draw_agent_stats(self, font_size=15, spacing=0): + """Showing agent information when paused""" + if self.show_all_stats: + font = pygame.font.Font(None, font_size) + for agent in self.agents: + a0 = agent.ALP0 if agent.ALP0 is not None else vf_params.ALP0 + b0 = agent.BET0 if agent.BET0 is not None else vf_params.BET0 + v0 = agent.V0 if agent.V0 is not None else vf_params.V0 + status = [ + f"ID: {agent.id}", + f"A0: {a0}", + f"B0: {b0}", + f"V0: {v0}", + ] + for i, stat_i in enumerate(status): + text = font.render(stat_i, True, colors.BLACK) + self.screen.blit( + text, + (agent.position[0] + 2 * agent.radius, + agent.position[1] + 2 * agent.radius + i * ( + font_size + spacing))) + + def add_new_agent(self, id, x, y, orient, with_proove=False, + behave_params=None): + """ + Adding a single new agent into agent sprites + """ + agent_proven = False + while not agent_proven: + agent = VFAgent( + id=id, + radius=self.agent_radii, + position=(x, y), + orientation=orient, + env_size=(self.WIDTH, self.HEIGHT), + color=colors.BLUE, + v_field_res=self.v_field_res, + FOV=self.agent_fov, + window_pad=self.window_pad, + pooling_time=self.pooling_time, + pooling_prob=self.pooling_prob, + consumption=self.agent_consumption, + vision_range=self.vision_range, + visual_exclusion=self.visual_exclusion, + patchwise_exclusion=self.patchwise_exclusion, + behave_params=None + ) + self.agents.add(agent) + agent_proven = True + + def create_agents(self): + """ + Creating agents according to how the simulation class was initialized + """ + for i in range(self.N): + orient = np.random.uniform(0, 2 * np.pi) + dist = np.random.rand() * ( + self.HEIGHT / 2 - 2 * self.window_pad - self.agent_radii) + x = np.cos(orient) * dist + self.WIDTH / 2 + y = np.sin(orient) * dist + self.HEIGHT / 2 + if not self.heterogen_agents: + # create agents according to environment variables homogeneously + self.add_new_agent(i, x, y, orient) + else: + self.add_new_agent( + i, x, y, orient, + behave_params=self.agent_behave_param_list[i]) + + def save_ori_pos_history(self): + """Saving orientation and position history of agents to visualize paths""" + if self.ori_memory is None: + self.ori_memory = np.zeros((len(self.agents), self.memory_length)) + self.pos_memory = np.zeros((len(self.agents), 2, self.memory_length)) + try: + self.ori_memory = np.roll(self.ori_memory, 1, axis=-1) + self.pos_memory = np.roll(self.pos_memory, 1, axis=-1) + self.ori_memory[:, 0] = np.array([ag.orientation for ag in self.agents]) + self.pos_memory[:, 0, 0] = np.array([ag.position[0]+ag.radius for ag in self.agents]) + self.pos_memory[:, 1, 0] = np.array([ag.position[1]+ag.radius for ag in self.agents]) + except: + self.ori_memory = None + self.pos_memory = None + + def draw_agent_paths_vf(self): + if self.ori_memory is not None: + path_length = self.memory_length + cmap = colmaps.get_cmap('jet') + transparency = 0.5 + transparency = int(transparency * 255) + big_colors = cmap(self.ori_memory / (2 * np.pi)) * 255 + # setting alpha + surface = pygame.Surface((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + surface.fill(colors.BACKGROUND) + surface.set_colorkey(colors.WHITE) + surface.set_alpha(255) + try: + for ai, agent in enumerate(self.agents): + subsurface = pygame.Surface((self.WIDTH + self.window_pad, self.HEIGHT + self.window_pad)) + subsurface.fill(colors.BACKGROUND) + subsurface.set_colorkey(colors.WHITE) + subsurface.set_alpha(transparency) + for t in range(2, path_length, 1): + point2 = self.pos_memory[ai, :, t] + color = big_colors[ai, t] + # pygame.draw.line(surface1, color, point1, point2, 4) + pygame.draw.circle(subsurface, color, point2, max(2, int(self.agent_radii / 2))) + surface.blit(subsurface, (0, 0)) + self.screen.blit(surface, (0, 0)) + except IndexError as e: + pass + + + def draw_frame(self, stats, stats_pos): + """Drawing environment, agents and every other visualization in each timestep""" + self.screen.fill(colors.BACKGROUND) + self.draw_lines_to_follow() + self.draw_walls() + if self.show_path_history: + self.draw_agent_paths_vf() + self.agents.draw(self.screen) + if self.show_vision_range: + self.draw_visual_fields() + self.draw_framerate() + self.draw_agent_stats() + + if self.show_vis_field: + # showing visual fields of the agents + self.show_visual_fields(stats, stats_pos) + + def step_sim(self): + """Stepping a single time step in the simulation loop""" + events = pygame.event.get() + # Carry out interaction according to user activity + self.interact_with_event(events) + # deciding if vis field needs to be shown in this timestep + self.turned_on_vfield = self.decide_on_vis_field_visibility( + self.turned_on_vfield) + + if not self.is_paused: + # Update agents according to current visible obstacles + self.agents.update(self.agents) + # move to next simulation timestep (only when not paused) + self.t += 1 + # Simulation is paused + else: + # Still calculating visual fields + for ag in self.agents: + # DEBUG: updating visual projections also when paused + ag.update_social_info(self.agents) + + if self.show_path_history: + self.save_ori_pos_history() + + # Draw environment and agents + if self.with_visualization: + self.draw_frame(self.stats, self.stats_pos) + pygame.display.flip() + # Monitoring with IFDB (also when paused) + if self.save_in_ifd: + ifdb.save_agent_data( + self.ifdb_client, self.agents, self.t, + exp_hash=self.ifdb_hash, + batch_size=self.write_batch_size) + ifdb.save_resource_data( + self.ifdb_client, self.rescources, self.t, + exp_hash=self.ifdb_hash, + batch_size=self.write_batch_size) + elif self.save_in_ram: + # saving data in ram for data processing, only when not paused + if not self.is_paused: + ifdb.save_agent_data_RAM(self.agents, self.t) + ifdb.save_resource_data_RAM(self.rescources, self.t) + # Moving time forward + if self.t % 500 == 0 or self.t == 1: + print( + f"{datetime.now().strftime('%Y-%m-%d_%H-%M-%S.%f')}" + f" t={self.t}") + print(f"Simulation FPS: {self.clock.get_fps()}") + self.clock.tick(self.framerate) + + def prepare_start(self): + """Preparing start method""" + print(f"Running simulation start method!") + # Creating N agents in the environment + print("Creating agents!") + self.create_agents() + # Creating surface to show visual fields + print("Creating visual field graph!") + self.stats, self.stats_pos = self.create_vis_field_graph() + # local var to decide when to show visual fields + self.turned_on_vfield = 0 + + + def start(self): + start_time = datetime.now() + self.prepare_start() + # Main Simulation loop until dedicated simulation time + print("Starting main simulation loop!") + while self.t < self.T: + self.step_sim() + end_time = datetime.now() + print( + f"{datetime.now().strftime('%Y-%m-%d_%H-%M-%S.%f')}" + f" Total simulation time: ", + (end_time - start_time).total_seconds()) + # Saving data from IFDB when simulation time is over + if self.agent_behave_param_list is not None: + if self.agent_behave_param_list[0].get( + "evo_summary_path") is not None: + pop_num = self.generate_evo_summary() + else: + pop_num = None + if self.save_csv_files: + if self.save_in_ifd or self.save_in_ram: + ifdb.save_ifdb_as_csv(exp_hash=self.ifdb_hash, + use_ram=self.save_in_ram, + as_zar=self.use_zarr, + save_extracted_vfield=False, + pop_num=pop_num) + env_saver.save_env_vars([self.env_path], "env_params.json", + pop_num=pop_num) + else: + raise Exception( + "Tried to save simulation data as csv file due to env " + "configuration, " + "but IFDB/RAM logging was turned off. Nothing to save! " + "Please turn on IFDB/RAM logging" + " or turn off CSV saving feature.") + end_save_time = datetime.now() + print( + f"{datetime.now().strftime('%Y-%m-%d_%H-%M-%S.%f')} " + f"Total saving time:", + (end_save_time - end_time).total_seconds()) + + pygame.quit() + # sys.exit() diff --git a/abm/projects/visual_flocking/vf_utils/attraction_repulsion_map.py b/abm/projects/visual_flocking/vf_utils/attraction_repulsion_map.py new file mode 100644 index 00000000..9be78286 --- /dev/null +++ b/abm/projects/visual_flocking/vf_utils/attraction_repulsion_map.py @@ -0,0 +1,102 @@ +import numpy as np + +from abm.projects.visual_flocking.vf_simulation.vf_sims import VFSimulation + +# VF_GAM=0.1 +# VF_V0=0 +# VF_ALP0=1 +# VF_ALP1=0.09 +# VF_ALP2=0 +# VF_BET0=1 +# VF_BET1=0.09 +# VF_BET2=0 + +width = 800 +height = 800 +undersample = 2 +radius = 5 + +foc_ag_pos = (int(width/2)-radius, int(height/2)-radius) +dv_matrix = np.zeros((int(width/undersample), int(height/undersample))) +dpsi_matrix = np.zeros((int(width/undersample), int(height/undersample))) +ablob_matrix = np.zeros((int(width/undersample), int(height/undersample))) +aedge_matrix = np.zeros((int(width/undersample), int(height/undersample))) +bblob_matrix = np.zeros((int(width/undersample), int(height/undersample))) +bedge_matrix = np.zeros((int(width/undersample), int(height/undersample))) + +sim = VFSimulation(N=2, + T=width*height+1, + v_field_res=2000, + agent_fov=1, + framerate=2000, + with_visualization=0, + width=width, + height=height, + show_vis_field=True, + show_vis_field_return=False, + vision_range=2000, + visual_exclusion=False, + show_vision_range=True, + use_ifdb_logging=False, + use_ram_logging=True, + save_csv_files=False, + use_zarr=True, + parallel=False, + window_pad=0, + agent_behave_param_list=None, + collide_agents=False) + +sim.prepare_start() +for w in range(0, width, undersample): + print(f"progress: {w}/{width}") + for h in range(0, height, undersample): + agent0 = list(sim.agents)[0] + agent1 = list(sim.agents)[1] + agent0.position = np.array(foc_ag_pos, dtype=np.float64) + agent0.orientation = np.pi / 2 + agent0.velocity = 0 + agent0.verbose_supcalc = True + agent1.position = np.array((w, h), dtype=np.float64) + sim.step_sim() + dv_matrix[int(w/undersample), int(h/undersample)] = agent0.dv + dpsi_matrix[int(w/undersample), int(h/undersample)] = agent0.dphi + ablob_matrix[int(w/undersample), int(h/undersample)] = agent0.ablob + aedge_matrix[int(w / undersample), int(h / undersample)] = agent0.aedge + bblob_matrix[int(w/undersample), int(h/undersample)] = agent0.bblob + bedge_matrix[int(w / undersample), int(h / undersample)] = agent0.bedge + +dv_matrix[dv_matrix>0.2] = 0.2 +dv_matrix[dv_matrix<-0.2] = -0.2 +ablob_matrix[ablob_matrix>0.2] = 0.2 +ablob_matrix[ablob_matrix<-0.2] = -0.2 +aedge_matrix[aedge_matrix>0.2] = 0.2 +aedge_matrix[aedge_matrix<-0.2] = -0.2 + +dpsi_matrix[dpsi_matrix>0.2] = 0.2 +dpsi_matrix[dpsi_matrix<-0.2] = -0.2 +bblob_matrix[bblob_matrix>0.2] = 0.2 +bblob_matrix[bblob_matrix<-0.2] = -0.2 +bedge_matrix[bedge_matrix>0.2] = 0.2 +bedge_matrix[bedge_matrix<-0.2] = -0.2 + +import matplotlib.pyplot as plt +import matplotlib + +cmap_vel = matplotlib.colors.LinearSegmentedColormap.from_list( + 'mycmap', ["red", "white", "green"]) +cmap_ori = matplotlib.colors.LinearSegmentedColormap.from_list( + 'mycmap', ["orange", "white", "blue"]) +plt.figure() +plt.imshow(dv_matrix.T, cmap=cmap_vel) +plt.show() +plt.figure() +plt.imshow(dpsi_matrix.T, cmap=cmap_ori) +plt.show() +plt.imshow(ablob_matrix.T, cmap=cmap_vel) +plt.show() +plt.imshow(aedge_matrix.T, cmap=cmap_vel) +plt.show() +plt.imshow(bblob_matrix.T, cmap=cmap_ori) +plt.show() +plt.imshow(bedge_matrix.T, cmap=cmap_ori) +plt.show() diff --git a/abm/projects/visual_flocking/vf_utils/blobs_vs_edges.py b/abm/projects/visual_flocking/vf_utils/blobs_vs_edges.py new file mode 100644 index 00000000..766f64d0 --- /dev/null +++ b/abm/projects/visual_flocking/vf_utils/blobs_vs_edges.py @@ -0,0 +1,58 @@ +import matplotlib.pyplot as plt +import numpy as np + +# creating figure and axes +fig, ax = plt.subplots(2, 2, sharex=True, sharey=True) + +# creating x axis +axis_len = 2000 +phi = np.arange(-np.pi, np.pi, (2*np.pi)/2000) + +# creating visual blobs +blob_widths = [600, 150] +offsets = [0, int(axis_len/4)] +centers = [int(axis_len/2) + offset for offset in offsets] + +# define trigonometric functions +trigs = [np.cos, np.sin] + +# define colors +attr_colors = ["green", "blue"] +repul_colors = ["red", "orange"] + +for casei, case in enumerate(["front-back", "left-right"]): + for disti, dis in enumerate(["near", "far"]): + # choosing axis + showax = ax[casei, disti] + plt.axes(showax) + + # creating blobs + VPF = np.zeros(axis_len) + blob_width = blob_widths[disti] + center = centers[casei] + bstart = center - int(blob_width/2) + bend = center + int(blob_width/2) + VPF[bstart:bend] = 1 + + # creating edges + VPF_edge = np.zeros(axis_len) + VPF_edge[bstart] = 1 + VPF_edge[bend] = 1 + + # modulating with trig function + VPF_mod = VPF * trigs[casei](phi) + VPF_edge_mod = VPF_edge * trigs[casei](phi) + + # plotting + plt.fill_between(phi, np.zeros(axis_len), VPF_mod, alpha=0.5, color=repul_colors[casei]) + plt.plot(phi, VPF, c="gray", ls="-") + plt.plot(phi, trigs[casei](phi), c="gray", ls="--") + + # plt.plot(phi, VPF_edge_mod, c="green") + plt.arrow(phi[bstart], 0, 0, trigs[casei](phi[bstart]), width=0.05, length_includes_head=True, edgecolor=attr_colors[casei], facecolor=attr_colors[casei]) + plt.arrow(phi[bend], 0, 0, trigs[casei](phi[bend]), width=0.05, length_includes_head=True, edgecolor=attr_colors[casei], facecolor=attr_colors[casei]) + # plt.plot(phi, VPF_mod) + + +plt.subplots_adjust(hspace=0.2, wspace=0) +plt.show() \ No newline at end of file diff --git a/abm/projects/visual_flocking/visual_flocking.env b/abm/projects/visual_flocking/visual_flocking.env new file mode 100644 index 00000000..b7907fd5 --- /dev/null +++ b/abm/projects/visual_flocking/visual_flocking.env @@ -0,0 +1,121 @@ +# Basic simulation parameters + +# Number of agents and resc. patches +N = 5 +N_RESOURCES = 1 + +# Simulation time +T = 1000 + +# Framerate (only matter if visualization is turned on) +INIT_FRAMERATE = 25 + +# Visualization ON/OFF +WITH_VISUALIZATION = 1 + +# Resolution (px) of agents' visual projection fields +VISUAL_FIELD_RESOLUTION = 1200 + +# Enviornment size (width and height in px) +ENV_WIDTH = 600 +ENV_HEIGHT = 600 + +# Agent and Resource sizes +RADIUS_AGENT = 10 +RADIUS_RESOURCE = 20 +MIN_RESOURCE_PER_PATCH = 400 +MAX_RESOURCE_PER_PATCH = 450 +REGENERATE_PATCHES = 1 +PATCH_BORDER_OVERLAP = 1 +MAX_SPEED = 1.5 + +# Exploitation +# max consumption of an agent per timestep +AGENT_CONSUMPTION = 1 +# resource patch quality +MIN_RESOURCE_QUALITY = 0.05 +MAX_RESOURCE_QUALITY = 0.3 +TELEPORT_TO_MIDDLE = 0 +GHOST_WHILE_EXPLOIT = 1 +AGENT_AGENT_COLLISION = 1 +PATCHWISE_SOCIAL_EXCLUSION =1 + +# Vision +# fov is symmetric and defined with percent of pi. e.g +# if 0.6 then fov is (-0.6*pi, 0.6*pi). 1 is full 360 degree vision +AGENT_FOV = 0.45 +VISION_RANGE = 200 +VISUAL_EXCLUSION = 1 + +# Visualize agents' visual field on screen by default or when return is pressed +SHOW_VISUAL_FIELDS = 0 +SHOW_VISUAL_FIELDS_RETURN = 1 +SHOW_VISION_RANGE = 0 + +# Time needed to acquire info from env and initial probability for pooling process +POOLING_TIME = 0 +POOLING_PROBABILITY = 0.05 + +# Monitoring Related parameters (True only tested on Ubuntu machines with +# initialized grafana and ifdb instances as in readme) +# log data in influxdb +USE_IFDB_LOGGING = 0 +# logs data in RAM +USE_RAM_LOGGING = 1 +# using zarr compression +USE_ZARR_FORMAT = 1 +# saves csv files (if ifdb logging is used, json files otherwise) +SAVE_CSV_FILES = 0 + +# Decision process parameters +DEC_TW = 0.5 +DEC_EPSW = 3 +DEC_GW = 0.085 +DEC_BW = 0 +DEC_WMAX = 1 + +DEC_TU = 0.5 +DEC_EPSU = 3 +DEC_GU = 0.085 +DEC_BU = 0 +DEC_UMAX = 1 + +DEC_SWU = 0.25 +DEC_SUW = 0.01 + +DEC_TAU = 10 +DEC_FN = 2 +DEC_FR = 1 + +# Movement parameters +# Exploration +MOV_EXP_VEL_MIN = 1 +MOV_EXP_VEL_MAX = 1 +MOV_EXP_TH_MIN = -0.3 +MOV_EXP_TH_MAX = 0.3 + +# Relocation +MOV_REL_DES_VEL = 1 +MOV_REL_TH_MAX = 0.5 +# Length of memory for signaling target +MEMORY_DEPTH = 1 + +# Exploitation +CONS_STOP_RATIO = 0.08 + +##### Visual Flocking ##### +# In case the app is started using only .env +# environment (headless or without playground) +# this will enable to start the simulation +# otherwise running the project specific +# app file is disabled to avoid mix up between projects +APP_VERSION = VisualFlocking + +VF_GAM = 0.05 +VF_V0 = 1 +VF_ALP0 = 3 +VF_ALP1 = 0.0012 +VF_ALP2 = 0 +VF_BET0 = 5 +VF_BET1 = 0.0012 +VF_BET2 = 0 \ No newline at end of file diff --git a/abm/replay/replay.py b/abm/replay/replay.py index 4b9b9688..72f84280 100644 --- a/abm/replay/replay.py +++ b/abm/replay/replay.py @@ -16,6 +16,7 @@ from abm.contrib import playgroundtool as pgt import shutil import cv2 +from matplotlib import cm as colmaps root_abm_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) @@ -68,7 +69,7 @@ def __init__(self, data_folder_path, undersample=1, t_start=None, t_end=None, co if self.project_version=="Base": self.coll_resc_z = self.experiment.agent_summary['collresource'] self.resc_left_z = self.experiment.res_summary['resc_left'] - self.resc_quality_z = self.experiment.res_summary['quality'] + self.resc_quality_z = self.experiment.res_summary.get('quality', None) elif self.project_version=="CooperativeSignaling": self.meter_z = self.experiment.agent_summary['meter'] self.sig_z = self.experiment.agent_summary['signalling'] @@ -262,6 +263,26 @@ def __init__(self, data_folder_path, undersample=1, t_start=None, t_end=None, co onClick=lambda: self.take_snapshot(), # Function to call when clicked on borderThickness=1 ) + self.NNdist_b = Button( + # Mandatory Parameters + self.screen, # Surface to place button on + self.button_start_x_2 + int(self.slider_width / 2), # X-coordinate of top left corner + button_start_y, # Y-coordinate of top left corner + int(self.slider_width / 2), # Width + self.button_height, # Height + + # Optional Parameters + text='Avg. NN Dist.', # Text to display + fontSize=20, # Size of font + margin=20, # Minimum distance between text/image and edge of button + inactiveColour=colors.GREY, + onClick=lambda: self.on_print_NNdist(), # Function to call when clicked on + borderThickness=1 + ) + if self.project_version in ["Base", "CooperativeSignaling"]: + self.NNdist_b.disable() + self.NNdist_b.hide() + # Plotting Button Line button_start_y += self.button_height @@ -281,6 +302,27 @@ def __init__(self, data_folder_path, undersample=1, t_start=None, t_end=None, co onClick=lambda: self.on_print_efficiency(), # Function to call when clicked on borderThickness=1 ) + + if self.experiment.env.get("APP_VERSION") == "VisualFlocking": + self.plot_efficiency.disable() + self.plot_efficiency.hide() + self.plot_polar = Button( + # Mandatory Parameters + self.screen, # Surface to place button on + self.slider_start_x, # X-coordinate of top left corner + button_start_y, # Y-coordinate of top left corner + int(self.slider_width / 2), # Width + self.button_height, # Height + + # Optional Parameters + text='Plot Polarization', # Text to display + fontSize=20, # Size of font + margin=20, # Minimum distance between text/image and edge of button + inactiveColour=colors.GREY, + onClick=lambda: self.on_print_polarization(), # Function to call when clicked on + borderThickness=1 + ) + self.plot_rel_time = Button( # Mandatory Parameters self.screen, # Surface to place button on @@ -297,6 +339,26 @@ def __init__(self, data_folder_path, undersample=1, t_start=None, t_end=None, co onClick=lambda: self.on_print_reloc_time(), # Function to call when clicked on borderThickness=1 ) + if self.experiment.env.get("APP_VERSION") == "VisualFlocking": + self.plot_rel_time.disable() + self.plot_rel_time.hide() + self.plot_coll = Button( + # Mandatory Parameters + self.screen, # Surface to place button on + self.button_start_x_2, # X-coordinate of top left corner + button_start_y, # Y-coordinate of top left corner + int(self.slider_width / 2), # Width + self.button_height, # Height + + # Optional Parameters + text='Plot Collisions', # Text to display + fontSize=20, # Size of font + margin=20, # Minimum distance between text/image and edge of button + inactiveColour=colors.GREY, + onClick=lambda: self.on_print_collision_times(), # Function to call when clicked on + borderThickness=1 + ) + if len(list(self.experiment.varying_params.keys())) in [3, 4]: self.collapse_dropdown = Dropdown( self.screen, @@ -467,7 +529,10 @@ def load_data_to_ram(self): if self.project_version=="Base": self.coll_resc = self.coll_resc_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] self.resc_left = self.resc_left_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] - self.resc_quality = self.resc_quality_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] + if self.resc_quality_z is not None: + self.resc_quality = self.resc_quality_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] + else: + self.resc_quality = None elif self.project_version=="CooperativeSignaling": self.meter = self.meter_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] self.sig = self.sig_z[self.index][:, (self.t_slice) * cs:(self.t_slice + 1) * cs] @@ -479,6 +544,46 @@ def on_print_reloc_time(self): self.experiment.set_collapse_param(self.collapse_dropdown.getSelected()) self.experiment.plot_mean_relocation_time() + def on_print_collision_times(self, with_read_collapse_param=True, used_batches=None): + if with_read_collapse_param: + if len(list(self.experiment.varying_params.keys())) in [3, 4]: + self.experiment.set_collapse_param(self.collapse_dropdown.getSelected()) + if self.t_start is None: + t_start = 0 + else: + t_start = self.t_start + if self.t_end is None: + t_end = self.T - 1 + else: + t_end = self.t_end + if self.T > 1000: + undersample = int(self.T / 1000) + print(f"Experiment longer than 1000 timesteps! To calculate iid reducing timesteps to 1000 with undersampling rate {undersample}.") + else: + undersample = 1 + fig, ax, cbar = self.experiment.plot_mean_collision_time(t_start=t_start, t_end=t_end, + from_script=self.from_script, + used_batches=used_batches, + undersample=undersample) + return fig, ax, cbar + + def on_print_polarization(self, with_read_collapse_param=True, used_batches=None): + if with_read_collapse_param: + if len(list(self.experiment.varying_params.keys())) in [3, 4]: + self.experiment.set_collapse_param(self.collapse_dropdown.getSelected()) + if self.t_start is None: + t_start = 0 + else: + t_start = self.t_start + if self.t_end is None: + t_end = self.T - 1 + else: + t_end = self.t_end + fig, ax, cbar = self.experiment.plot_mean_polarization(t_start=t_start, t_end=t_end, + from_script=self.from_script, + used_batches=used_batches) + return fig, ax, cbar + def on_print_efficiency(self, with_read_collapse_param=True, used_batches=None): """print mean search efficiency""" if with_read_collapse_param: @@ -510,6 +615,20 @@ def on_print_iid(self, with_read_collapse_param=True, used_batches=None): fig, ax, cbar = self.experiment.plot_mean_iid(from_script=self.from_script, undersample=undersample) return fig, ax, cbar + def on_print_NNdist(self, with_read_collapse_param=True, used_batches=None): + """print mean inter-individual distance""" + if with_read_collapse_param: + if len(list(self.experiment.varying_params.keys())) in [3, 4]: + self.experiment.set_collapse_param(self.collapse_dropdown.getSelected()) + if self.T > 1000: + undersample = int(self.T / 1000) + print( + f"Experiment longer than 1000 timesteps! To calculate iid reducing timesteps to 1000 with undersampling rate {undersample}.") + else: + undersample = 1 + fig, ax, cbar = self.experiment.plot_mean_NN_dist(from_script=self.from_script, undersample=undersample) + return fig, ax, cbar + def on_set_t_start(self): """Setting starting timestep for plotting and calculations""" if self.t_start is None: @@ -614,7 +733,6 @@ def draw_separator(self): def draw_frame(self, events): """Drawing environment, agents and every other visualization in each timestep""" self.screen.fill(colors.BACKGROUND) - self.draw_separator() pygame_widgets.update(events) self.framerate = self.framerate_slider.getValue() @@ -710,10 +828,22 @@ def update_frame_data(self): res_posx = self.res_pos_x[:, t_ind] res_posy = self.res_pos_y[:, t_ind] - if self.project_version=="Base": + if self.project_version in ["Base", "CooperativeSignaling"]: coll_resc = self.coll_resc[:, t_ind] - resc_left = self.resc_left[:, t_ind] - resc_quality = self.resc_quality[:, t_ind] + if self.project_version == "Base": + if self.resc_quality is not None: + resc_quality = self.resc_quality[:, t_ind] + else: + resc_quality = [0 for i in range(len(res_posx))] + resc_left = self.resc_left[:, t_ind] + else: + resc_quality = None + resc_left = None + elif self.project_version == "VisualFlocking": + coll_resc = [0 for i in range(len(posx))] + resc_left = [0 for i in range(len(posx))] + resc_quality = [0 for i in range(len(posx))] + res_unit = self.env["MIN_RESOURCE_PER_PATCH"] if res_unit == "----TUNED----": @@ -744,19 +874,47 @@ def update_frame_data(self): indexalongdim = slider.getValue() res_radius = self.varying_params[radius_keyword][indexalongdim] + # getting FOV if it is varying + if self.env["AGENT_FOV"] == "----TUNED----": + var_keys = sorted(list(self.varying_params.keys())) + dimnum = var_keys.index("AGENT_FOV") + slider = self.varying_sliders[dimnum] + indexalongdim = slider.getValue() + self.current_fov = float(self.varying_params["AGENT_FOV"][indexalongdim]) + else: + self.current_fov = float(self.env["AGENT_FOV"]) + + if self.project_version=="Base": self.draw_resources(res_posx, res_posy, max_units, resc_left, resc_quality, res_radius) elif self.project_version=="CooperativeSignaling": self.draw_resources(res_posx, res_posy, [1], [1], [1], res_radius) if self.show_paths: - self.draw_agent_paths(self.posx[:, max(0, t_ind - self.path_length):t_ind], - self.posy[:, max(0, t_ind - self.path_length):t_ind], - radius, - modes=self.agmodes[:, max(0, t_ind - self.path_length):t_ind]) + if self.experiment.env.get("APP_VERSION", "") == "VisualFlocking": + self.draw_agent_paths_vf(self.posx[:, max(0, t_ind - self.path_length):t_ind], + self.posy[:, max(0, t_ind - self.path_length):t_ind], + radius, + self.orientation[:, max(0, t_ind - self.path_length):t_ind]) + else: + self.draw_agent_paths(self.posx[:, max(0, t_ind - self.path_length):t_ind], + self.posy[:, max(0, t_ind - self.path_length):t_ind], + radius, + modes=self.agmodes[:, max(0, t_ind - self.path_length):t_ind]) + if self.experiment.clustering_data is not None: + t_ind_cl = int(t_ind / 25) + clusters_idx = tuple(list(self.index) + [slice(None), t_ind_cl]) + clusters = self.experiment.clustering_ids[clusters_idx] + print(clusters) + else: + clusters = None + + if self.project_version=="Base": self.draw_agents(posx, posy, orientation, mode, coll_resc, radius) elif self.project_version == "CooperativeSignaling": self.draw_agents(posx, posy, orientation, mode, [0 for i in range(len(posx))], radius) + elif self.project_version == "VisualFlocking": + self.draw_agents(posx, posy, orientation, mode, [0 for i in range(len(posx))], radius, clusters=clusters) num_agents = len(posx) if self.show_stats: @@ -770,9 +928,43 @@ def update_frame_data(self): if np.isnan(std_reloc): std_reloc = 0 time_dep_stats.append(f"Relocation Time (0-t): Mean:{mean_reloc:10.2f} ± {std_reloc:10.2f}") - end_pos = self.draw_agent_stat_summary([ai for ai in range(num_agents)], posx, posy, orientation, mode, - coll_resc, previous_metrics=time_dep_stats) - self.draw_resource_stat_summary(posx, posy, max_units, resc_left, resc_quality, end_pos) + if self.project_version in ["Base", "CooperativeSignaling"]: + end_pos = self.draw_agent_stat_summary([ai for ai in range(num_agents)], posx, posy, orientation, mode, + coll_resc, previous_metrics=time_dep_stats) + elif self.project_version == "VisualFlocking": + end_pos = self.draw_agent_stat_summary([ai for ai in range(num_agents)], posx, posy, orientation, mode, + [0 for i in range(len(posx))], previous_metrics=time_dep_stats) + if self.project_version == "Base": + self.draw_resource_stat_summary(posx, posy, max_units, resc_left, resc_quality, end_pos) + + def draw_agent_paths_vf(self, posx, posy, radius, orientations): + num_agents = posx.shape[0] + path_length = posx.shape[1] + cmap = colmaps.get_cmap('jet') + transparency = 0.5 + transparency = int(transparency * 255) + big_colors = cmap(orientations/(2 * np.pi))*255 + # setting alpha + surface = pygame.Surface((self.WIDTH+self.window_pad, self.HEIGHT+self.window_pad)) + surface.fill(colors.BACKGROUND) + surface.set_colorkey(colors.WHITE) + surface.set_alpha(255) + try: + for ai in range(num_agents): + subsurface = pygame.Surface((self.WIDTH+self.window_pad, self.HEIGHT+self.window_pad)) + subsurface.fill(colors.BACKGROUND) + subsurface.set_colorkey(colors.WHITE) + subsurface.set_alpha(transparency) + for t in range(1, path_length, 1): + #point1 = (posx[ai, t - 1] + radius, posy[ai, t - 1] + radius) + point2 = (posx[ai, t] + radius, posy[ai, t] + radius) + color = big_colors[ai, t] + # pygame.draw.line(surface1, color, point1, point2, 4) + pygame.draw.circle(subsurface, color, point2, int(self.radius/2)) + surface.blit(subsurface, (0, 0)) + self.screen.blit(surface, (0, 0)) + except IndexError as e: + pass def draw_agent_paths(self, posx, posy, radius, modes=None): path_cols = [colors.BLUE, colors.GREEN, colors.PURPLE, colors.RED] @@ -823,11 +1015,30 @@ def draw_res_patch(self, id, posx, posy, max_unit, resc_left, resc_quality, radi tbsize = text.get_size() self.screen.blit(text, (int(posx + radius - tbsize[0] / 2), int(posy + radius - tbsize[1] / 2))) - def draw_agents(self, posx, posy, orientation, mode, coll_resc, radius): + def draw_agents(self, posx, posy, orientation, mode, coll_resc, radius, clusters=None): """Drawing agents in arena according to data""" + if self.experiment.env.get("APP_VERSION") == "VisualFlocking": + cmap = colmaps.get_cmap('jet') + if clusters is None: + colors = [np.array(cmap(ori / (2 * np.pi)))[0:-1] for ori in orientation] + # rescaling color for pygame + for col in colors: + col[0:3] *= 255 + else: + # defining colors according to cluster id_s in the clusters array + colors = [np.array(cmap(cl / (np.max(clusters))))[0:-1] for cl in clusters] + # making the colors brighter + colors = [col * 0.5 + 0.5 for col in colors] + # rescaling color for pygame + for col in colors: + col[0:3] *= 255 + else: + colors = [self.mode_to_color(m) for m in mode] + num_agents = len(posx) for ai in range(num_agents): - self.draw_agent(ai, posx[ai], posy[ai], orientation[ai], mode[ai], coll_resc[ai], radius) + cid = clusters[ai] if clusters is not None else None + self.draw_agent(ai, posx[ai], posy[ai], orientation[ai], colors[ai], coll_resc[ai], radius, cluster_id=cid) if self.show_vfield: self.draw_vfield(ai, posx[ai], posy[ai], orientation[ai], radius) @@ -906,11 +1117,11 @@ def draw_resource_stat_summary(self, posx, posy, max_units, resc_left, resc_qual def draw_vfield(self, id, posx, posy, orientation, radius): """Drawing a single agent according to position and orientation""" - FOV_rat = float(self.env.get("AGENT_FOV", 1)) + FOV_rat = self.current_fov FOV = (-FOV_rat * np.pi, FOV_rat * np.pi) # Show limits of FOV - if FOV[1] < np.pi: + if 0 < FOV[1] < np.pi: # Center and radius of pie chart cx, cy, r = posx + radius, posy + radius, self.env.get("VISION_RANGE", 3 * radius) @@ -941,12 +1152,13 @@ def draw_vfield(self, id, posx, posy, orientation, radius): self.screen.blit(image, (0, 0)) - def draw_agent(self, id, posx, posy, orientation, mode, coll_resc, radius): + def draw_agent(self, id, posx, posy, orientation, agent_color, coll_resc, radius, cluster_id=None): """Drawing a single agent according to position and orientation""" + # todo: add magnify agent checkbox + radius = int(radius) image = pygame.Surface([radius * 2, radius * 2]) image.fill(colors.BACKGROUND) image.set_colorkey(colors.BACKGROUND) - agent_color = self.mode_to_color(mode) pygame.gfxdraw.filled_circle(image, radius, radius, radius-1, agent_color) pygame.gfxdraw.aacircle(image, radius, radius, radius - 1, colors.BLACK) @@ -957,7 +1169,13 @@ def draw_agent(self, id, posx, posy, orientation, mode, coll_resc, radius): if self.show_stats: font = pygame.font.Font(None, 16) - text = font.render(f"ID:{id}, R:{coll_resc:.2f}", True, colors.BLACK) + if self.env.get("APP_VERSION", "Base") in ["Base", "CooperativeSignaling"]: + text = font.render(f"ID:{id}, R:{coll_resc:.2f}", True, colors.BLACK) + elif self.experiment.env.get("APP_VERSION") == "VisualFlocking": + if cluster_id is None: + text = font.render(f"ID:{id}, ori:{orientation:.2f}", True, colors.BLACK) + else: + text = font.render(f"ID:{id}, ori:{orientation:.2f}, cl:{cluster_id}", True, colors.BLACK) self.screen.blit(text, (posx - 10, posy - 10)) def interact_with_event(self, event): diff --git a/abm/replay_exec.py b/abm/replay_exec.py new file mode 100644 index 00000000..79fd2ba1 --- /dev/null +++ b/abm/replay_exec.py @@ -0,0 +1,25 @@ +'''General replay script to browse and replay experiments from a folder.''' + +from abm.replay.replay import ExperimentReplay + +# opening file browser to choose folder of interest +from tkinter import Tk +from tkinter.filedialog import askdirectory +Tk().withdraw() + +# get current directory to set as starting directory for file browsing +import os +current_dir = os.getcwd() +folder_path = askdirectory(title="Choose folder with experiment data", initialdir=current_dir) + +if folder_path == "": + print("No folder selected, exiting...") + exit() + +# try: +loaded_experiment = ExperimentReplay(folder_path) +loaded_experiment.start() +# except Exception as e: +# print("Error during experiment loading:", e) +# print("Please check if the folder contains all necessary files to be considered as an experiment.") +# exit() \ No newline at end of file diff --git a/abm/simulation/isims.py b/abm/simulation/isims.py index 6852ebdf..c3d7f8fc 100644 --- a/abm/simulation/isims.py +++ b/abm/simulation/isims.py @@ -147,6 +147,17 @@ def __init__(self): inactiveColour=colors.GREY, borderThickness=1, onClick=lambda: self.start_stop_IFDB_logging()) self.function_buttons.append(self.IFDB_button) + # hiding ifdb button + self.IFDB_button.hide() + + function_button_start_x += self.function_button_width + self.function_button_pad + self.Snapshot_button = Button(self.screen, function_button_start_x, function_button_start_y, + self.function_button_width, + self.function_button_height, text='Take Snapshot', + fontSize=self.function_button_height - 2, + inactiveColour=colors.GREY, borderThickness=1, + onClick=lambda: self.take_snapshot()) + self.function_buttons.append(self.Snapshot_button) self.global_stats_start += 2 * self.function_button_height + self.function_button_pad + self.window_pad @@ -305,6 +316,33 @@ def __init__(self): self.sliders.append(self.SUMR_slider) self.slider_texts.append(self.SUMR_textbox) + def take_snapshot(self): + """Taking a single picture of the current status of the replay into an image""" + filename = f"{pad_to_n_digits(self.t, n=6)}.png" + path = os.path.join(self.image_save_path, filename) + pygame.image.save(self.screen, path) + # cropping image + img = cv2.imread(path) + src = img[0:self.vis_area_end_height, 0:self.vis_area_end_width] + # Convert image to image gray + tmp = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY) + + # Applying thresholding technique + _, alpha = cv2.threshold(tmp, 254, 255, cv2.THRESH_BINARY_INV) + + # Using cv2.split() to split channels + # of coloured image + b, g, r = cv2.split(src) + + # Making list of Red, Green, Blue + # Channels and alpha + rgba = [b, g, r, alpha] + + # Using cv2.merge() to merge rgba + # into a coloured/multi-channeled image + dst = cv2.merge(rgba, 4) + cv2.imwrite(path, dst) + def start_stop_IFDB_logging(self): """Start or stop IFDB logging in case of grafana interface is used""" self.save_in_ifd = not self.save_in_ifd @@ -321,6 +359,7 @@ def start_stop_IFDB_logging(self): def change_ghost_mode(self): """Changing ghost mdoe during exploutation""" self.ghost_mode = not self.ghost_mode + self.collide_agents = not self.ghost_mode if self.ghost_mode: self.ghost_mode_button.inactiveColour = colors.GREEN else: @@ -506,6 +545,10 @@ def interact_with_event(self, events): """Carry out functionality according to user's interaction""" super().interact_with_event(events) pygame_widgets.update(events) + for event in events: + if event.type == pygame.KEYDOWN and event.key == pygame.K_p: + print("Snapshot!") + self.take_snapshot() self.framerate = self.framerate_slider.getValue() self.N = self.N_slider.getValue() self.N_resc = self.NRES_slider.getValue() diff --git a/abm/simulation/sims.py b/abm/simulation/sims.py index 327a5a5d..3e592a2b 100644 --- a/abm/simulation/sims.py +++ b/abm/simulation/sims.py @@ -273,7 +273,7 @@ def draw_visual_fields(self): """Visualizing the range of vision for agents as opaque circles around the agents""" for agent in self.agents: # Show visual range - pygame.draw.circle(self.screen, colors.LIGHT_BLUE, agent.position + agent.radius, agent.vision_range, + pygame.draw.circle(self.screen, colors.LIGHT_BLUE, np.array(agent.position) + agent.radius, agent.vision_range, width=1) # Show limits of FOV @@ -567,11 +567,7 @@ def interact_with_event(self, events): sys.exit() # Change orientation with mouse wheel - if event.type == pygame.MOUSEWHEEL: - if event.y == -1: - event.y = 0 - for ag in self.agents: - ag.move_with_mouse(pygame.mouse.get_pos(), event.y, 1 - event.y) + self.handle_mouse_wheel_event(event) # Pause on Space if event.type == pygame.KEYDOWN and event.key == pygame.K_SPACE: @@ -590,22 +586,34 @@ def interact_with_event(self, events): self.framerate = self.framerate_orig # Continuous mouse events (move with cursor) - if pygame.mouse.get_pressed()[0]: - try: - for ag in self.agents: - ag.move_with_mouse(event.pos, 0, 0) - for res in self.rescources: - res.update_clicked_status(event.pos) - except AttributeError: - for ag in self.agents: - ag.move_with_mouse(pygame.mouse.get_pos(), 0, 0) - else: + self.handle_cursor_event(event) + + def handle_mouse_wheel_event(self, event): + """Handling event if mouse wheel is moved""" + if event.type == pygame.MOUSEWHEEL: + if event.y == -1: + event.y = 0 + for ag in self.agents: + ag.move_with_mouse(pygame.mouse.get_pos(), event.y, 1 - event.y) + + def handle_cursor_event(self, event): + """Handling event if cursor buttons are clicked""" + if pygame.mouse.get_pressed()[0]: + try: for ag in self.agents: - ag.is_moved_with_cursor = False - ag.draw_update() + ag.move_with_mouse(event.pos, 0, 0) for res in self.rescources: - res.is_clicked = False - res.draw_update() + res.update_clicked_status(event.pos) + except AttributeError: + for ag in self.agents: + ag.move_with_mouse(pygame.mouse.get_pos(), 0, 0) + else: + for ag in self.agents: + ag.is_moved_with_cursor = False + ag.draw_update() + for res in self.rescources: + res.is_clicked = False + res.draw_update() def decide_on_vis_field_visibility(self, turned_on_vfield): """Deciding f the visual field needs to be shown or not""" diff --git a/setup.py b/setup.py index cc2489f3..49e3e72a 100644 --- a/setup.py +++ b/setup.py @@ -23,7 +23,8 @@ 'influxdb<5.3.0', 'opencv-python', 'xvfbwrapper', - 'zarr' + 'zarr', + 'fastcluster' ], extras_require={ 'test': [ @@ -37,9 +38,11 @@ 'console_scripts': [ 'abm-start=abm.app:start', 'abm-start-coopsig=abm.app_collective_signaling:start', + 'abm-start-visflock=abm.app_visual_flocking:start', 'headless-abm-start=abm.app:start_headless', 'playground-start=abm.app:start_playground', - 'playground-start-coopsig=abm.app_collective_signaling:start_playground' + 'playground-start-coopsig=abm.app_collective_signaling:start_playground', + 'playground-start-visflock=abm.app_visual_flocking:start_playground' ] }, classifiers=[