diff --git a/config/env/20241110-neurotoxin.toml b/config/env/20241110-neurotoxin.toml new file mode 100644 index 0000000..58bdd24 --- /dev/null +++ b/config/env/20241110-neurotoxin.toml @@ -0,0 +1,42 @@ +n_initial_agents = 100 +n_max_agents = 240 +n_max_foods = 160 +n_food_sources = 2 +observe_food_label = true +food_num_fn = [ + ["linear", 20, 0.2, 120], + ["linear", 20, 0.2, 40], +] +food_loc_fn = ["uniform", "uniform"] +food_color = [[254, 2, 162, 255], [2, 254, 162, 255]] +food_energy_coef = [1.0] +agent_loc_fn = "uniform" +xlim = [0.0, 480.0] +ylim = [0.0, 480.0] +env_shape = "square" +neighbor_stddev = 100.0 +n_agent_sensors = 24 +sensor_length = 200.0 +sensor_range = "wide" +agent_radius = 10.0 +food_radius = 4.0 +foodloc_interval = 1000 +dt = 0.1 +linear_damping = 0.8 +angular_damping = 0.6 +max_force = 80.0 +min_force = -20.0 +init_energy = 80.0 +energy_capacity = 400.0 +force_energy_consumption = 1e-5 +basic_energy_consumption = 5e-4 +energy_share_ratio = 0.4 +n_velocity_iter = 6 +n_position_iter = 2 +n_physics_iter = 5 +max_place_attempts = 10 +n_max_food_regen = 10 +toxin_t0 = 5.0 +toxin_alpha = 1.0 +toxin_decay = 0.01 +toxin_delta = 10.0 \ No newline at end of file diff --git a/experiments/cf_simple.py b/experiments/cf_simple.py index b23e2cf..c106ce7 100644 --- a/experiments/cf_simple.py +++ b/experiments/cf_simple.py @@ -499,7 +499,7 @@ def replay( videopath: Path | None = None, start: int = 0, end: int | None = None, - cfconfig_path: Path = PROJECT_ROOT / "config/env/20231214-square.toml", + cfconfig_path: Path = DEFAULT_CFCONFIG, env_override: str = "", ) -> None: with cfconfig_path.open("r") as f: diff --git a/experiments/cf_toxin.py b/experiments/cf_toxin.py new file mode 100644 index 0000000..12c0db1 --- /dev/null +++ b/experiments/cf_toxin.py @@ -0,0 +1,688 @@ +"""Asexual reward evolution with Circle Foraging""" + +import dataclasses +import itertools +import json +from pathlib import Path +from typing import cast + +import chex +import equinox as eqx +import jax +import jax.numpy as jnp +import numpy as np +import optax +import typer +from serde import serde, toml + +from emevo import Env +from emevo import birth_and_death as bd +from emevo import genetic_ops as gops +from emevo import make +from emevo import reward_fn as rfn +from emevo.env import ObsProtocol as Obs +from emevo.env import StateProtocol as State +from emevo.eqx_utils import get_slice +from emevo.eqx_utils import where as eqx_where +from emevo.exp_utils import ( + BDConfig, + CfConfig, + FoodLog, + GopsConfig, + Log, + Logger, + LogMode, + SavedPhysicsState, + SavedProfile, + is_cuda_ready, +) +from emevo.rl import ppo_normal as ppo +from emevo.spaces import BoxSpace +from emevo.visualizer import SaveVideoWrapper + +PROJECT_ROOT = Path(__file__).parent.parent +DEFAULT_CFCONFIG = PROJECT_ROOT / "config/env/20241110-neurotoxin.toml" + + +@serde +@dataclasses.dataclass +class CfConfigWithToxin(CfConfig): + toxin_t0: float = 5.0 + toxin_alpha: float = 1.0 + toxin_decay: float = 0.01 + toxin_delta: float = 10.0 + + +@dataclasses.dataclass +class RewardExtractor: + act_space: BoxSpace + act_coef: float + _max_norm: jax.Array = dataclasses.field(init=False) + + def __post_init__(self) -> None: + self._max_norm = jnp.sqrt(jnp.sum(self.act_space.high**2, axis=-1)) + + def normalize_action(self, action: jax.Array) -> jax.Array: + scaled = self.act_space.sigmoid_scale(action) + norm = jnp.sqrt(jnp.sum(scaled**2, axis=-1, keepdims=True)) + return norm / self._max_norm + + def extract( + self, + ate_food: jax.Array, + action: jax.Array, + energy: jax.Array, + ) -> jax.Array: + del energy + act_input = self.act_coef * self.normalize_action(action) + return jnp.concatenate((ate_food.astype(jnp.float32), act_input), axis=1) + + +def serialize_weight(w: jax.Array) -> dict[str, jax.Array]: + wd = w.shape[0] + rd = {f"food_{i + 1}": rfn.slice_last(w, i) for i in range(wd - 1)} + rd["action"] = rfn.slice_last(w, wd - 1) + return rd + + +def exec_rollout( + state: State, + initial_obs: Obs, + env: Env, + network: ppo.NormalPPONet, + reward_fn: rfn.RewardFn, + hazard_fn: bd.HazardFunction, + birth_fn: bd.BirthFunction, + prng_key: jax.Array, + n_rollout_steps: int, +) -> tuple[State, ppo.Rollout, Log, FoodLog, SavedPhysicsState, Obs, jax.Array]: + def step_rollout( + carried: tuple[State, Obs], + key: jax.Array, + ) -> tuple[tuple[State, Obs], tuple[ppo.Rollout, Log, FoodLog, SavedPhysicsState]]: + act_key, hazard_key, birth_key = jax.random.split(key, 3) + state_t, obs_t = carried + obs_t_array = obs_t.as_array() + net_out = ppo.vmap_apply(network, obs_t_array) + actions = net_out.policy().sample(seed=act_key) + state_t1, timestep = env.step( + state_t, + env.act_space.sigmoid_scale(actions), # type: ignore + ) + obs_t1 = timestep.obs + energy = state_t.status.energy + rewards = reward_fn(timestep.info["n_ate_food"], actions, energy).reshape(-1, 1) + rollout = ppo.Rollout( + observations=obs_t_array, + actions=actions, + rewards=rewards, + terminations=jnp.zeros_like(rewards), + values=net_out.value, + means=net_out.mean, + logstds=net_out.logstd, + ) + # Birth and death + death_prob = hazard_fn(state_t1.status.age, state_t1.status.energy) + dead_nonzero = jax.random.bernoulli(hazard_key, p=death_prob) + dead = jnp.where( + # If the agent's energy is lower than 0, it should immediately die + state_t1.status.energy < 0.0, + jnp.ones_like(dead_nonzero), + dead_nonzero, + ) + state_t1d = env.deactivate(state_t1, dead) + birth_prob = birth_fn(state_t1d.status.age, state_t1d.status.energy) + possible_parents = jnp.logical_and( + jnp.logical_and( + jnp.logical_not(dead), + state.unique_id.is_active(), # type: ignore + ), + jax.random.bernoulli(birth_key, p=birth_prob), + ) + state_t1db, parents = env.activate(state_t1d, possible_parents) + log = Log( + dead=jnp.where(dead, state_t.unique_id.unique_id, -1), # type: ignore + n_got_food=timestep.info["n_ate_food"], + action_magnitude=actions, + energy_gain=timestep.info["energy_gain"], + consumed_energy=timestep.info["energy_consumption"], + energy=state_t1db.status.energy, + parents=parents, + rewards=rewards.ravel(), + unique_id=state_t1db.unique_id.unique_id, + additional_fields={"n_got_toxin": timestep.info["n_ate_toxin"]}, + ) + foodlog = FoodLog( + eaten=timestep.info["n_food_eaten"], + regenerated=timestep.info["n_food_regenerated"], + ) + phys = state_t.physics # type: ignore + phys_state = SavedPhysicsState( + circle_axy=phys.circle.p.into_axy(), + static_circle_axy=phys.static_circle.p.into_axy(), + circle_is_active=phys.circle.is_active, + static_circle_is_active=phys.static_circle.is_active, + static_circle_label=phys.static_circle.label, + ) + return (state_t1db, obs_t1), (rollout, log, foodlog, phys_state) + + (state, obs), (rollout, log, foodlog, phys_state) = jax.lax.scan( + step_rollout, + (state, initial_obs), + jax.random.split(prng_key, n_rollout_steps), + ) + next_value = ppo.vmap_value(network, obs.as_array()) + return state, rollout, log, foodlog, phys_state, obs, next_value + + +@eqx.filter_jit +def epoch( + state: State, + initial_obs: Obs, + env: Env, + network: ppo.NormalPPONet, + reward_fn: rfn.RewardFn, + hazard_fn: bd.HazardFunction, + birth_fn: bd.BirthFunction, + prng_key: jax.Array, + n_rollout_steps: int, + gamma: float, + gae_lambda: float, + adam_update: optax.TransformUpdateFn, + opt_state: optax.OptState, + minibatch_size: int, + n_optim_epochs: int, + entropy_weight: float, +) -> tuple[ + State, Obs, Log, FoodLog, SavedPhysicsState, optax.OptState, ppo.NormalPPONet +]: + keys = jax.random.split(prng_key, env.n_max_agents + 1) + env_state, rollout, log, foodlog, phys_state, obs, next_value = exec_rollout( + state, + initial_obs, + env, + network, + reward_fn, + hazard_fn, + birth_fn, + keys[0], + n_rollout_steps, + ) + batch = ppo.vmap_batch(rollout, next_value, gamma, gae_lambda) + opt_state, pponet = ppo.vmap_update( + batch, + network, + adam_update, + opt_state, + keys[1:], + minibatch_size, + n_optim_epochs, + 0.2, + entropy_weight, + ) + return env_state, obs, log, foodlog, phys_state, opt_state, pponet + + +def run_evolution( + *, + key: jax.Array, + env: Env, + n_initial_agents: int, + adam: optax.GradientTransformation, + gamma: float, + gae_lambda: float, + n_optim_epochs: int, + minibatch_size: int, + n_rollout_steps: int, + n_total_steps: int, + entropy_weight: float, + reward_fn: rfn.RewardFn, + hazard_fn: bd.HazardFunction, + birth_fn: bd.BirthFunction, + mutation: gops.Mutation, + xmax: float, + ymax: float, + logger: Logger, + save_interval: int, + debug_vis: bool, +) -> None: + key, net_key, reset_key = jax.random.split(key, 3) + obs_space = env.obs_space.flatten() + input_size = int(np.prod(obs_space.shape)) + act_size = int(np.prod(env.act_space.shape)) + + def initialize_net(key: chex.PRNGKey) -> ppo.NormalPPONet: + return ppo.vmap_net( + input_size, + 64, + act_size, + jax.random.split(key, env.n_max_agents), + ) + + pponet = initialize_net(net_key) + adam_init, adam_update = adam + + @eqx.filter_jit + def initialize_opt_state(net: eqx.Module) -> optax.OptState: + return jax.vmap(adam_init)(eqx.filter(net, eqx.is_array)) + + @eqx.filter_jit + def replace_net( + key: chex.PRNGKey, + flag: jax.Array, + pponet: ppo.NormalPPONet, + opt_state: optax.OptState, + ) -> tuple[ppo.NormalPPONet, optax.OptState]: + initialized = initialize_net(key) + pponet = eqx_where(flag, initialized, pponet) + opt_state = jax.tree_util.tree_map( + lambda a, b: jnp.where( + jnp.expand_dims(flag, tuple(range(1, a.ndim))), + b, + a, + ), + opt_state, + initialize_opt_state(pponet), + ) + return pponet, opt_state + + opt_state = initialize_opt_state(pponet) + env_state, timestep = env.reset(reset_key) + obs = timestep.obs + + if debug_vis: + visualizer = env.visualizer(env_state, figsize=(xmax * 2, ymax * 2)) + else: + visualizer = None + + for i in range(n_initial_agents): + logger.reward_fn_dict[i + 1] = get_slice(reward_fn, i) + logger.profile_dict[i + 1] = SavedProfile(0, 0, i + 1) + + for i, key_i in enumerate(jax.random.split(key, n_total_steps // n_rollout_steps)): + epoch_key, init_key = jax.random.split(key_i) + old_state = env_state + # Use `with jax.disable_jit():` here for debugging + env_state, obs, log, foodlog, phys_state, opt_state, pponet = epoch( + env_state, + obs, + env, + pponet, + reward_fn, + hazard_fn, + birth_fn, + epoch_key, + n_rollout_steps, + gamma, + gae_lambda, + adam_update, + opt_state, + minibatch_size, + n_optim_epochs, + entropy_weight, + ) + print(jnp.sum(log.additional_fields["n_got_toxin"])) + + if visualizer is not None: + visualizer.render(env_state.physics) # type: ignore + visualizer.show() + is_active = env_state.unique_id.is_active() + popl = int(jnp.sum(is_active)) + avg_e = float(jnp.mean(env_state.status.energy[is_active])) + if popl > 0: + print(f"Population: {popl} Avg. Energy: {avg_e}") + + # Extinct? + n_active = jnp.sum(env_state.unique_id.is_active()) # type: ignore + if n_active == 0: + print(f"Extinct after {i + 1} epochs") + break + + # Save dead agents + log_with_step = log.with_step(i * n_rollout_steps) + log_death = log_with_step.filter_death() + ages = old_state.status.age[log_death.slots] + logger.save_agents( + pponet, + log_death.log.dead, + log_death.slots, + ages + log_death.step - i * n_rollout_steps, + ) + # Save alive agents + saved = jnp.logical_and( + env_state.status.age > 0, + ((env_state.status.age // n_rollout_steps) % save_interval) == 0, + ) + (saved_slots,) = jnp.nonzero(saved) + logger.save_agents( + pponet, + env_state.unique_id.unique_id[saved_slots], + saved_slots, + env_state.status.age[saved_slots], + prefix="intermediate", + ) + # Initialize network and adam state for new agents + log_birth = log_with_step.filter_birth() + is_new = jnp.zeros(env.n_max_agents, dtype=bool).at[log_birth.slots].set(True) + if jnp.any(is_new): + pponet, opt_state = replace_net(init_key, is_new, pponet, opt_state) + + # Mutation + reward_fn = rfn.mutate_reward_fn( + key, + logger.reward_fn_dict, + reward_fn, + mutation, + log_birth.log.parents, + log_birth.log.unique_id, + log_birth.slots, + ) + # Update profile + for step, uid, parent in zip( + log_birth.step, + log_birth.log.unique_id, + log_birth.log.parents, + ): + ui = uid.item() + logger.profile_dict[ui] = SavedProfile(step.item(), parent.item(), ui) + + # Push log and physics state + logger.push_foodlog(foodlog) + logger.push_log(log_with_step.filter_active()) + logger.push_physstate(phys_state) + + # Save logs before exiting + logger.finalize() + is_active = env_state.unique_id.is_active() + logger.save_agents( + pponet, + env_state.unique_id.unique_id[is_active], + jnp.arange(len(is_active))[is_active], + env_state.status.age[is_active], + ) + + +app = typer.Typer(pretty_exceptions_show_locals=False) + + +@app.command() +def evolve( + seed: int = 1, + adam_lr: float = 3e-4, + adam_eps: float = 1e-7, + gamma: float = 0.999, + gae_lambda: float = 0.95, + n_optim_epochs: int = 10, + minibatch_size: int = 256, + n_rollout_steps: int = 1024, + n_total_steps: int = 1024 * 10000, + act_reward_coef: float = 0.01, + entropy_weight: float = 0.001, + cfconfig_path: Path = DEFAULT_CFCONFIG, + bdconfig_path: Path = PROJECT_ROOT / "config/bd/20240318-mild-slope.toml", + gopsconfig_path: Path = PROJECT_ROOT / "config/gops/20240326-cauthy-002.toml", + min_age_for_save: int = 0, + save_interval: int = 100000000, # No saving by default + env_override: str = "", + birth_override: str = "", + hazard_override: str = "", + gops_params_override: str = "", + logdir: Path = Path("./log"), + log_mode: LogMode = LogMode.REWARD_LOG_STATE, + log_interval: int = 1000, + savestate_interval: int = 1000, + debug_vis: bool = False, + force_gpu: bool = True, +) -> None: + if force_gpu and not is_cuda_ready(): + raise RuntimeError("Detected some problem in CUDA!") + + # Load config + with cfconfig_path.open("r") as f: + cfconfig = toml.from_toml(CfConfigWithToxin, f.read()) + with bdconfig_path.open("r") as f: + bdconfig = toml.from_toml(BDConfig, f.read()) + with gopsconfig_path.open("r") as f: + gopsconfig = toml.from_toml(GopsConfig, f.read()) + + # Apply overrides + cfconfig.apply_override(env_override) + bdconfig.apply_birth_override(birth_override) + bdconfig.apply_hazard_override(hazard_override) + gopsconfig.apply_params_override(gops_params_override) + + # Load models + birth_fn, hazard_fn = bdconfig.load_models() + mutation = gopsconfig.load_model() + # Make env + env = make("CircleForaging-v1", **dataclasses.asdict(cfconfig)) + key, reward_key = jax.random.split(jax.random.PRNGKey(seed)) + reward_extracor = RewardExtractor( + act_space=env.act_space, # type: ignore + act_coef=act_reward_coef, + ) + reward_fn_instance = rfn.LinearReward( + key=reward_key, + n_agents=cfconfig.n_max_agents, + n_weights=cfconfig.n_food_sources, # Because one of the foods is toxin + std=gopsconfig.init_std, + mean=gopsconfig.init_mean, + extractor=reward_extracor.extract, + serializer=serialize_weight, + **gopsconfig.init_kwargs, + ) + + logger = Logger( + logdir=logdir, + mode=log_mode, + log_interval=log_interval, + savestate_interval=savestate_interval, + min_age_for_save=min_age_for_save, + ) + run_evolution( + key=key, + env=env, + n_initial_agents=cfconfig.n_initial_agents, + adam=optax.adam(adam_lr, eps=adam_eps), + gamma=gamma, + gae_lambda=gae_lambda, + n_optim_epochs=n_optim_epochs, + minibatch_size=minibatch_size, + n_rollout_steps=n_rollout_steps, + n_total_steps=n_total_steps, + entropy_weight=entropy_weight, + reward_fn=reward_fn_instance, + hazard_fn=hazard_fn, + birth_fn=birth_fn, + mutation=cast(gops.Mutation, mutation), + xmax=cfconfig.xlim[1], + ymax=cfconfig.ylim[1], + logger=logger, + save_interval=save_interval, + debug_vis=debug_vis, + ) + + +@app.command() +def replay( + physstate_path: Path, + backend: str = "pyglet", # Use "headless" for headless rendering + videopath: Path | None = None, + start: int = 0, + end: int | None = None, + cfconfig_path: Path = DEFAULT_CFCONFIG, + env_override: str = "", +) -> None: + with cfconfig_path.open("r") as f: + cfconfig = toml.from_toml(CfConfig, f.read()) + # For speedup + cfconfig.n_initial_agents = 1 + cfconfig.apply_override(env_override) + phys_state = SavedPhysicsState.load(physstate_path) + env = make("CircleForaging-v1", **dataclasses.asdict(cfconfig)) + env_state, _ = env.reset(jax.random.PRNGKey(0)) + end_index = end if end is not None else phys_state.circle_axy.shape[0] + visualizer = env.visualizer( + env_state, + figsize=(cfconfig.xlim[1] * 2, cfconfig.ylim[1] * 2), + backend=backend, + ) + if videopath is not None: + visualizer = SaveVideoWrapper(visualizer, videopath, fps=60) + for i in range(start, end_index): + phys = phys_state.set_by_index(i, env_state.physics) + env_state = dataclasses.replace(env_state, physics=phys) + visualizer.render(env_state.physics) + visualizer.show() + visualizer.close() + + +@app.command() +def widget( + physstate_path: Path, + start: int = 0, + end: int | None = None, + cfconfig_path: Path = DEFAULT_CFCONFIG, + log_path: Path | None = None, + self_terminate: bool = False, + profile_and_rewards_path: Path | None = None, + cm_fixed_minmax: str = "", + env_override: str = "", + scale: float = 2.0, + force_cpu: bool = False, +) -> None: + from emevo.analysis.qt_widget import CFEnvReplayWidget, start_widget + + if force_cpu: + jax.config.update("jax_default_device", jax.devices("cpu")[0]) + + with cfconfig_path.open("r") as f: + cfconfig = toml.from_toml(CfConfigWithToxin, f.read()) + + # For speedup + cfconfig.n_initial_agents = 1 + cfconfig.apply_override(env_override) + phys_state = SavedPhysicsState.load(physstate_path) + env = make("CircleForaging-v1", **dataclasses.asdict(cfconfig)) + end = phys_state.circle_axy.shape[0] if end is None else end + if log_path is None: + log_ds = None + step_offset = 0 + else: + import pyarrow.dataset as ds + + log_ds = ds.dataset(log_path) + step_offset = log_ds.scanner(columns=["step"]).head(1)["step"][0].as_py() + + if profile_and_rewards_path is None: + profile_and_rewards = None + else: + import pyarrow.parquet as pq + + profile_and_rewards = pq.read_table(profile_and_rewards_path) + + if len(cm_fixed_minmax) > 0: + cm_fixed_minmax_dict = json.loads(cm_fixed_minmax) + else: + cm_fixed_minmax_dict = {} + + start_widget( + CFEnvReplayWidget, + xlim=int(cfconfig.xlim[1]), + ylim=int(cfconfig.ylim[1]), + env=env, + saved_physics=phys_state, + start=start, + end=end, + log_ds=log_ds, + step_offset=step_offset, + self_terminate=self_terminate, + profile_and_rewards=profile_and_rewards, + cm_fixed_minmax=cm_fixed_minmax_dict, + scale=scale, + ) + + +@app.command() +def vis_policy( + physstate_path: Path, + policy_path: list[Path], + subtitle: list[str] | None = None, + agent_index: int | None = None, + cfconfig_path: Path = DEFAULT_CFCONFIG, + fig_unit: float = 4.0, + scale: float = 1.0, +) -> None: + from emevo.analysis.policy import draw_cf_policy + + with cfconfig_path.open("r") as f: + cfconfig = toml.from_toml(CfConfigWithToxin, f.read()) + + cfconfig.n_initial_agents = 1 + # Load env state + phys_state = SavedPhysicsState.load(physstate_path) + env = make("CircleForaging-v1", **dataclasses.asdict(cfconfig)) + key = jax.random.PRNGKey(0) + env_state, _ = env.reset(key) + loaded_phys = phys_state.set_by_index(..., env_state.physics) + env_state = dataclasses.replace(env_state, physics=loaded_phys) + # agent_index + if agent_index is None: + file_name = physstate_path.stem + if "slot" in file_name: + agent_index = int(file_name[file_name.index("slot") + 4 :]) + else: + print("Set --agent-index") + return + # Load agents + input_size = int(np.prod(env.obs_space.flatten().shape)) + act_size = int(np.prod(env.act_space.shape)) + ref_net = ppo.NormalPPONet(input_size, 64, act_size, key) + names, net_params = [], [] + for policy_path_i, name in itertools.zip_longest( + policy_path, + [] if subtitle is None else subtitle, + ): + pponet = eqx.tree_deserialise_leaves(policy_path_i, ref_net) + # Append only params of the network, excluding functions (etc. tanh). + net_params.append(eqx.filter(pponet, eqx.is_array)) + names.append(policy_path_i.stem if name is None else name) + net_params = jax.tree.map(lambda *args: jnp.stack(args), *net_params) + network = eqx.combine(net_params, ref_net) + # Get obs + n_agents = cfconfig.n_max_agents + zero_action = jnp.zeros((n_agents, *env.act_space.shape)) + _, timestep = env.step(env_state, zero_action) + obs_array = timestep.obs.as_array() + obs_i = obs_array[agent_index] + + @eqx.filter_vmap(in_axes=(eqx.if_array(0), None)) + def evaluate(network: ppo.NormalPPONet, obs: jax.Array) -> ppo.Output: + return network(obs) + + # Get output + output = evaluate(network, obs_i) + # Make visualizer + visualizer = env.visualizer( + env_state, + figsize=(cfconfig.xlim[1] * scale, cfconfig.ylim[1] * scale), + sensor_index=agent_index, + sensor_width=0.004, + sensor_color=np.array([0.0, 0.0, 0.0, 0.3], dtype=np.float32), + ) + visualizer.render(env_state.physics) + visualizer.show() + max_force = max(cfconfig.max_force, -cfconfig.min_force) + rot = env_state.physics.circle.p.angle[agent_index].item() + policy_mean = env.act_space.sigmoid_scale(output.mean) + draw_cf_policy( + names, + np.array(policy_mean), + rotation=rot, + fig_unit=fig_unit, + max_force=max_force, + ) + + +if __name__ == "__main__": + app() diff --git a/notebooks/toxin.ipynb b/notebooks/toxin.ipynb new file mode 100644 index 0000000..7079412 --- /dev/null +++ b/notebooks/toxin.ipynb @@ -0,0 +1,240 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "0dadbf8d-d3eb-42b8-a9c1-265e61f6edd8", + "metadata": {}, + "outputs": [], + "source": [ + "import dataclasses\n", + "from typing import Any, Literal\n", + "\n", + "import ipywidgets as widgets\n", + "import numpy as np\n", + "from emevo import birth_and_death as bd\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.figure import Figure\n", + "from matplotlib.lines import Line2D\n", + "from matplotlib.text import Text\n", + "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", + "\n", + "from emevo.plotting import (\n", + " vis_birth,\n", + " vis_expected_n_children,\n", + " vis_hazard,\n", + " vis_lifetime,\n", + " show_params_text,\n", + ")\n", + "\n", + "%matplotlib ipympl" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c4b6ebde-ea34-4a3e-92b3-964ac39c8452", + "metadata": {}, + "outputs": [], + "source": [ + "def make_slider(\n", + " vmin: float,\n", + " vmax: float,\n", + " logscale: bool = True,\n", + " n_steps: int = 400,\n", + ") -> widgets.FloatSlider | widgets.FloatLogSlider:\n", + " if logscale:\n", + " logmin = np.log10(vmin)\n", + " logmax = np.log10(vmax)\n", + " logstep = (logmax - logmin) / n_steps\n", + " return widgets.FloatLogSlider(\n", + " min=logmin,\n", + " max=logmax,\n", + " step=logstep,\n", + " value=10 ** ((logmax + logmin) / 2.0),\n", + " base=10,\n", + " readout_format=\".3e\",\n", + " )\n", + " else:\n", + " str_vmin = str(min(abs(vmin), abs(vmax)))\n", + " dot = str_vmin.find(\".\")\n", + " if dot != -1:\n", + " format_n = len(str_vmin[dot + 1: ].rstrip(\"0\"))\n", + " readout_format=f\".{format_n + 1}e\"\n", + " else:\n", + " readout_format = \".2f\"\n", + " return widgets.FloatSlider(\n", + " min=vmin,\n", + " max=vmax,\n", + " step=(vmax - vmin) / n_steps,\n", + " value=(vmax + vmin) / 2,\n", + " readout_format=readout_format,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "8d94990f-99ac-4dc0-a5c9-1b19806b8885", + "metadata": {}, + "outputs": [], + "source": [ + "def savefig_widgets(fig: Figure) -> list:\n", + " text = widgets.Text(\n", + " value=\"figure.png\",\n", + " description=\"Filename:\",\n", + " disabled=False,\n", + " )\n", + " button = widgets.Button(description=\"Save File\")\n", + " output = widgets.Output()\n", + "\n", + " def on_button_clicked(b):\n", + " filename = text.value\n", + " if any([filename.endswith(ext) for ext in [\".png\", \".svg\", \".pdf\"]]):\n", + " fig.savefig(filename)\n", + " else:\n", + " with output:\n", + " print(\"Enter valid file name!\")\n", + " \n", + " button.on_click(on_button_clicked)\n", + " return [text, button, output]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "8c0388cd-0f78-4094-8129-a448bd2446e0", + "metadata": {}, + "outputs": [], + "source": [ + "def make_toxin_widget(**kwargs) -> widgets.VBox:\n", + " fig = plt.figure(figsize=(8, 6))\n", + " ax = fig.add_subplot(111)\n", + " ax.set_title(\"Motor output decreased by toxin\")\n", + " \n", + " @dataclasses.dataclass\n", + " class State:\n", + " line: Line2D | None = None\n", + "\n", + " state = State()\n", + "\n", + " x = np.linspace(0, 10, 1000)\n", + " def update_figure(alpha, t0):\n", + " y = 1.0 / (1.0 + alpha * np.exp(t0 - x))\n", + " if state.line is None:\n", + " ax.grid(True, which=\"major\")\n", + " ax.set_xlabel(\"Energy\", fontsize=12)\n", + " ax.set_ylabel(\"Decrease Ratio\", fontsize=12)\n", + " else:\n", + " state.line.remove()\n", + " \n", + " state.line = ax.plot(x, y, color=\"xkcd:bluish purple\")[0]\n", + " fig.canvas.draw()\n", + " fig.canvas.flush_events()\n", + "\n", + " sliders = {key: make_slider(*range_) for key, range_ in kwargs.items()}\n", + " interactive = widgets.interactive(update_figure, **sliders)\n", + " return widgets.VBox(savefig_widgets(fig) + [interactive])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6db2c2eb-c63d-439e-8030-e374dfeb3d5f", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "07346a69d84645ca8825c414e26e52c4", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Text(value='figure.png', description='Filename:'), Button(description='Save File', style=Button…" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "9a3ebfe213f64499848c0a6839d70f5a", + "version_major": 2, + "version_minor": 0 + }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrXElEQVR4nO3dd3hUZd7G8Xtmkkw6nQChhN4lSBMEBKQoiGLDTlHRFdgVsCy6AqICr6KIuqygIijKLlZERSEioEivgnRpAgYIJYHUycx5/xgyEpJAgGROTvL9XFcucsrM+U2eJMydpxybYRiGAAAAAMAP7GYXAAAAAKDkIIAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAAAA8BsCCAAAAAC/IYAAAC5o5syZstls2rdvn9mlWMKAAQMUExNz0fNiYmJ00003FX5BfrBv3z7ZbDbNnDnT7FIAWAABBECByHqTarPZtGzZshzHDcNQtWrVZLPZLvtN1/jx4zV37twrrNRc8+fP1/PPP++Xa6WkpOj555/XkiVL/HI9FD2HDx/W888/r40bN5pdCgD4EEAAFKjg4GDNnj07x/6lS5fq4MGDcjqdl/3cxSWAjB071i/XSklJ0dixYwkgJdjhw4c1duzYQg8gNWrUUGpqqh544IFCvQ6A4oEAAqBA9ezZU59++qkyMzOz7Z89e7ZatGihSpUqmVRZ7pKTk80uAbmgXazFZrMpODhYDofD7FIAWAABBECBuueee3T8+HHFxcX59mVkZOizzz7Tvffem+tjkpOT9cQTT6hatWpyOp2qX7++Xn31VRmG4TvHZrMpOTlZH3zwgW+o14ABA3zHN2zYoBtvvFGRkZEKDw/X9ddfr5UrV2a7TtYwsaVLl2rw4MGqWLGiqlatesHXc/ToUT300EOKiopScHCwmjVrpg8++CDbOUuWLJHNZsvR03D+uPgBAwZoypQpvteT9XHuua+++qpef/111ahRQyEhIbruuuu0ZcuWbM/bqVMnderUKUet58492LdvnypUqCBJGjt2rO9aFxv+9dtvv6lLly4KCQlR1apV9dJLL8nj8eR67nfffacOHTooLCxMERER6tWrl3777bcc523fvl19+/ZVhQoVFBISovr16+tf//qX7/jzzz8vm82mrVu36t5771WZMmXUvn173/GPPvpILVq0UEhIiMqWLau7775bf/zxR7Zr/Pzzz7rzzjtVvXp1OZ1OVatWTcOHD1dqamq28+Lj4zVw4EBVrVpVTqdTlStX1i233JJjfkt+X9vcuXPVpEkTBQcHq0mTJvryyy8v+PXNzcKFCxUbG6vg4GA1atRIX3zxhe/Ynj17ZLPZ9Prrr+d43PLly2Wz2fTf//431+ddsmSJWrVqJUkaOHCg73vg3Hkan376qe9rW758ed1///06dOiQ7/iYMWNkt9u1aNGibM/9yCOPKCgoSJs2bZKU+xyQAQMGKDw8XIcOHVKfPn0UHh6uChUq6Mknn5Tb7b7krxOA4iPA7AIAFC8xMTFq27at/vvf/+rGG2+U5H0zl5iYqLvvvltvvvlmtvMNw9DNN9+sxYsX66GHHlJsbKwWLFigp556SocOHfK98Zo1a5YefvhhtW7dWo888ogkqXbt2pK8b5o7dOigyMhIPf300woMDNS0adPUqVMnLV26VG3atMl2zcGDB6tChQoaPXr0Bf/Snpqaqk6dOmn37t0aOnSoatasqU8//VQDBgzQqVOn9Pjjj1/S1+bRRx/V4cOHFRcXp1mzZuV6zocffqjTp09ryJAhSktL0xtvvKEuXbpo8+bNioqKyve1KlSooLfffluPPfaYbr31Vt12222SpKuuuirPx8THx6tz587KzMzUyJEjFRYWpnfeeUchISE5zp01a5b69++vHj166OWXX1ZKSorefvtttW/fXhs2bPAFoV9//VUdOnRQYGCgHnnkEcXExOj333/X119/rXHjxmV7zjvvvFN169bV+PHjfeFz3LhxGjVqlPr27auHH35Yx44d01tvvaWOHTtqw4YNKl26tCTvG+mUlBQ99thjKleunFavXq233npLBw8e1Keffuq7xu23367ffvtNf//73xUTE6OjR48qLi5OBw4c8NWc39e2cOFC3X777WrUqJEmTJig48eP+8JNfu3atUt33XWX/va3v6l///6aMWOG7rzzTn3//ffq1q2batWqpWuvvVYff/yxhg8fnu2xH3/8sSIiInTLLbfk+twNGzbUCy+8oNGjR+uRRx5Rhw4dJEnt2rWT5A3kAwcOVKtWrTRhwgQdOXJEb7zxhn755Rff1/a5557T119/rYceekibN29WRESEFixYoHfffVcvvviimjVrdsHX53a71aNHD7Vp00avvvqqfvjhB7322muqXbu2HnvssXx/nQAUMwYAFIAZM2YYkow1a9YY//73v42IiAgjJSXFMAzDuPPOO43OnTsbhmEYNWrUMHr16uV73Ny5cw1JxksvvZTt+e644w7DZrMZu3fv9u0LCwsz+vfvn+Paffr0MYKCgozff//dt+/w4cNGRESE0bFjxxw1tm/f3sjMzLzoa5o8ebIhyfjoo498+zIyMoy2bdsa4eHhRlJSkmEYhrF48WJDkrF48eJsj9+7d68hyZgxY4Zv35AhQ4zcfvVmnRsSEmIcPHjQt3/VqlWGJGP48OG+fdddd51x3XXX5XiO/v37GzVq1PBtHzt2zJBkjBkz5qKv1TAMY9iwYYYkY9WqVb59R48eNUqVKmVIMvbu3WsYhmGcPn3aKF26tDFo0KBsj4+PjzdKlSqVbX/Hjh2NiIgIY//+/dnO9Xg8vs/HjBljSDLuueeebOfs27fPcDgcxrhx47Lt37x5sxEQEJBtf9b32rkmTJhg2Gw237VPnjxpSDImTpyY59fgUl5bbGysUblyZePUqVO+fQsXLjQkZWuHvNSoUcOQZHz++ee+fYmJiUblypWN5s2b+/ZNmzbNkGRs27bNty8jI8MoX758rj8P51qzZk2O78Gsx1esWNFo0qSJkZqa6tv/zTffGJKM0aNH+/Zt3rzZCAoKMh5++GHj5MmTRnR0tNGyZUvD5XL5zsnte71///6GJOOFF17Idu3mzZsbLVq0uGDdAIo3hmABKHB9+/ZVamqqvvnmG50+fVrffPNNnsOv5s+fL4fDoX/84x/Z9j/xxBMyDEPffffdBa/ldru1cOFC9enTR7Vq1fLtr1y5su69914tW7ZMSUlJ2R4zaNCgfI1Vnz9/vipVqqR77rnHty8wMFD/+Mc/dObMGS1duvSiz3Gp+vTpo+joaN9269at1aZNG82fP7/Ar3W++fPn65prrlHr1q19+ypUqKD77rsv23lxcXE6deqU7rnnHiUkJPg+HA6H2rRpo8WLF0uSjh07pp9++kkPPvigqlevnu05soaenetvf/tbtu0vvvhCHo9Hffv2zXadSpUqqW7dur7rSMrWS5OcnKyEhAS1a9dOhmFow4YNvnOCgoK0ZMkSnTx5MtevQX5f259//qmNGzeqf//+KlWqlO/x3bp1U6NGjfL+Ip+nSpUquvXWW33bkZGR6tevnzZs2KD4+HhJ3p+n4OBgffzxx77zFixYoISEBN1///35vta51q5dq6NHj2rw4MEKDg727e/Vq5caNGigb7/91revSZMmGjt2rN577z316NFDCQkJ+uCDDxQQkL9BFOe3a4cOHbRnz57LqhtA8cAQLAAFrkKFCuratatmz56tlJQUud1u3XHHHbmeu3//flWpUkURERHZ9jds2NB3/EKOHTumlJQU1a9fP8exhg0byuPx6I8//lDjxo19+2vWrJmv17F//37VrVtXdnv2v9Xkt7bLUbdu3Rz76tWrp08++aTAr3W+/fv35xiuJinH13bXrl2SpC5duuT6PJGRkZLke5PZpEmTfF3//HbZtWuXDMPI9WsiecNglgMHDmj06NGaN29ejnCRmJgoSXI6nXr55Zf1xBNPKCoqStdcc41uuukm9evXz7c4Qn5fW1bb51Zb/fr1tX79+ou+XkmqU6dOjjBWr149Sd55FZUqVVLp0qXVu3dvzZ49Wy+++KIk7/Cr6OjoPOu8mKz6c/u5adCgQY6ltJ966in973//0+rVqzV+/Ph8h6zg4GDfXKQsZcqUyTMAAigZCCAACsW9996rQYMGKT4+XjfeeKNvrH5RkNuchiuR21/zJRXaRFubzZZtgn5hX+98WZPSZ82aleuqZvn9y/j5zm8Xj8cjm82m7777Ltceq/DwcEne192tWzedOHFC//znP9WgQQOFhYXp0KFDGjBgQLZJ9MOGDVPv3r01d+5cLViwQKNGjdKECRP0448/qnnz5oX22q5Uv3799Omnn2r58uVq2rSp5s2bp8GDB+cIx4Vlz549vnC2efPmfD+OVbEA5IYAAqBQ3HrrrXr00Ue1cuVKzZkzJ8/zatSooR9++EGnT5/O1guyfft23/Esub3Rr1ChgkJDQ7Vjx44cx7Zv3y673a5q1apd1muoUaOGfv31V3k8nmxv9M6vrUyZMpKkU6dOZXt8bj0keYWVLFlv8s61c+fObHfWLlOmTK5DWM6/3sWudb4aNWrkev3zv7ZZk/8rVqyorl275vl8WUPizl/FK79q164twzBUs2ZNX69AbjZv3qydO3fqgw8+UL9+/Xz7z12J7fznfeKJJ/TEE09o165dio2N1WuvvaaPPvoo368tq+3z8/W6kN27d8swjGxttXPnTknK1uY33HCDKlSooI8//lht2rRRSkpKvu65kdf3QFb9O3bsyNGLsmPHjmw/dx6PRwMGDFBkZKSGDRum8ePH64477vAtbAAAl4o5IAAKRXh4uN5++209//zz6t27d57n9ezZU263W//+97+z7X/99ddls9l8K2lJUlhYWI43+Q6HQ927d9dXX32VbSnVI0eOaPbs2Wrfvr1v2Myl6tmzp+Lj47MFqMzMTL311lsKDw/XddddJ8n7Zs7hcOinn37K9vj//Oc/OZ4zLCxMUs6wkmXu3LnZlkFdvXq1Vq1ale3rULt2bW3fvl3Hjh3z7du0aZN++eWXbM8VGhp6wWudr2fPnlq5cqVWr17t23fs2LFscw8kqUePHoqMjNT48ePlcrlyPE9WXRUqVFDHjh31/vvv68CBA9nOya0H53y33XabHA6Hxo4dm+N8wzB0/PhxSX/9lf3ccwzD0BtvvJHtMSkpKUpLS8u2r3bt2oqIiFB6evolvbbKlSsrNjZWH3zwgW+Il+QNPVu3br3oa8ty+PDhbEv3JiUl6cMPP1RsbGy2HpiAgADdc889+uSTTzRz5kw1bdr0giuaZcnr+61ly5aqWLGipk6d6nvtknfFum3btqlXr16+fZMmTdLy5cv1zjvv6MUXX1S7du302GOPKSEhId+vEwDORQ8IgELTv3//i57Tu3dvde7cWf/617+0b98+NWvWTAsXLtRXX32lYcOG+f4iLUktWrTQDz/8oEmTJqlKlSqqWbOm2rRpo5deeklxcXFq3769Bg8erICAAE2bNk3p6el65ZVXLrv+Rx55RNOmTdOAAQO0bt06xcTE6LPPPtMvv/yiyZMn+3psSpUqpTvvvFNvvfWWbDabateurW+++UZHjx7N8ZwtWrSQJP3jH/9Qjx495HA4dPfdd/uO16lTR+3bt9djjz2m9PR0TZ48WeXKldPTTz/tO+fBBx/UpEmT1KNHDz300EM6evSopk6dqsaNG2ebcB8SEqJGjRppzpw5qlevnsqWLasmTZrkOSfj6aef1qxZs3TDDTfo8ccf9y3Dm9UTlCUyMlJvv/22HnjgAV199dW6++67VaFCBR04cEDffvutrr32Wl+gfPPNN9W+fXtdffXVeuSRR1SzZk3t27dP33777UXvzl27dm299NJLeuaZZ7Rv3z716dNHERER2rt3r7788ks98sgjevLJJ9WgQQPVrl1bTz75pA4dOqTIyEh9/vnnOeYZ7Ny5U9dff7369u2rRo0aKSAgQF9++aWOHDnia4NLeW0TJkxQr1691L59ez344IM6ceKE3nrrLTVu3Fhnzpy54GvLUq9ePT300ENas2aNoqKi9P777+vIkSOaMWNGjnP79eunN998U4sXL9bLL7+cr+evXbu2SpcuralTpyoiIkJhYWFq06aNatasqZdfflkDBw7Uddddp3vuuce3DG9MTIxvyd9t27Zp1KhRGjBggO8PCTNnzlRsbKwGDx7sl7lJAIohk1bfAlDMnLsM74WcvwyvYXiXPh0+fLhRpUoVIzAw0Khbt64xceLEbEu1GoZhbN++3ejYsaMREhJiSMq2BOn69euNHj16GOHh4UZoaKjRuXNnY/ny5ZdV47mOHDliDBw40ChfvrwRFBRkNG3aNMeSpobhXfL29ttvN0JDQ40yZcoYjz76qLFly5YcS5NmZmYaf//7340KFSoYNpvNtyRv1jKmEydONF577TWjWrVqhtPpNDp06GBs2rQpx/U++ugjo1atWkZQUJARGxtrLFiwIMcyvIZhGMuXLzdatGhhBAUF5WtJ3l9//dW47rrrjODgYCM6Otp48cUXjenTp2dbhjfL4sWLjR49ehilSpUygoODjdq1axsDBgww1q5dm+28LVu2GLfeeqtRunRpIzg42Khfv74xatQo3/GsZXiPHTuWa02ff/650b59eyMsLMwICwszGjRoYAwZMsTYsWOH75ytW7caXbt2NcLDw43y5csbgwYNMjZt2pTt65+QkGAMGTLEaNCggREWFmaUKlXKaNOmjfHJJ5/kuGZ+X9vnn39uNGzY0HA6nUajRo2ML774Itd2yE3Wz8KCBQuMq666ynA6nUaDBg2MTz/9NM/HNG7c2LDb7dmWar6Yr776ymjUqJEREBCQ4/txzpw5RvPmzQ2n02mULVvWuO+++3zPnZmZabRq1cqoWrVqtqWGDcMw3njjDUOSMWfOHMMw8l6GNywsLEc9We0NoOSyGUY++sEBAIVq3759qlmzpiZOnKgnn3zS7HJQRDVv3lxly5bNcWdyALAS5oAAAGABa9eu1caNG7NNtAcAK2IOCAAARdiWLVu0bt06vfbaa6pcubLuuusus0sCgCtCDwgAAEXYZ599poEDB8rlcum///1vtjuXA4AVMQcEAAAAgN/QAwIAAADAbwggAAAAAPyGSehFiMfj0eHDhxURESGbzWZ2OQAAADiPYRg6ffq0qlSpIrudv+VfDgJIEXL48GFVq1bN7DIAAABwEX/88YeqVq1qdhmWRAApQiIiIiR5v6EjIyML/Xoul0sLFy5U9+7dFRgYWOjXQ8GjDa2PNrQ22s/6aEPr83cbJiUlqVq1ar73bbh0BJAiJGvYVWRkpN8CSGhoqCIjI/mla1G0ofXRhtZG+1kfbWh9ZrUhw+UvHwPXAAAAAPgNAQQAAACA3xBAAAAAAPgNAQQAAACA3xBAAAAAAPgNAQQAAACA3xBAAAAAAPgNAQQAAACA3xBAAAAAAPgNAQQAAACA3xBAAAAAAPgNAQQAAACA3xBAAAAAAPgNASQPP/30k3r37q0qVarIZrNp7ty5F33MkiVLdPXVV8vpdKpOnTqaOXNmodcJAAAAWAkBJA/Jyclq1qyZpkyZkq/z9+7dq169eqlz587auHGjhg0bpocfflgLFiwo5EoBAAAA6wgwu4Ci6sYbb9SNN96Y7/OnTp2qmjVr6rXXXpMkNWzYUMuWLdPrr7+uHj16FFaZAAAAgKUQQArIihUr1LVr12z7evTooWHDhplTEAAAgB8ZhnHO5+ceUO77L+HcbA/Ldo6hTJcht4tBPVZCACkg8fHxioqKyrYvKipKSUlJSk1NVUhISI7HpKenKz093bedlJQkSXK5XHK5XIVb8NnrnPsvrIc2tD7a0NqKavsZhiGPW/K4Jbdb8rgNGZ6s7b+Oec4eO/9zd27HPDkfl3We4Tn7/Ib3TaPhPvuvx5BhSB6PZHiy9nnrMzw593vOnp9931/P49vO7XnOe5wMQ4bO2T77uYy/3thm1XgqsbF2zTspm2x/naNzzr3g9tnrePL7OCPXOi66rexv8JX7p3m/2Tcu/dy8zs8tRJgtuHRT3dDTPz+HRe3n3YoIICaaMGGCxo4dm2P/woULFRoa6rc64uLi/HYtFA7a0PpoQ+syDGnBdz/InWmX59wPt00el9273312n+vsft853m3DbZPHY5ORbdvu/dd9dr/Hluu24bHJ47af8xw2GR7+GnxpwnX6T4/ZReAK+ev3aEpKil+uU5wRQApIpUqVdOTIkWz7jhw5osjIyFx7PyTpmWee0YgRI3zbSUlJqlatmrp3767IyMhCrVfyJvi4uDh169ZNgYGBhX49FDza0PpoQ/9xZxpKSzaUnmooI/WvfzNS9de+tKx955yTpuznpxlyZRjKzJBc6YZc6YYkm9kvL19sdsnu8H44HDbZHX/ty9r2fpzzuT37tiOXYza793ObXbLZJJvddt6291+73ZZzn0Oy2bLvt5/zPDabZHOcPdf3OJtvn81+zv6s822Szv6bbVt/7ZfNJpsktydTGzduVPPmsQoICMj2OCn7c0m287bPeS3ZnvvcbVuOa+f+3Odu23Ic932HnfOtZjvn82z7lcd+28XOteX2afbrXOzaefwo5HZ+/s7No6ZzTnFluvTjj6v99ns0a8QKLh8BpIC0bdtW8+fPz7YvLi5Obdu2zfMxTqdTTqczx/7AwEC/vhHx9/VQ8GhD66MNL86daSjltEcpiYbOJHqUetqjtGRDaWc8Sj1jKPWcz9POeLzbyR6lnvYGj4y0who3kv0NUkCQFOi0eT+CvP/mtS8g0KaAIJscAd7PHQGSI9CmgACbHIFZn3v/9Z0TmMv5gZIj4Oy/gTY5HN5jdodkc9jkyAoadmsEJX9yuVw6eOaUmnYI5WfQolwumxyBht9+j/J9cuUIIHk4c+aMdu/e7dveu3evNm7cqLJly6p69ep65plndOjQIX344YeSpL/97W/697//raeffloPPvigfvzxR33yySf69ttvzXoJAFBkGYY3FJw+4VHScY+SEz1KPuVRctLZzxM9Sk40lJJ4dt8pb7AoCEHBNu9HiE3OUJucId5tZ4hNQedtO0PsCgo9ZzvU5gsSQU6bZM/U0p9/1A09r1doWJAcgX/99RoAkDsCSB7Wrl2rzp07+7azhkr1799fM2fO1J9//qkDBw74jtesWVPffvuthg8frjfeeENVq1bVe++9xxK8AEoUj9vQ6ZMenTrq9oWL0yc8On3craTztl0Zl3eN0AibQkvZFRphV3CYTSHhNgWH2xUSZlPI2X3B4XaFhNsUcs7nweF2BYfa5AgouIDgchlyhrsUEm5XQCDBAwDygwCSh06dOmVfbeI8ud3lvFOnTtqwYUMhVgUA5jEMQ8mJhk4ddXs/jrh18qjn7L/efYnHPPK48/+czlCbIsvaFVbGrvBSdoVG2hVW2q6wSNvZf+0KK/XXR0hEwQYIAID/EUAAAD6ZLkMn491KOOTW8UOZOn7YrYTDbh0/5Nbxw+58zaOw2aXIcnaVKu9QRFm7IsraFVnOroiyjrP//vXhDGG1JgAoaQggAFDCGIahxGMeHdmfqSP7MnVkf6aO/eENGCePuH33TMhLeBm7Sle0q0xFh0pHOVS6osO7ffbzyHJ2eikAAHkigABAMeXxGDrxp1t/7vkraBzd79aR/ZlKT8m7JyPQKZWrEqDy0Q6Vi3aofBWHylUJULloh8pEORToJFwAAC4fAQQAioGMdEPxezJ1eLdLh3Zl6tBulw7vzjto2B1S+WiHomoEqGKNAFWs7lD5aG/IiCxnZyUnAEChIYAAgMW40g0d2uXS/q0uHdjmDRxHD2TmOnQqIEiKiglQVI2zH2c/L1/VwapNAABTEEAAoAjzeAwlHHRr/29ZgSNDh3dnyp2Z89zw0nZVqRug6DoBiq4bqCp1AlSxegDzMQAARQoBBACKEHemoUO7MvX7xgz9vjFDe3/NUMrpnMOowsvYVaNRoKo3DFS1+oGqUjeAoVMAAEsggACAidyZhhIPhuvHj1O179fT2rvZpfTU7IEjMEiqWj/QGzgaef8tU8lB2AAAWBIBBAD8LOFQpravytCO1enatT5D6SmNJaX4jgeH21TrqiDViQ1SrdhAVa0XyDAqAECxQQABgEKWnurRzrUZ2rEqXdtXZ+j44ey3Cg8IcalByzDVvdqp2rFBqlwrQHYHgQMAUDwRQACgEJw+4dZvv6Rr88/p2rU2Xa6Mv47ZHVLNpoGq39qpOlc7tGn3QvXq1VOBgYHmFQwAgJ8QQACggBz7I1O//pSmLcvStX+LS8Y5UznKVnaoUVun6rcOUp2rgxQcapckuVwu/fq7SQUDAGACAggAXIET8W5tWJSqDYvSdGhn9rVxq9YPUNMOwWrSwanKtQKYNA4AgAggAHDJko67tXFxmjYsStO+zS7ffrtDqnN1kDd0tHeqdEWHiVUCAFA0EUAAIB8yXYa2Lk/Xqm9TtW1luu+u4zabVDs2SM27Buuq64IVXtpubqEAABRxBBAAuIDDv7u0en6q1i5IVfKpvyZ11GgUqOZdgxXbJVilytPTAQBAfhFAAOA8GemG1selavncFP2x/a95HZHl7Gp1Y4ha9wxRxer8+gQA4HLwPygAnHX8cKZ+mZuqVd+kKCXJ29thd0hN2jvVuleIGrR2ckNAAACuEAEEQIlmGIZ2rs3Qz5+laOvydN/SuWUrO9SuT4ja9AxVeBnmdQAAUFAIIABKJHemoU2L0/Tj7GQd2vXXMKv6rYPU/rZQNWrr5G7kAAAUAgIIgBIlI83Q6m9TtPh/KTrxp1uSFBRsU+teIWp/W6iiavBrEQCAwsT/tABKhLQUj37+LEVLP0n2rWYVVtqmjneE6dpbQxVWimFWAAD4AwEEQLGWnurRsi9StHh2spITvcGjbGWHOt8dqta9QhUUzDArAAD8iQACoFjKSDe04qsU/TArWWdOeu8aWKGqQ90Hhqv59cGsZgUAgEkIIACKFY/b0Or5qfp++hklJniDR9nKDvUYGKYW3UMIHgAAmIwAAqDY2L4qXfOmnNafe7yrWpWuaFf3AeFq3ZPgAQBAUUEAAWB5h393ad6U09qxOkOSFBphU/eB4bq2T6gCgggeAAAUJQQQAJaVkuTR/HdOa/m8VBkeyREgdbg9VF37hyssklWtAAAoigggACzH4zG05rtUff32ad+Sus06B+umv4WrfDS/1gAAKMr4nxqApRza5dLnk5K0d7NLkhQVE6A7RkSoztVOkysDAAD5QQABYAkZ6Ya+f++0lsxJkeGRgkJsumFguDr2DWWCOQAAFkIAAVDk/b4pQ3MmJOrYQbck73CrPn+PUOmKDpMrAwAAl4oAAqDISk/x6JtpZ7Ts8xRJUqkKdt35VKQatws2uTIAAHC5CCAAiqRd69P1vwlJOvGnt9fjmt4hunlIhELCWd0KAAArI4AAKFIyXYbmv3tGS/6bLMPw3sW879ORqt+KSeYAABQHBBAARcaR/Zn6aOwpHdzpvZN525tDdMvQCDlD6fUAAKC4IIAAMJ1hGFrxVarmvpUkV7oUVsqmu/5ZSk07MtcDAIDihgACwFSpZzz634RE/bo0XZJUr1WQ7v1XKZUqzwpXAAAURwQQAKY5vNulGc+dUsJBtxwB0k1/i1DHvqGy27mvBwAAxRUBBIAp1nyfqk8nJsqVLpWJsmvAS2VUvWGg2WUBAIBCRgAB4FeZGYa+fDNJy+emSpIatA7S/WNKK6wUE80BACgJCCAA/CYxwa33nzmlA9tcstmk7gPD1b1/mOwOhlwBAFBSEEAA+MUfO1yaPvKkEo95FBpp0/2jS6vhNdzbAwCAkoYAAqDQbVqSpo9fPCVXuhQV49DDL5dR+Wh+/QAAUBLxDgBAoTEMQ3EfJOu7985Ikhq0CVK/saUVEs58DwAASioCCIBC4c409L//S9Ta79MkSR37hurmwRFyBDDfAwCAkowAAqDApad4NPO5U9q+OkN2h3T7iEi1uyXU7LIAAEARQAABUKBOn3DrnadO6uCOTAUF29TvhVJq3C7Y7LIAAEARQQABUGCO/ZGpaU+c1PHDboWVtmnQK2VUo1GQ2WUBAIAihAACoEAc3u3S28NP6sxJj8pVcejR18qoQjV+xQAAgOx4dwDgih3Y5tK0ESeUctpQdN0APfpaGUWUdZhdFgAAKIIIIACuyJ5NGXrnqZNKTzFUo3GgHn21jEIiWGYXAADkjgAC4LLtWJOu6SNPypUu1WkepIdeLq3gUMIHAADIGwEEwGXZtiJd7//rpDIzvDcYHDi+jIKc3OMDAABcGAEEwCXbseav8NG0g1P9xpZWQBDhAwAAXBwBBMAl2b3eO+wqM0Nq0sGp/i+W5u7mAAAg3xisDSDf9mzK0Lv/PCVXutSwrVP9xxI+AADApSGAAMiXfVu8q11lpBqq1ypIA19i2BUAALh0BBAAF/XnHpfeedK71G6dq4P00IQyCmTCOQAAuAwEEAAXdCLerWkjTir1jKGYJoF6+P9KKyiY8AEAAC4PAQRAns6c9Gjq8BNKTPCoUs0APfxKGTm5zwcAALgCvJMAkKv0FI/effqkjv3hVpkoux59rYzCIvmVAQAArgzvJgDkkOkyNOO5UzqwzaWwUjY9OqmsSld0mF0WAAAoBgggALIxDENzXk7UjtUZCgqxadDEMoqqwS2DAABAwSCAAMjmhw+Ttfb7NNkd0sCXSqtGoyCzSwIAAMUIAQSAz4ZFqZr/7hlJ0u3DI9WgjdPkigAAQHFDAAEgyXujwdnjEiVJne4KVbs+oSZXBAAAiiMCCAAdP5yp6SNPKTNDatLeqd6DI8wuCQAAFFMEEKCES0vx6N2nT+nMKY+i6wXo/tGlZHdwo0EAAFA4CCBACebxGJr9UqKO7MtUqfJ2PfwyNxoEAACFi3caQAm2aFayNv+ULkegNOCl0ipdgXt9AACAwkUAAUqorSvS9d173hWv7hgRqZgmLLcLAAAKHwEEKIGOHczUrLGnZBhS21tCdE1vVrwCAAD+QQC5gClTpigmJkbBwcFq06aNVq9enee5LpdLL7zwgmrXrq3g4GA1a9ZM33//vR+rBfInPcWj9585pbQzhmKaBuq2YZFmlwQAAEoQAkge5syZoxEjRmjMmDFav369mjVrph49eujo0aO5nv/cc89p2rRpeuutt7R161b97W9/06233qoNGzb4uXIgb4Zh6JNXkhS/N1OR5ewa8GJpBQSy4hUAAPAfAkgeJk2apEGDBmngwIFq1KiRpk6dqtDQUL3//vu5nj9r1iw9++yz6tmzp2rVqqXHHntMPXv21GuvvebnyoG8rfw6Vet/SJPdIfV/sbRKlWfSOQAA8C8CSC4yMjK0bt06de3a1bfPbrera9euWrFiRa6PSU9PV3BwcLZ9ISEhWrZsWaHWCuTX4d0ufTk5SZLU85Fw1bqKSecAAMD/AswuoChKSEiQ2+1WVFRUtv1RUVHavn17ro/p0aOHJk2apI4dO6p27dpatGiRvvjiC7nd7jyvk56ervT0dN92UpL3zaHL5ZLL5SqAV3JhWdfwx7VQOPLbhmkphmaMOiVXhtSgTaDa3xFEuxcR/BxaG+1nfbSh9fm7DfleuXIEkALyxhtvaNCgQWrQoIFsNptq166tgQMH5jlkS5ImTJigsWPH5ti/cOFChYb6b1WiuLg4v10LheNCbWgY0rZvayvhj/IKCs9Q2Zbr9P33mX6sDvnBz6G10X7WRxtan7/aMCUlxS/XKc4IILkoX768HA6Hjhw5km3/kSNHVKlSpVwfU6FCBc2dO1dpaWk6fvy4qlSpopEjR6pWrVp5XueZZ57RiBEjfNtJSUmqVq2aunfvrsjIwl+ZyOVyKS4uTt26dVNgYGChXw8FLz9tuPrbNC3dliy7XXpofHnVbNrdz1XiQvg5tDbaz/poQ+vzdxtmjVjB5SOA5CIoKEgtWrTQokWL1KdPH0mSx+PRokWLNHTo0As+Njg4WNHR0XK5XPr888/Vt2/fPM91Op1yOp059gcGBvr1l6C/r4eCl1cbHtmfqa/eSpbknfdR72ru91FU8XNobbSf9dGG1uevNuT75MoRQPIwYsQI9e/fXy1btlTr1q01efJkJScna+DAgZKkfv36KTo6WhMmTJAkrVq1SocOHVJsbKwOHTqk559/Xh6PR08//bSZLwMlWKbL0EdjT8mVLtVrFaTO94aZXRIAAAABJC933XWXjh07ptGjRys+Pl6xsbH6/vvvfRPTDxw4ILv9r0XE0tLS9Nxzz2nPnj0KDw9Xz549NWvWLJUuXdqkV4CSbsH7Z3RwZ6ZCI22691+lZLdzvw8AAGA+AsgFDB06NM8hV0uWLMm2fd1112nr1q1+qAq4uD2bMrToY+/QqzufKsX9PgAAQJHBfUCAYiYt2aOPX0qU4ZFa9wxRbOfgiz8IAADATwggQDHzxeQknfjTrbKVHbr18QizywEAAMiGAAIUI5uWpGnNd2my2aX7niul4DB+xAEAQNHCuxOgmDhzyqPPXvWuTX79fWGq1SzI5IoAAAByIoAAxcSXk5N05pRHlWoGqMfAcLPLAQAAyBUBBCgGfluWofU/eIde3fNsKQUEseQuAAAomliGF7A4V6pDX0w/I0nqfE+YqjfkDq0AAKDoogcEsLjdi2vo9AlDFWs4dMODDL0CAABFGwEEsLBtKzJ05LcKstmke54ppUAnQ68AAEDRRgABLCot2aMvJnnvdt7+jmDFNGHVKwAAUPQRQACLmv/uGSUmeBRcOk09Hgw1uxwAAIB8IYAAFvTHdpeWfZEiSarXba+Cghl6BQAArIFVsACL8bgNfTIxUYZHiu0SpNIxSWaXBAAAkG/0gAAWs+zLFB3ckangcJt6DwkzuxwAAIBLQgABLOTUMbfmv+O958dNj0Yooiw/wgAAwFp49wJYyNw3kpSeYqhGo0C1vSXE7HIAAAAuGQEEsIitK9K1aUm67A7pzqciZbcz8RwAAFgPAQSwgMwMQ1++4Z1s3vGOUEXXDTS5IgAAgMtDAAEsYOknyUo46FZEObt6PBhudjkAAACXjQACFHGJCW7FfeC943nvv0UoOIwfWwAAYF28kwGKuG/ePq30VEM1GgeqRY9gs8sBAAC4IgQQoAjbtyVDaxekSZJuG8bEcwAAYH0EEKCI8ngMff66d+J5m14hqt6QiecAAMD6CCBAEbV6fqr3judhNvV6lInnAACgeCCAAEVQWrJH30713vG8x4PhiijrMLkiAACAgkEAAYqgRR8l68wpjypWd6jD7aFmlwMAAFBgCCBAEXPyiFtL53iX3b3psQg5Aph4DgAAig8CCFDEfPfeabkypFrNAtWkvdPscgAAAAoUAQQoQg7tcmnt995ld28eEiGbjd4PAABQvBBAgCLCMAzNm3JahiE1vz5YNRoFmV0SAABAgSOAAEXE9lUZ2rk2Q45AsewuAAAotgggQBHgcRv6+j+nJUkdbgtVuSoBJlcEAABQOAggQBGw5rtU/bknU6ERNnXtT+8HAAAovggggMlc6Ya+f99708Gu/cIVFsmPJQAAKL54pwOYbPlXKTp11KNSFexqfxs3HQQAAMUbAQQwUXqqRz/M8t50sPuAcAU6WXYXAAAUbwQQwEQ/f5aiMyc9Kh/tUJteIWaXAwAAUOgIIIBJUk979OPH3t6PGx4KlyOA3g8AAFD8EUAAkyz+b7JSzxiqXCtAza8PNrscAAAAvyCAACY4fdKtpZ+mSJJufDhcdge9HwAAoGQggAAm+GFWsjJSDVVvGKgmHZxmlwMAAOA3BBDAz04dc2v5XG/vR89HwmWz0fsBAABKDgII4Gc/fpyszAypVrNA1WsZZHY5AAAAfkUAAfwoMcGtlfO8vR89BtL7AQAASh4CCOBHi/+bLFeGFNM0UHVb0PsBAABKHgII4CenT/w194PeDwAAUFIRQAA/WfzfZLnSpRqNAlW/Fb0fAACgZCKAAH5w5qRHv3yZKknqTu8HAAAowQgggB8smZOsjDRD1RoEqOE19H4AAICSiwACFLLkRI+Wfe6d+9F9AL0fAACgZCOAAIVs6SfJSk81FF03QI2v5a7nAACgZCOAAIUoLeWv3o9u/en9AAAAIIAAhWjFV6lKPWOoYnWHmnak9wMAAIAAAhSSzAxDS+ckS5K63Bsmu53eDwAAAAIIUEjWLkxVYoJHpSrY1aJ7iNnlAAAAFAkEEKAQeNyGfvzY2/vR6a4wBQTR+wEAACARQIBCsfnndB37w63QCJuuuZneDwAAgCwEEKCAGYahRR+dkSS1vz1UwaH8mAEAAGQpdu+MUlNTlZqaanYZKMF2rcvQH9szFeiUOtweZnY5AAAARUqxCCAHDhzQwIEDFRUVpfDwcIWHhysqKkoPPvig9u/fb3Z5KGEWnZ370eamUIWXKRY/YgAAAAUmwOwCrtT27dvVvn17nTp1St26dVPDhg19+z/88EN9/fXXWrZsmerXr29ypSgJDu50aeeaDNkdUue7Q80uBwAAoMixfAAZOXKk7Ha7NmzYoKZNm2Y7tmXLFl1//fUaOXKkvvzyS5MqREmSdd+P2M7BKlvZ8j9eAAAABc7y40OWLl2qf/zjHznChyQ1adJEQ4cO1ZIlS/xfGEqcxAS3NixKkyRddzdzPwAAAHJj+QDicrkUEpL3MqehoaFyuVx+rAgl1bLPU+TOlGpeFajqDQLNLgcAAKBIsnwAad68ud577z0lJibmOJaUlKTp06fr6quvNqEylCTpqR4t/ypFkvfGgwAAAMid5Qepjx07VjfccIMaNGiggQMHql69epKkHTt26IMPPtDx48c1ZcoUk6tEcbf2+zSlJBkqV8WhJu2dZpcDAABQZFk+gHTp0kXz58/XU089pf/7v//Ldiw2NlazZs1S586dTaoOJYHHY2jpJ97J5x37hsrusJlcEQAAQNFl+QAiSV27dtWGDRsUHx/vu+9HjRo1VKlSJZMrQ0mwbUW6jv3hVnC4TW165j0fCQAAAMUkgGSpVKkSoQN+t2SOd+5H25tD5Qy1/LQqAACAQmW5APLhhx9Kkh544AHZbDbf9sX069evMMtCCXVol0u713tvPNjhdm48CAAAcDGWCyADBgyQzWbT3XffraCgIA0YMOCij7HZbAQQFIqsuR+xnYNVJsphcjUAAABFn+UCyN69eyVJQUFB2bYBfztz0uO78WDHO+n9AAAAyA/LBZAaNWpccBvwl5XfpCgzQ6reMFA1GgeZXQ4AAIAlWH7GbK1atTRv3rw8j3/zzTeqVavWZT33lClTFBMTo+DgYLVp00arV6++4PmTJ09W/fr1FRISomrVqmn48OFKS0u7rGujaHNnGvrlS+/k8/a30fsBAACQX5YPIPv27dOZM2fyPH7mzBnf0ryXYs6cORoxYoTGjBmj9evXq1mzZurRo4eOHj2a6/mzZ8/WyJEjNWbMGG3btk3Tp0/XnDlz9Oyzz17ytVH0/fZLuk4d9Si8tF2xXYLNLgcAAMAyLB9AJO8k87ysWbNGpUuXvuTnnDRpkgYNGqSBAweqUaNGmjp1qkJDQ/X+++/nev7y5ct17bXX6t5771VMTIy6d++ue+6556K9JrCmnz/39n5c0ztEgU5uPAgAAJBflgwgb7zxhmrVqqVatWrJZrNp2LBhvu1zP8qVK6fJkyerZ8+el/T8GRkZWrdunbp27erbZ7fb1bVrV61YsSLXx7Rr107r1q3zBY49e/Zo/vz5l3xtFH1/7vlr6d12fRh+BQAAcCksNwldkipWrKjGjRtL8g7Bio6OVnR0dLZzbDabwsLC1KJFCw0ePPiSnj8hIUFut1tRUVHZ9kdFRWn79u25Pubee+9VQkKC2rdvL8MwlJmZqb/97W8XHIKVnp6u9PR033ZSUpIkyeVyyeVyXVLNlyPrGv64VnHy02feIX+Nrw1SeFmPXC6PabXQhtZHG1ob7Wd9tKH1+bsN+V65cjbDMAyzi7gSnTt31nPPPafrr7++wJ7z8OHDio6O1vLly9W2bVvf/qefflpLly7VqlWrcjxmyZIluvvuu/XSSy+pTZs22r17tx5//HENGjRIo0aNyvU6zz//vMaOHZtj/+zZsxUayl/WiyJXmkMrpjaXx+VQs7u2qkz102aXBAAA/CglJUX33nuvEhMTFRkZaXY5lmT5AFIYMjIyFBoaqs8++0x9+vTx7e/fv79OnTqlr776KsdjOnTooGuuuUYTJ0707fvoo4/0yCOP6MyZM7Lbc452y60HpFq1akpISPDLN7TL5VJcXJy6deumwMDAQr9ecfDzZ6n6ekqKKtV0aPj0Uhecf+QPtKH10YbWRvtZH21off5uw6SkJJUvX54AcgUsOQQrNy6XS9u3b1diYqI8npxDYjp27Jjv5woKClKLFi20aNEiXwDxeDxatGiRhg4dmutjUlJScoQMh8N7Z+y8Mp7T6ZTT6cyxPzAw0K+/BP19PavyeAyt+MobGDvcHua7GWZRQBtaH21obbSf9dGG1uevNuT75MpZPoB4PB4988wz+s9//qOUlJQ8z3O73Zf0vCNGjFD//v3VsmVLtW7dWpMnT1ZycrIGDhwoSerXr5+io6M1YcIESVLv3r01adIkNW/e3DcEa9SoUerdu7cviMDadq7NUMJBt4LDbWrRnaV3AQAALoflA8j48eM1ceJEPfroo2rfvr0eeOABvfzyyypdurT+85//yGaz6ZVXXrnk573rrrt07NgxjR49WvHx8YqNjdX333/vm5h+4MCBbD0ezz33nGw2m5577jkdOnRIFSpUUO/evTVu3LgCe60w14qvvAG3ZY8QOUMtuYAcAACA6SwfQGbOnKm+ffvq7bff1vHjxyVJLVq0UJcuXdS/f3+1bdtWP/74Y7YldfNr6NCheQ65WrJkSbbtgIAAjRkzRmPGjLnk66DoS0xwa8sy7/CrdreEmFwNAACAdVn+z7gHDx5Uly5dJMk3nyItLU2Sdy7H/fffr1mzZplWH4qH1d+myuOWYpoGqnItxn4CAABcLssHkHLlyunMGe99GcLDwxUZGak9e/ZkO+fkyZNmlIZiwuM2tOJr7/CrdjezPDIAAMCVsPwQrObNm2vNmjW+7c6dO2vy5Mlq3ry5PB6P3nzzTTVr1szECmF1O1Zn6GS8R6ERNjXrwuRzAACAK2H5HpBHHnkk2/00xo0bp1OnTqljx4667rrrlJSUpNdee83kKmFly+ednXx+Q4iCnObe9wMAAMDqLN8DcvPNN+vmm2/2bTdq1Ei///67lixZIofDoXbt2qls2bImVggrO3XMra3LsyafM/wKAADgSlk+gOSmVKlSuuWWW3zbP/300yXdiBDIsuob7+TzWs0CFRVTLH9cAAAA/MryQ7AuZN68ebr22mvVuXNns0uBBXnchlZmTT6n9wMAAKBAWDaAxMXF6aabblLDhg3Vrl07vf76675jc+fOVZMmTXTrrbdq165d3JsDl2XbynSdOupRWCmbrrqOyecAAAAFwZJjSubPn6/evXvLMAyVL19eu3fv1qpVq3T06FGlpKTorbfeUu3atTVlyhQNGDBAwcG8ecSlW/5VqiSp1Q0hCmTyOQAAQIGwZAB55ZVXVKVKFcXFxalBgwZKTEzU3Xffrddff102m03//ve/9eijj8rhcJhdKizq5BG3tq30Tj5vy/ArAACAAmPJIVgbNmzQY489pgYNGkjyTjp/6aWXlJGRoWeffVaDBw8mfOCKrPkuVYZHqh0bqIrVLZnTAQAAiiRLBpDTp0+rRo0a2fZlbbdq1cqMklCMeDyGVs/3Dr9qcxO9HwAAAAXJkgFEkmw2W67bQUFBZpSDYmTPxgwdP+yWM9SmZp2YPwQAAFCQLDu25MMPP9TKlSt922lpab75H3Pnzs12rs1m0xtvvOHnCmFVq7719n5c3TVYQcFMPgcAAChIlg0gCxcu1MKFC3PsPz98SAQQ5F/qGY82LUmTJLXpFWJyNQAAAMWPJQOIx+MxuwQUUxsWpcmVLkXFBKh6o0CzywEAACh2LDsHBCgMq88Ov2rTMyTHPCMAAABcOQIIcFb83kzt3+qS3SG1vIHJ5wAAAIWBAAKctXp+iiSpUTunIspyHxkAAIDCQAABJLkzDa35nsnnAAAAhY0AAkjauiJdZ056FFHWrobXOM0uBwAAoNgigAD6a/J5yxtC5Ahg8jkAAEBhKVYB5M8//9SmTZuUnJxsdimwkKTjbm1dkS7Ju/oVAAAACk+xCCBfffWVGjRooKpVq+rqq6/WqlWrJEkJCQlq3rx5rjcnBLKsW5gmj1uq0ThQUTGWvDUOAACAZVg+gHz99de67bbbVL58eY0ZM0aGYfiOlS9fXtHR0ZoxY4aJFaKoW7vAO/yq1Q30fgAAABQ2yweQF154QR07dtSyZcs0ZMiQHMfbtm2rDRs2mFAZrODwbpcO786UI0CKvZ57fwAAABQ2yweQLVu2qG/fvnkej4qK0tGjR/1YEawkq/ejUTunwiIt/+MAAABQ5Fn+HVdoaOgFJ53v2bNH5cqV82NFsAqP29C6hd57fzD8CgAAwD8sH0A6d+6sDz74QJmZmTmOxcfH691331X37t1NqAxF3c51GUo67lFopE0N23LvDwAAAH+wfAAZN26cDh48qFatWmnatGmy2WxasGCBnnvuOTVt2lSGYWjMmDFml4kiKGv4VfMuwQoI5N4fAAAA/mD5AFK/fn0tW7ZM5cqV06hRo2QYhiZOnKjx48eradOm+vnnnxUTE2N2mShi0lM82rzUe++Plgy/AgAA8JticdODxo0b64cfftDJkye1e/dueTwe1apVSxUqVDC7NBRRvy5NV0aaoQpVHarRONDscgAAAEqMYhFAspQpU0atWrUyuwxYQNbwqxY9QmSzMfwKAADAXyw/BGvRokWaOHFitn3vv/++qlevrqioKA0fPlxut9uk6lAUnTrm1q51GZKkFt259wcAAIA/WT6APP/889q0aZNve/PmzXr00UdVoUIFderUSW+++aZeffVVEytEUbNuYaoMQ6p5VaDKRxerTkAAAIAiz/IBZNu2bWrZsqVve9asWYqMjNTPP/+sOXPmaNCgQfrwww9NrBBFiWEYWreAe38AAACYxfIBJDk5WZGRkb7t77//XjfccINCQ0MlSa1atdL+/fvNKg9FzOHdmfpzT6YCgqTYzgy/AgAA8DfLB5Bq1appzZo1kqTdu3dry5Yt2W48eOLECTmd3GQOXusWeiefN2rnVEiE5b/9AQAALMfyA+Dvu+8+vfDCCzp06JB+++03lSlTRrfccovv+Lp161SvXj0TK0RR4fEY2rDIO/yqRTeGXwEAAJjB8gHkX//6lzIyMjR//nxVr15dM2fOVOnSpSV5ez+WLFmixx9/3NwiUSTs3ezSqaMeBYfZ1PAaesUAAADMYPkAEhAQoHHjxmncuHE5jpUtW1bx8fEmVIWiaH2cd/jVVdcFK9DJvT8AAADMwCB4lAjuTEObFnuHX13djcnnAAAAZrF8D4gkpaWl6fPPP9f69euVmJgoj8eT7bjNZtP06dNNqg5Fwc41GUpONBRexq46zYPMLgcAAKDEsnwA2b9/vzp37qx9+/apdOnSSkxMVNmyZXXq1Cm53W6VL19e4eHhZpcJk63/wTv8KrZLsBwBDL8CAAAwi+WHYD311FNKTEzUypUrtXPnThmGoTlz5ujMmTN6+eWXFRISogULFphdJkyUkWZo80/pkqQWDL8CAAAwleUDyI8//qjBgwerdevWstu9L8cwDDmdTj311FO6/vrrNWzYMHOLhKm2Lk9XeqqhspUdqtE40OxyAAAASjTLB5CUlBTFxMRIkiIjI2Wz2ZSYmOg73rZtWy1btsyk6lAUZA2/an59sGw2hl8BAACYyfIBpHr16jp48KAk75K80dHRWrlype/41q1bFRzMsJuSKvW0R1tXeIdfsfoVAACA+Sw/Cb1Lly766quvNGbMGEnSgAEDNGHCBJ08eVIej0ezZs1Sv379TK4SZvn1pzS5XVKlmgGqUpvhVwAAAGazfAAZOXKk1qxZo/T0dDmdTj377LM6fPiwPvvsMzkcDt17772aNGmS2WXCJOt/OHvvj670fgAAABQFlg8g1atXV/Xq1X3bwcHBeu+99/Tee++ZWBWKgtMn3Nq1LkOS1JwAAgAAUCRYfg7I+RITE+V2u80uA0XAxsVpMjxSjUaBKh9t+awNAABQLBSLALJ27VrdcMMNCg0NVbly5bR06VJJUkJCgm655RYtWbLE3AJhio0/eodfxV5P7wcAAEBRYfkAsnz5crVv3167du3S/fffL4/H4ztWvnx5JSYmatq0aSZWCDMkJri191eXJCm2MwEEAACgqLB8AHn22WfVsGFDbd26VePHj89xvHPnzlq1apUJlcFMm5akyTCkmCaBKl3RYXY5AAAAOMvyAWTNmjUaOHCgnE5nrjeZi46OVnx8vAmVwUybFp8dftWF3g8AAICixPIBJDAwMNuwq/MdOnRI4eHhfqwIZjt3+FWzTgQQAACAosTyAeSaa67RZ599luux5ORkzZgxQ9ddd52fq4KZfMOvmjL8CgAAoKixfAAZO3as1q5dq169eum7776TJG3atEnvvfeeWrRooWPHjmnUqFEmVwl/8g2/YvI5AABAkWP5myO0adNG8+fP12OPPaZ+/fpJkp544glJUu3atTV//nxdddVVZpYIP2L4FQAAQNFm6QBiGIZOnz6tdu3aaceOHdq4caN27dolj8ej2rVrq0WLFrlOTEfxxfArAACAos3SASQjI0Nly5bV+PHj9fTTTys2NlaxsbFmlwUT+W4+yPArAACAIsnSc0CcTqcqVaokp9NpdikoAk4dc2vfZoZfAQAAFGWWDiCSNGDAAH344YfKyMgwuxSY7FeGXwEAABR5lh6CJUlNmzbV3Llz1bhxYw0YMEAxMTEKCQnJcd5tt91mQnXwp42sfgUAAFDkWT6A3HPPPb7P81pu12azye12+6skmODUMVa/AgAAsALLB5DFixebXQKKgF+XeHs/ajL8CgAAoEizfADhLueQ/lr9qhnDrwAAAIo0y09CP3HihH799dc8j2/evFknT570Y0Xwt8QEt/ay+hUAAIAlWD6ADB8+XI888kiexx999FE9+eSTfqwI/rb5p3RJUkwThl8BAAAUdZYPID/++KNuvvnmPI/37t1bP/zwgx8rgr/9utQ7/OqqjvR+AAAAFHWWDyDHjh1T+fLl8zxerlw5HT161I8VwZ+SEz36faP3HjBNr+OGlAAAAEWd5QNI5cqVtWHDhjyPr1u3ThUqVPBjRfCn335Jk8ctVakToPLRll9TAQAAoNizfADp06ePpk+frnnz5uU49tVXX2nGjBm69dZbL+u5p0yZopiYGAUHB6tNmzZavXp1nud26tRJNpstx0evXr0u69rIn1+Xeud/XHUdw68AAACswPJ/Mn7++ef1ww8/6NZbb1WzZs3UpEkTSdKWLVu0adMmNWzYUGPHjr3k550zZ45GjBihqVOnqk2bNpo8ebJ69OihHTt2qGLFijnO/+KLL5SRkeHbPn78uJo1a6Y777zz8l8cLigtxaMda7ICCMOvAAAArMDyPSClSpXSypUr9dxzz8nlcumzzz7TZ599JpfLpVGjRmnVqlUqXbr0JT/vpEmTNGjQIA0cOFCNGjXS1KlTFRoaqvfffz/X88uWLatKlSr5PuLi4hQaGkoAKUTbV2YoM0OqUNWhSjUtn6UBAABKhGLxri0sLExjx469rJ6O3GRkZGjdunV65plnfPvsdru6du2qFStW5Os5pk+frrvvvlthYWEFUhNyylr9qmnHYNlsNpOrAQAAQH4UiwCS5c8//9TRo0dVp06dK3rjn5CQILfbraioqGz7o6KitH379os+fvXq1dqyZYumT59+wfPS09OVnp7u205KSpIkuVwuuVyuy6j80mRdwx/XKmiuDEO/LfcGkEbXOiz5GgqCldsQXrShtdF+1kcbWp+/25DvlStXLALIV199pX/+85/atWuXJCkuLk5dunRRQkKCunXrpjFjxqhPnz5+q2f69Olq2rSpWrdufcHzJkyYkGuvzcKFCxUaGlpY5eUQFxfnt2sVlOO/l1ZGan0FhWdo8944bdlndkXmsmIbIjva0NpoP+ujDa3PX22YkpLil+sUZ5YPIF9//bVuu+02tW3bVvfee6+ef/5537Hy5csrOjpaM2bMuKQAUr58eTkcDh05ciTb/iNHjqhSpUoXfGxycrL+97//6YUXXrjodZ555hmNGDHCt52UlKRq1aqpe/fuioyMzHe9l8vlcikuLk7dunVTYGBgoV+vIH36yhlJ6WrZLVK9evU0uxzTWLkN4UUbWhvtZ320ofX5uw2zRqzg8lk+gLzwwgvq2LGjFi9erOPHj2cLIJLUtm1bTZs27ZKeMygoSC1atNCiRYt8wcXj8WjRokUaOnToBR/76aefKj09Xffff/9Fr+N0OuV05ly9KTAw0K+/BP19vSvlzjS0dbl3xbHYziGWqr2wWK0NkRNtaG20n/XRhtbnrzbk++TKWX4VrC1btqhv3755Ho+KirqsO6GPGDFC7777rj744ANt27ZNjz32mJKTkzVw4EBJUr9+/bJNUs8yffp09enTR+XKlbvkayJ/9vyaoeREQ2GlbKp1VZDZ5QAAAOASWL4HJDQ0VMnJyXke37Nnz2WFgbvuukvHjh3T6NGjFR8fr9jYWH3//fe+iekHDhyQ3Z49v+3YsUPLli3TwoULL/l6yL+smw82aR8sRwCrXwEAAFiJ5QNI586d9cEHH2jYsGE5jsXHx+vdd9/VTTfddFnPPXTo0DyHXC1ZsiTHvvr168swjMu6FvLH4zG0+aezy+9y80EAAADLsfwQrHHjxungwYNq1aqVpk2bJpvNpgULFui5555T06ZNZRiGxowZY3aZKCB/bHMp8ZhHzlCb6rUggAAAAFiN5QNI/fr1tWzZMpUrV06jRo2SYRiaOHGixo8fr6ZNm+rnn39WTEyM2WWigPz6k3f4VaO2TgU6GX4FAABgNZYfgiVJjRs31g8//KCTJ09q9+7d8ng8qlWrlipUqGB2aShgvuFXHen9AAAAsKJiEUCylClTRq1atTK7DBSSI/szdewPtxwBUsNrCCAAAABWZOkAkp6ero8++kgLFy7U77//rtOnTysiIkJ16tTRDTfcoHvvvVdBQSzTWlxsWebt/ah7dZCCwyw/ehAAAKBEsmwA2bx5s2655Rbt379fhmGoVKlSCg8P19GjR7V+/Xp9+umnGjdunObNm6eGDRuaXS4KwJZl3vkfjdsHm1wJAAAALpcl/4x85swZ3XzzzTpy5IjGjRunP/74QydPnsz270svvaTDhw+rd+/eF7xPCKzh9Am39m9xSZKatGf4FQAAgFVZMoDMmDFDBw4c0LfffquRI0cqOjo62/Ho6Gg988wz+vrrr7V3717NnDnTnEJRYH5bni7DkKrWD1Dpig6zywEAAMBlsmQA+fbbb9W9e3d16tTpgud16dJF3bp109dff+2fwlBofjs7/KppB4ZfAQAAWJklA8jmzZsvGj6ydOnSRZs3by7cglCoMtIM7Vhzdv7HtQy/AgAAsDJLBpATJ06oUqVK+To3KipKJ06cKOSKUJh2rEmXK10qU8muKnUsu24CAAAAZNEAkp6ersDAwHydGxAQoIyMjEKuCIVpy8/e3o8m7YNls3H3cwAAACuz7J+T9+3bp/Xr11/0vL179/qhGhQWj9vQ1uVZ8z8YfgUAAGB1lg0go0aN0qhRoy56nmEY/NXcwvb95tKZUx6FhNtUqxk3lQQAALA6SwaQGTNmmF0C/CTr5oMN2zrlCCBIAgAAWJ0lA0j//v3NLgF+suXnNEncfBAAAKC4sOQkdJQMRw9k6tgfbjkCpIbXEEAAAACKAwIIiqzNZ3s/6lwdpOAwvlUBAACKA97VocjKmv/RpD13PwcAACguCCAokk6fcGv/Fpck5n8AAAAUJwQQFElbl6fLMKSq9QJUuqLD7HIAAABQQAggKJJ8w686MPwKAACgOCGAoMjJSDe0Y03W/A+GXwEAABQnBBAUObvXZ8iVLpWuaFeVOpa8VQ0AAADyQABBkbP1F+/yu43aOWWzcfdzAACA4oQAgiLFMAz9ttw7/KpRO4ZfAQAAFDcEEBQpf/6eqVNHPQp0SnVbEEAAAACKGwIIipStK7y9H3VbOBXkZPgVAABAcUMAQZGSNfyqMcOvAAAAiiUCCIqMM6c8vrufN2xLAAEAACiOCCAoMrat9N79vEqdAJWJ4u7nAAAAxREBBEXGVoZfAQAAFHsEEBQJ7kxD21ez/C4AAEBxRwBBkbD31wylnTEUXtqu6g0DzS4HAAAAhYQAgiIha/WrBtcEye5g+V0AAIDiigCCIuGv+R/BJlcCAACAwkQAgemOHczU0QNu2R1S/dZBZpcDAACAQkQAgemyej9qxwYpJJxvSQAAgOKMd3swXVYAYfUrAACA4o8AAlOlJXv0+8YMSVIj7n4OAABQ7BFAYKodazLkzpQqVHWoYvUAs8sBAABAISOAwFS+4VfX0vsBAABQEhBAYBqPx9DWFVnL7xJAAAAASgICCEzzx3aXzpz0KDjMpppXsfwuAABASUAAgWmyhl/VbxWkgEDufg4AAFASEEBgmqzhVyy/CwAAUHIQQGCKpONuHdyRKUlqeA0BBAAAoKQggMAU21d77/1RtX6AIso6TK4GAAAA/kIAgSm2r/QOv6L3AwAAoGQhgMDv3JmGdqwmgAAAAJREBBD43YFtLqWcNhQaYVONRoFmlwMAAAA/IoDA77adXf2qfmun7A6W3wUAAChJCCDwu23M/wAAACixCCDwq6Tjbh3c6V1+t0Eb7n4OAABQ0hBA4FcsvwsAAFCyEUDgVyy/CwAAULIRQOA37kxD21cRQAAAAEoyAgj8Zv9Wl1LPGAqNZPldAACAkooAAr/JGn5VvxXL7wIAAJRUBBD4DcvvAgAAgAACv2D5XQAAAEgEEPhJ1uRzlt8FAAAo2Qgg8IttK733/2jUluFXAAAAJRkBBIXOnWlox2rmfwAAAIAAAj84d/nd6g1ZfhcAAKAkI4Cg0LH8LgAAALIQQFDotrL8LgAAAM4igKBQJR1361DW8rvXsPwuAABASUcAQaHKWn63WoMARZRh+V0AAICSjgCCQpW1/C7DrwAAACARQFCIWH4XAAAA5yOAoNAc2HZ2+d0Ilt8FAACAFwEEhWb72d6Peiy/CwAAgLMIIBcwZcoUxcTEKDg4WG3atNHq1asveP6pU6c0ZMgQVa5cWU6nU/Xq1dP8+fP9VG3Rs2OVd/5Hg9asfgUAAACvALMLKKrmzJmjESNGaOrUqWrTpo0mT56sHj16aMeOHapYsWKO8zMyMtStWzdVrFhRn332maKjo7V//36VLl3a/8UXAclJHh3Y7pIk1W/N/A8AAAB4EUDyMGnSJA0aNEgDBw6UJE2dOlXffvut3n//fY0cOTLH+e+//75OnDih5cuXKzDQO98hJibGnyUXKTvXpMvwSJVqBqh0RZbfBQAAgBdDsHKRkZGhdevWqWvXrr59drtdXbt21YoVK3J9zLx589S2bVsNGTJEUVFRatKkicaPHy+32+2vsouUHasZfgUAAICc6AHJRUJCgtxut6KiorLtj4qK0vbt23N9zJ49e/Tjjz/qvvvu0/z587V7924NHjxYLpdLY8aMyfUx6enpSk9P920nJSVJklwul1wuVwG9mrxlXaOgr2UYhu8GhHVaOPzyWkqqwmpD+A9taG20n/XRhtbn7zbke+XKEUAKiMfjUcWKFfXOO+/I4XCoRYsWOnTokCZOnJhnAJkwYYLGjh2bY//ChQsVGhpa2CX7xMXFFejzJSeEKDHhKtkDPNpxeJF2zzcK9PmRU0G3IfyPNrQ22s/6aEPr81cbpqSk+OU6xRkBJBfly5eXw+HQkSNHsu0/cuSIKlWqlOtjKleurMDAQDkcf813aNiwoeLj45WRkaGgoJxDkZ555hmNGDHCt52UlKRq1aqpe/fuioyMLKBXkzeXy6W4uDh169bNN2+lIPz0SaqkFNVp7lTvW24ssOdFToXVhvAf2tDaaD/row2tz99tmDViBZePAJKLoKAgtWjRQosWLVKfPn0keXs4Fi1apKFDh+b6mGuvvVazZ8+Wx+OR3e6dWrNz505Vrlw51/AhSU6nU05nzhWiAgMD/fpLsKCvt2vdaUlSw2uC+WXuJ/7+nkHBow2tjfazPtrQ+vzVhnyfXDkmoedhxIgRevfdd/XBBx9o27Zteuyxx5ScnOxbFatfv3565plnfOc/9thjOnHihB5//HHt3LlT3377rcaPH68hQ4aY9RJMkZFuaM/GrAnoLL8LAACA7OgBycNdd92lY8eOafTo0YqPj1dsbKy+//5738T0AwcO+Ho6JKlatWpasGCBhg8frquuukrR0dF6/PHH9c9//tOsl2CK3zdkyJUhla5oV1QMy+8CAAAgOwLIBQwdOjTPIVdLlizJsa9t27ZauXJlIVdVtG1f7V39qn5rp2w2m8nVAAAAoKhhCBYK1I6zAYThVwAAAMgNAQQF5uQRt47sc8tml+q15AaEAAAAyIkAggKT1ftRo1GgQiP51gIAAEBOvEtEgcm6+3l9hl8BAAAgDwQQFAh3pqGda7OW32X4FQAAAHJHAEGB+GO7S6lnDIVG2FS9ITfoAQAAQO4IICgQWcOv6rZ0yu5g+V0AAADkjgCCArF99dnhV20YfgUAAIC8EUBwxVKSPDqwzSWJ+38AAADgwggguGI712bI8EhRMQEqXdFhdjkAAAAowggguGLbs+5+zvArAAAAXAQBBFfEMAzfDQgZfgUAAICLIYDgihzZ59apox4FBkm1YukBAQAAwIURQHBFsoZf1YoNUpCT5XcBAABwYQQQXJEdqxh+BQAAgPwjgOCyZaQb+n2j9/4f9Vsz/AoAAAAXRwDBZduzKUOuDKlUBbsq1QwwuxwAAABYAAEEl+3c4Vc2G/M/AAAAcHEEEFy2HWsZfgUAAIBLQwDBZUk67tafv2fKZpPqtmACOgAAAPKHAILLsvNs70d03QCFl+bbCAAAAPnDO0dclqy7n9dn+V0AAABcAgIILplhGL4ekPqtmP8BAACA/COA4JLF781U0nGPAp1SzaYEEAAAAOQfAQSXbMcab+9H7WZBCghi+V0AAADkHwEEl2znGu/8j3qtmP8BAACAS0MAwSXJzDD0+0aXJO7/AQAAgEtHAMEl2bslQxlphiLK2VW5VoDZ5QAAAMBiCCC4JDvPzv+o1zJINhvzPwAAAHBpCCC4JDvOzv+o35L5HwAAALh0BBDkW3KiRwd3ZEqS6nH/DwAAAFwGAgjybde6DBmGVLlWgEqVd5hdDgAAACyIAIJ82+FbfpfeDwAAAFweAgjyxTAM5n8AAADgihFAkC8JB906Ge+RI1CqFRtodjkAAACwKAII8mXH2eV3azYJkjOEbxsAAABcHt5JIl92Mv8DAAAABYAAgotyZxratd7bA1K/NfM/AAAAcPkIILioA9tcSks2FBppU9W6AWaXAwAAAAsjgOCisuZ/1G3hlN1hM7kaAAAAWBkBBBeVNf+jPvM/AAAAcIUIILigtGSP9m91SSKAAAAA4MoRQHBBuzdkyOOWKlR1qGxl5n8AAADgyhBAcEFZ8z/qtWL1KwAAAFw5AgguaMdq5n8AAACg4BBAkKeT8W4d+8Mtu0OqczUBBAAAAFeOAII87Vjr7f2o3jBQIeF8qwAAAODK8a4Sedp5dv4Hw68AAABQUAggyJXHY2jn2R4QJqADAACgoBBAkKtDOzOVnGjIGWpTjUaBZpcDAACAYoIAglztOHv387pXB8kRYDO5GgAAABQXBBDkaufarPt/MP8DAAAABYcAghwy0gzt+TVrAjrzPwAAAFBwCCDIYc+mDLldUpkouypUc5hdDgAAAIoRAghyyJr/Ua+VUzYb8z8AAABQcAggyGEH9/8AAABAISGAIJuk4279+XumbDapbgvmfwAAAKBgEUCQTdbqV9F1AxRemm8PAAAAFCzeYSKbrPkf9VvT+wEAAICCRwCBj2EY2sn8DwAAABQiAgh84vdmKum4R4FOqWZTAggAAAAKHgEEPlmrX9VuFqSAIJbfBQAAQMEjgMBn5zn3/wAAAAAKAwEEkqTMDEO/b3RJYv4HAAAACg8BBJKkfVsylJFmKKKsXZVrB5hdDgAAAIopAggk/TX/o16rINlszP8AAABA4SCAQNI59/9oyfwPAAAAFB4CCJSc6NHBHZmSvD0gAAAAQGEhgEC/b3DJMKRKNQNUqrzD7HIAAABQjBFAoJ1rz65+1ZreDwAAABQuAkgJZxjSrnVnAwjzPwAAAFDICCAlXOopp07Ge+QIlGrFBppdDgAAAIo5AkgJd3JfKUlSzSZBcobw7QAAAIDCxTvOC5gyZYpiYmIUHBysNm3aaPXq1XmeO3PmTNlstmwfwcHBfqz28mQFEFa/AgAAgD8QQPIwZ84cjRgxQmPGjNH69evVrFkz9ejRQ0ePHs3zMZGRkfrzzz99H/v37/djxZfO7TZ08kCkJKl+a+Z/AAAAoPARQPIwadIkDRo0SAMHDlSjRo00depUhYaG6v3338/zMTabTZUqVfJ9REVF+bHiS3dwe6bcGQEKjbSpat0As8sBAABACUAAyUVGRobWrVunrl27+vbZ7XZ17dpVK1asyPNxZ86cUY0aNVStWjXdcsst+u233/xR7mXLWn63TvNA2R02k6sBAABAScCfvXORkJAgt9udowcjKipK27dvz/Ux9evX1/vvv6+rrrpKiYmJevXVV9WuXTv99ttvqlq1aq6PSU9PV3p6um87KSlJkuRyueRyuQro1eRt55oMSVKt5na/XA8FL6vdaD/rog2tjfazPtrQ+vzdhnyvXDkCSAFp27at2rZt69tu166dGjZsqGnTpunFF1/M9TETJkzQ2LFjc+xfuHChQkNDC61WSfK4bYr/o5kkp/5MWaH58zMK9XooXHFxcWaXgCtEG1ob7Wd9tKH1+asNU1JS/HKd4owAkovy5cvL4XDoyJEj2fYfOXJElSpVytdzBAYGqnnz5tq9e3ee5zzzzDMaMWKEbzspKUnVqlVT9+7dFRkZeXnFX4LuPTI075NluvmO6xQYyD1ArMjlcikuLk7dunWjDS2KNrQ22s/6aEPr83cbZo1YweUjgOQiKChILVq00KJFi9SnTx9Jksfj0aJFizR06NB8PYfb7dbmzZvVs2fPPM9xOp1yOnOuPhUYGOi3X4KhZdP8ej0UDtrQ+mhDa6P9rI82tD5/tSHfJ1eOAJKHESNGqH///mrZsqVat26tyZMnKzk5WQMHDpQk9evXT9HR0ZowYYIk6YUXXtA111yjOnXq6NSpU5o4caL279+vhx9+2MyXAQAAABQpBJA83HXXXTp27JhGjx6t+Ph4xcbG6vvvv/dNTD9w4IDs9r8WETt58qQGDRqk+Ph4lSlTRi1atNDy5cvVqFEjs14CAAAAUOQQQC5g6NCheQ65WrJkSbbt119/Xa+//rofqgIAAACsi/uAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvwkwuwD8xTAMSVJSUpJfrudyuZSSkqKkpCQFBgb65ZooWLSh9dGG1kb7WR9taH3+bsOs92lZ79tw6QggRcjp06clSdWqVTO5EgAAAFzI6dOnVapUKbPLsCSbQXwrMjwejw4fPqyIiAjZbLZCv15SUpKqVaumP/74Q5GRkYV+PRQ82tD6aENro/2sjza0Pn+3oWEYOn36tKpUqSK7ndkMl4MekCLEbreratWqfr9uZGQkv3Qtjja0PtrQ2mg/66MNrc+fbUjPx5UhtgEAAADwGwIIAAAAAL8hgJRgTqdTY8aMkdPpNLsUXCba0PpoQ2uj/ayPNrQ+2tB6mIQOAAAAwG/oAQEAAADgNwQQAAAAAH5DAAEAAADgNwQQAAAAAH5DACnBpkyZopiYGAUHB6tNmzZavXq12SUhnyZMmKBWrVopIiJCFStWVJ8+fbRjxw6zy8Jl+r//+z/ZbDYNGzbM7FJwCQ4dOqT7779f5cqVU0hIiJo2baq1a9eaXRbyye12a9SoUapZs6ZCQkJUu3Ztvfjii2JtnqLrp59+Uu/evVWlShXZbDbNnTs323HDMDR69GhVrlxZISEh6tq1q3bt2mVOsbggAkgJNWfOHI0YMUJjxozR+vXr1axZM/Xo0UNHjx41uzTkw9KlSzVkyBCtXLlScXFxcrlc6t69u5KTk80uDZdozZo1mjZtmq666iqzS8ElOHnypK699loFBgbqu+++09atW/Xaa6+pTJkyZpeGfHr55Zf19ttv69///re2bduml19+Wa+88oreeusts0tDHpKTk9WsWTNNmTIl1+OvvPKK3nzzTU2dOlWrVq1SWFiYevToobS0ND9XiothGd4Sqk2bNmrVqpX+/e9/S5I8Ho+qVaumv//97xo5cqTJ1eFSHTt2TBUrVtTSpUvVsWNHs8tBPp05c0ZXX321/vOf/+ill15SbGysJk+ebHZZyIeRI0fql19+0c8//2x2KbhMN910k6KiojR9+nTfvttvv10hISH66KOPTKwM+WGz2fTll1+qT58+kry9H1WqVNETTzyhJ598UpKUmJioqKgozZw5U3fffbeJ1eJ89ICUQBkZGVq3bp26du3q22e329W1a1etWLHCxMpwuRITEyVJZcuWNbkSXIohQ4aoV69e2X4WYQ3z5s1Ty5Ytdeedd6pixYpq3ry53n33XbPLwiVo166dFi1apJ07d0qSNm3apGXLlunGG280uTJcjr179yo+Pj7b79NSpUqpTZs2vLcpggLMLgD+l5CQILfbraioqGz7o6KitH37dpOqwuXyeDwaNmyYrr32WjVp0sTscpBP//vf/7R+/XqtWbPG7FJwGfbs2aO3335bI0aM0LPPPqs1a9boH//4h4KCgtS/f3+zy0M+jBw5UklJSWrQoIEcDofcbrfGjRun++67z+zScBni4+MlKdf3NlnHUHQQQACLGzJkiLZs2aJly5aZXQry6Y8//tDjjz+uuLg4BQcHm10OLoPH41HLli01fvx4SVLz5s21ZcsWTZ06lQBiEZ988ok+/vhjzZ49W40bN9bGjRs1bNgwValShTYEChlDsEqg8uXLy+Fw6MiRI9n2HzlyRJUqVTKpKlyOoUOH6ptvvtHixYtVtWpVs8tBPq1bt05Hjx7V1VdfrYCAAAUEBGjp0qV68803FRAQILfbbXaJuIjKlSurUaNG2fY1bNhQBw4cMKkiXKqnnnpKI0eO1N13362mTZvqgQce0PDhwzVhwgSzS8NlyHr/wnsbayCAlEBBQUFq0aKFFi1a5Nvn8Xi0aNEitW3b1sTKkF+GYWjo0KH68ssv9eOPP6pmzZpml4RLcP3112vz5s3auHGj76Nly5a67777tHHjRjkcDrNLxEVce+21OZa+3rlzp2rUqGFSRbhUKSkpstuzvw1yOBzyeDwmVYQrUbNmTVWqVCnbe5ukpCStWrWK9zZFEEOwSqgRI0aof//+atmypVq3bq3JkycrOTlZAwcONLs05MOQIUM0e/ZsffXVV4qIiPCNby1VqpRCQkJMrg4XExERkWO+TlhYmMqVK8c8HosYPny42rVrp/Hjx6tv375avXq13nnnHb3zzjtml4Z86t27t8aNG6fq1aurcePG2rBhgyZNmqQHH3zQ7NKQhzNnzmj37t2+7b1792rjxo0qW7asqlevrmHDhumll15S3bp1VbNmTY0aNUpVqlTxrZSFIsRAifXWW28Z1atXN4KCgozWrVsbK1euNLsk5JOkXD9mzJhhdmm4TNddd53x+OOPm10GLsHXX39tNGnSxHA6nUaDBg2Md955x+yScAmSkpKMxx9/3KhevboRHBxs1KpVy/jXv/5lpKenm10a8rB48eJc/+/r37+/YRiG4fF4jFGjRhlRUVGG0+k0rr/+emPHjh3mFo1ccR8QAAAAAH7DHBAAAAAAfkMAAQAAAOA3BBAAAAAAfkMAAQAAAOA3BBAAAAAAfkMAAQAAAOA3BBAAAAAAfkMAAQAAAOA3BBAAKIFmzpwpm82W58fKlSvNLhEAUEwFmF0AAMA8L7zwgmrWrJljf506dUyoBgBQEhBAAKAEu/HGG9WyZUtTa0hOTlZYWJipNQAA/IchWACAXO3bt082m02vvvqq3nnnHdWuXVtOp1OtWrXSmjVrcpy/fft23XHHHSpbtqyCg4PVsmVLzZs3L9s5WUO/li5dqsGDB6tixYqqWrWq7/iUKVNUq1YthYSEqHXr1vr555/VqVMnderUSZJ05swZhYWF6fHHH89x/YMHD8rhcGjChAkF+4UAABQoekAAoARLTExUQkJCtn02m03lypXzbc+ePVunT5/Wo48+KpvNpldeeUW33Xab9uzZo8DAQEnSb7/9pmuvvVbR0dEaOXKkwsLC9Mknn6hPnz76/PPPdeutt2a7xuDBg1WhQgWNHj1aycnJkqS3335bQ4cOVYcOHTR8+HDt27dPffr0UZkyZXwhJTw8XLfeeqvmzJmjSZMmyeFw+J7zv//9rwzD0H333VcoXysAQMEggABACda1a9cc+5xOp9LS0nzbBw4c0K5du1SmTBlJUv369XXLLbdowYIFuummmyRJjz/+uKpXr641a9bI6XRK8oaM9u3b65///GeOAFK2bFktWrTIFyAyMjI0atQotWrVSj/++KMCArz/PV111VUaMGBAtl6Sfv366eOPP1ZcXJxuuOEG3/6PPvpIHTt2VPXq1QviSwMAKCQMwQKAEmzKlCmKi4vL9vHdd99lO+euu+7yhQ9J6tChgyRpz549kqQTJ07oxx9/VN++fXX69GklJCQoISFBx48fV48ePbRr1y4dOnQo23MOGjQoW+/F2rVrdfz4cQ0aNMgXPiTpvvvuy3ZtyRuaqlSpoo8//ti3b8uWLfr11191//33X+FXBABQ2OgBAYASrHXr1hedhH5+j0JWIDh58qQkaffu3TIMQ6NGjdKoUaNyfY6jR48qOjrat33+ylv79++XlHP1rYCAAMXExGTbZ7fbdd999+ntt99WSkqKQkND9fHHHys4OFh33nnnBV8LAMB8BBAAwAWd21NxLsMwJEkej0eS9OSTT6pHjx65nnt+sAgJCbmimvr166eJEydq7ty5uueeezR79mzddNNNKlWq1BU9LwCg8BFAAABXpFatWpKkwMDAXOeU5EeNGjUkeXtTOnfu7NufmZmpffv26aqrrsp2fpMmTdS8eXN9/PHHqlq1qg4cOKC33nrrMl8BAMCfmAMCALgiFStWVKdOnTRt2jT9+eefOY4fO3bsos/RsmVLlStXTu+++64yMzN9+z/++GPfUK/zPfDAA1q4cKEmT56scuXK6cYbb7z8FwEA8Bt6QACgBPvuu++0ffv2HPvbtWsnuz3/f6OaMmWK2rdvr6ZNm2rQoEGqVauWjhw5ohUrVujgwYPatGnTBR8fFBSk559/Xn//+9/VpUsX9e3bV/v27dPMmTNVu3Zt2Wy2HI+599579fTTT+vLL7/UY4895lsSGABQtBFAAKAEGz16dK77Z8yY4bv5X340atRIa9eu1dixYzVz5kwdP35cFStWVPPmzfO8xvmGDh0qwzD02muv6cknn1SzZs00b948/eMf/1BwcHCO86OiotS9e3fNnz9fDzzwQL5rBQCYy2ZkzSIEAKCI8Xg8qlChgm677Ta9++67OY7feuut2rx5s3bv3m1CdQCAy8EcEABAkZCWlqbz/yb24Ycf6sSJE7n2xvz555/69ttv6f0AAIuhBwQAUCQsWbJEw4cP15133qly5cpp/fr1mj59uho2bKh169YpKChIkrR371798ssveu+997RmzRr9/vvvqlSpksnVAwDyizkgAIAiISYmRtWqVdObb76pEydOqGzZsurXr5/+7//+zxc+JGnp0qUaOHCgqlevrg8++IDwAQAWQw8IAAAAAL9hDggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPAbAggAAAAAvyGAAAAAAPCb/wdI6TOgAgVidAAAAABJRU5ErkJggg==", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "make_toxin_widget(\n", + " alpha=(0.1, 10.0),\n", + " t0=(-10, 10, False),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6de1fd60-7356-4151-be34-695721fd52ba", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2d514ca-e92c-4198-84ed-942432d80fe1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca152029-a095-4607-a368-8ec4d6886f8d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "emevo-lab", + "language": "python", + "name": "emevo-lab" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 9eaa378..df233d9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -94,5 +94,4 @@ dev-dependencies = [ "seaborn >= 0.12", "typer >= 0.12", "tqdm >= 4.6", - "pillow>=10.4.0", ] diff --git a/src/emevo/env.py b/src/emevo/env.py index ab51d0b..ce110b1 100644 --- a/src/emevo/env.py +++ b/src/emevo/env.py @@ -52,34 +52,6 @@ def update(self, energy_delta: jax.Array, capacity: float | None = 100.0) -> Sel return replace(self, energy=jnp.clip(energy, max=capacity)) -@chex.dataclass -class StatusWithToxin(Status): - toxin: jax.Array - - def deactivate(self, flag: jax.Array) -> Self: - return replace( - self, - age=jnp.where(flag, 0, self.age), - toxin=jnp.where(flag, 0, self.toxin), - ) - - def update( - self, - energy_delta: jax.Array, - toxin_delta: jax.Array, - capacity: float | None = 100.0, - toxin_capacity: float | None = 10.0, - ) -> Self: - """Update energy and toxin.""" - energy = self.energy + energy_delta - toxin = self.toxin + toxin_delta - return replace( - self, - energy=jnp.clip(energy, max=capacity), - toxin=jnp.clip(toxin, max=toxin_capacity), - ) - - def init_status(max_n: int, init_energy: float) -> Status: return Status( age=jnp.zeros(max_n, dtype=jnp.int32), diff --git a/src/emevo/environments/__init__.py b/src/emevo/environments/__init__.py index 63c523e..04a4fc7 100644 --- a/src/emevo/environments/__init__.py +++ b/src/emevo/environments/__init__.py @@ -11,6 +11,6 @@ register( "CircleForaging-v1", - "emevo.environments.circle_foraging_with.CircleForaging", - "Phyjax2d circle foraging environment", + "emevo.environments.circle_foraging_with_neurotoxin.CircleForagingWithNeurotoxin", + "Phyjax2d circle foraging environment with neuro toxin", ) diff --git a/src/emevo/environments/circle_foraging.py b/src/emevo/environments/circle_foraging.py index 6d2cef7..54bab6d 100644 --- a/src/emevo/environments/circle_foraging.py +++ b/src/emevo/environments/circle_foraging.py @@ -6,7 +6,7 @@ import warnings from collections.abc import Callable, Iterable from dataclasses import replace -from typing import Any, Literal, NamedTuple +from typing import Any, Generic, Literal, NamedTuple, TypeVar import chex import jax @@ -87,8 +87,11 @@ def as_array(self) -> jax.Array: ) +S = TypeVar("S", bound=Status) + + @chex.dataclass -class CFState: +class CFState(Generic[S]): physics: StateDict solver: VelocitySolver food_num: list[FoodNumState] @@ -97,7 +100,7 @@ class CFState: key: chex.PRNGKey step: jax.Array unique_id: UniqueID - status: Status + status: S n_born_agents: jax.Array @property @@ -407,7 +410,7 @@ def _set_b2a( def _make_food_energy_coef_array( food_energy_coef: Iterable[float | tuple[float, ...]], -) -> tuple[bool,]: +) -> jax.Array: has_tuple = any([isinstance(fec, tuple) for fec in food_energy_coef]) if has_tuple: length = [len(fec) if isinstance(fec, tuple) else 1 for fec in food_energy_coef] @@ -457,7 +460,7 @@ def __init__( agent_radius: float = 10.0, food_radius: float = 4.0, foodloc_interval: int = 1000, - fec_intervals: tuple[int, ...] = 1, + fec_intervals: tuple[int, ...] = (1,), dt: float = 0.1, linear_damping: float = 0.8, angular_damping: float = 0.6, @@ -501,12 +504,10 @@ def __init__( self._fec_intervals = jnp.array(fec_intervals, dtype=jnp.int32) self._food_num_fns, self._initial_foodnum_states = [], [] if n_food_sources > 1: - assert isinstance(food_loc_fn, list | tuple) and n_food_sources == len( - food_loc_fn - ) - assert isinstance(food_num_fn, list | tuple) and n_food_sources == len( - food_num_fn - ) + assert isinstance(food_loc_fn, list | tuple) + assert n_food_sources == len(food_loc_fn) + assert isinstance(food_num_fn, list | tuple) + assert n_food_sources == len(food_num_fn) else: food_loc_fn, food_num_fn = [food_loc_fn], [food_num_fn] # type: ignore for maybe_loc_fn in food_loc_fn: # type: ignore @@ -812,9 +813,9 @@ def _get_selected_sensor( def step( self, - state: CFState, + state: CFState[Status], action: ArrayLike, - ) -> tuple[CFState, TimeStep[CFObs]]: + ) -> tuple[CFState[Status], TimeStep[CFObs]]: # Add force act = jax.vmap(self.act_space.clip)(jnp.array(action)) f1_raw = jax.lax.slice_in_dim(act, 0, 1, axis=-1) @@ -921,7 +922,7 @@ def step( def activate( self, - state: CFState, + state: CFState[Status], is_parent: jax.Array, ) -> tuple[CFState, jax.Array]: N = self.n_max_agents @@ -985,7 +986,7 @@ def activate( ) return new_state, parent_id - def deactivate(self, state: CFState, flag: jax.Array) -> CFState: + def deactivate(self, state: CFState, flag: jax.Array) -> CFState[Status]: expanded_flag = jnp.expand_dims(flag, axis=1) p_xy = jnp.where(expanded_flag, NOWHERE, state.physics.circle.p.xy) p = replace(state.physics.circle.p, xy=p_xy) @@ -999,7 +1000,7 @@ def deactivate(self, state: CFState, flag: jax.Array) -> CFState: status = state.status.deactivate(flag) return replace(state, physics=physics, unique_id=unique_id, status=status) - def reset(self, key: chex.PRNGKey) -> tuple[CFState, TimeStep[CFObs]]: + def reset(self, key: chex.PRNGKey) -> tuple[CFState[Status], TimeStep[CFObs]]: physics, agent_loc, food_loc, food_num = self._initialize_physics_state(key) N = self.n_max_agents unique_id = init_uniqueid(self._n_initial_agents, N) @@ -1175,7 +1176,7 @@ def _remove_and_regenerate_foods( def visualizer( self, - state: CFState, + state: CFState[Status], figsize: tuple[float, float] | None = None, sensor_index: int | None = None, backend: str = "pyglet", @@ -1193,7 +1194,7 @@ def visualizer( figsize=figsize, backend=backend, sensor_fn=( - self._get_sensors + self._get_sensors # type: ignore if sensor_index is None else lambda stated: self._get_selected_sensor(stated, sensor_index) ), diff --git a/src/emevo/environments/circle_foraging_with_neurotoxin.py b/src/emevo/environments/circle_foraging_with_neurotoxin.py new file mode 100644 index 0000000..e2c8ff5 --- /dev/null +++ b/src/emevo/environments/circle_foraging_with_neurotoxin.py @@ -0,0 +1,215 @@ +"""Various hand-coded rules for the effect of food""" + +from dataclasses import replace +from typing import Any + +import chex +import jax +import jax.numpy as jnp +from jax.typing import ArrayLike + +from emevo.env import Status, TimeStep +from emevo.environments.circle_foraging import ( + CFObs, + CFState, + CircleForaging, + get_tactile, + init_uniqueid, + nstep, +) + +Self = Any + + +@chex.dataclass +class StatusWithToxin(Status): + toxin: jax.Array + + def deactivate(self, flag: jax.Array) -> Self: + return replace( + self, + age=jnp.where(flag, 0, self.age), + toxin=jnp.where(flag, 0, self.toxin), + ) + + +def init_status(max_n: int, init_energy: float) -> StatusWithToxin: + return StatusWithToxin( + age=jnp.zeros(max_n, dtype=jnp.int32), + energy=jnp.ones(max_n, dtype=jnp.float32) * init_energy, + toxin=jnp.zeros(max_n, dtype=jnp.float32), + ) + + +class CircleForagingWithNeurotoxin(CircleForaging): + def __init__( + self, + *args, + toxin_t0: float = 5.0, + toxin_alpha: float = 1.0, + toxin_delta: float = 10.0, + toxin_decay: float = 0.1, + **kwargs, + ) -> None: + super().__init__(*args, **kwargs) + self._toxin_t0 = toxin_t0 + self._toxin_alpha = toxin_alpha + self._toxin_decay = toxin_decay + self._toxin_delta = toxin_delta + assert self._n_food_sources - 1 == self._food_energy_coef.shape[1] + + def step( # type: ignore + self, + state: CFState[StatusWithToxin], + action: ArrayLike, + ) -> tuple[CFState, TimeStep[CFObs]]: + # Compute action decay ratio by toxin + toxin_decay_rate = 1.0 / ( + 1.0 + self._toxin_alpha * jnp.exp(self._toxin_t0 - state.status.toxin) + ) + toxin_decay = jnp.expand_dims(1.0 - toxin_decay_rate, axis=1) + # Add force + act = jax.vmap(self.act_space.clip)(jnp.array(action)) * toxin_decay + f1_raw = jax.lax.slice_in_dim(act, 0, 1, axis=-1) + f2_raw = jax.lax.slice_in_dim(act, 1, 2, axis=-1) + f1 = jnp.concatenate((jnp.zeros_like(f1_raw), f1_raw), axis=1) + f2 = jnp.concatenate((jnp.zeros_like(f2_raw), f2_raw), axis=1) + circle = state.physics.circle + circle = circle.apply_force_local(self._act_p1, f1) + circle = circle.apply_force_local(self._act_p2, f2) + stated = replace(state.physics, circle=circle) + # Step physics simulator + stated, solver, nstep_contacts = nstep( + self._n_physics_iter, + self._physics, + stated, + state.solver, + ) + # Gather circle contacts + contacts = jnp.max(nstep_contacts, axis=0) + c2c = self._physics.get_contact_mat("circle", "circle", contacts) + c2sc = self._physics.get_contact_mat("circle", "static_circle", contacts) + seg2c = self._physics.get_contact_mat("segment", "circle", contacts) + # Get tactile obs + food_tactile, ft_raw = self._food_tactile( + stated.static_circle.label, + stated.circle, + stated.static_circle, + c2sc, + ) + ag_tactile, _ = get_tactile( + self._n_tactile_bins, + stated.circle, + stated.circle, + c2c, + ) + wall_tactile, _ = get_tactile( + self._n_tactile_bins, + stated.circle, + stated.segment, + seg2c.transpose(), + ) + collision = jnp.concatenate( + (ag_tactile > 0, food_tactile > 0, wall_tactile > 0), + axis=1, + ) + # Gather sensor obs + sensor_obs = self._sensor_obs(stated=stated) + # energy_delta = food - coef * |force| + force_norm = jnp.sqrt(f1_raw**2 + f2_raw**2).ravel() + energy_consumption = ( + self._force_energy_consumption * force_norm + self._basic_energy_consumption + ) + n_ate = jnp.sum(food_tactile[:, :, self._foraging_indices], axis=-1) + # toxin and foods + n_ate_foods = n_ate[:, : self._n_food_sources - 1] # (N-agents, N-foods) + n_ate_toxin = n_ate[:, self._n_food_sources - 1] # (N-agents,) + energy_gain = jnp.sum(n_ate_foods * self._food_energy_coef, axis=1) + energy_delta = energy_gain - energy_consumption + # Remove and regenerate foods + key, food_key = jax.random.split(state.key) + eaten = jnp.sum(ft_raw[:, :, :, self._foraging_indices], axis=(0, 3)) > 0 + stated, food_num, food_loc, n_regen = self._remove_and_regenerate_foods( + food_key, + eaten, # (N_FOOD, N_LABEL) + stated, + state.step, + state.food_num, + state.food_loc, + ) + status = state.status.update( + energy_delta=energy_delta, + capacity=self._energy_capacity, + ) + toxin = jnp.clip( + status.toxin + n_ate_toxin * self._toxin_delta - self._toxin_decay, + min=0.0, + ) + status = replace(status, toxin=toxin) + # Construct obs + obs = CFObs( + sensor=sensor_obs.reshape(-1, self._n_sensors, self._n_obj), + collision=collision, + angle=stated.circle.p.angle, + velocity=stated.circle.v.xy, + angular_velocity=stated.circle.v.angle, + energy=status.energy, + ) + timestep = TimeStep( + encount=c2c, + obs=obs, + info={ + "energy_gain": energy_gain, + "energy_consumption": energy_consumption, + "n_food_regenerated": n_regen, + "n_food_eaten": jnp.sum(eaten, axis=0), # (N_LABEL,) + "n_ate_food": n_ate_foods, # (N_AGENT, N_LABEL - 1) + "n_ate_toxin": n_ate_toxin, # (N_AGENT,) + }, + ) + state = CFState( + physics=stated, + solver=solver, + food_num=food_num, + agent_loc=state.agent_loc, + food_loc=food_loc, + key=key, + step=state.step + 1, + unique_id=state.unique_id, + status=status.step(state.unique_id.is_active()), + n_born_agents=state.n_born_agents, + ) + return state, timestep + + def reset( # type: ignore + self, + key: chex.PRNGKey, + ) -> tuple[CFState[StatusWithToxin], TimeStep[CFObs]]: + physics, agent_loc, food_loc, food_num = self._initialize_physics_state(key) + N = self.n_max_agents + unique_id = init_uniqueid(self._n_initial_agents, N) + status = init_status(N, self._init_energy) + state = CFState( + physics=physics, + solver=self._physics.init_solver(), + agent_loc=agent_loc, + food_loc=food_loc, + food_num=food_num, + key=key, + step=jnp.array(0, dtype=jnp.int32), + unique_id=unique_id, + status=status, + n_born_agents=jnp.array(self._n_initial_agents, dtype=jnp.int32), + ) + sensor_obs = self._sensor_obs(stated=physics) + obs = CFObs( + sensor=sensor_obs.reshape(-1, self._n_sensors, self._n_obj), + collision=jnp.zeros((N, self._n_obj, self._n_tactile_bins), dtype=bool), + angle=physics.circle.p.angle, + velocity=physics.circle.v.xy, + angular_velocity=physics.circle.v.angle, + energy=state.status.energy, + ) + # They shouldn't encount now + timestep = TimeStep(encount=jnp.zeros((N, N), dtype=bool), obs=obs) + return state, timestep diff --git a/src/emevo/exp_utils.py b/src/emevo/exp_utils.py index e5b103d..6a5bab9 100644 --- a/src/emevo/exp_utils.py +++ b/src/emevo/exp_utils.py @@ -165,6 +165,7 @@ class Log: parents: jax.Array rewards: jax.Array unique_id: jax.Array + additional_fields: dict[str, jax.Array] = dataclasses.field(default_factory=dict) def with_step(self, from_: int) -> LogWithStep: if self.parents.ndim == 2: @@ -214,9 +215,10 @@ def filter_death(self) -> Any: def to_flat_dict(self) -> dict[str, jax.Array]: d = dataclasses.asdict(self.log) # type: ignore + additional = d.pop("additional_fields") d["step"] = self.step d["slots"] = self.slots - return d + return d | additional @dataclasses.dataclass diff --git a/tests/test_tree.py b/tests/test_tree.py index 092b539..40b2c82 100644 --- a/tests/test_tree.py +++ b/tests/test_tree.py @@ -60,7 +60,7 @@ def test_from_iter(treedef: list[tuple[int, int]]) -> None: assert preorder == [0, 1, 3, 4, 5, 8, 9, 2, 6, 7] postorder = list(map(operator.attrgetter("index"), tree.traverse(preorder=False))) assert postorder == [3, 4, 8, 9, 5, 1, 6, 7, 2, 0] - assert tree.root.n_total_children == 10 + assert tree.root.n_descendants == 10 def test_split(treedef: list[tuple[int, int]]) -> None: