From c9ed1730f8a73ebbc08f89d056b9c1f15aa4963d Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 4 May 2023 16:00:49 -0700 Subject: [PATCH 1/5] initial ReaddySimpleProcess runs, todo viz and test output --- env.yml | 4 +- .../processes/readdy_actin_process.py | 50 +--- .../processes/readdy_simple_process.py | 231 ++++++++++++++++++ vivarium_models/util.py | 49 ++++ 4 files changed, 284 insertions(+), 50 deletions(-) create mode 100644 vivarium_models/processes/readdy_simple_process.py diff --git a/env.yml b/env.yml index fb055f5..bd90a8e 100644 --- a/env.yml +++ b/env.yml @@ -3,8 +3,8 @@ name: vivarium-models channels: - conda-forge dependencies: - - python=3.8 - - conda-forge::readdy + - python=3.9 + - readdy==2.0.9 - pip - pip: - "vivarium_models @ git+https://github.com/simularium/vivarium-models.git" diff --git a/vivarium_models/processes/readdy_actin_process.py b/vivarium_models/processes/readdy_actin_process.py index ba92cd1..a7ceadd 100644 --- a/vivarium_models/processes/readdy_actin_process.py +++ b/vivarium_models/processes/readdy_actin_process.py @@ -12,7 +12,7 @@ ActinAnalyzer, ) from simularium_readdy_models import ReaddyUtil -from vivarium_models.util import create_monomer_update, format_monomer_results +from vivarium_models.util import create_monomer_update, format_monomer_results, monomer_ports_schema from vivarium_models.library.scan import Scan NAME = "ReaDDy_actin" @@ -69,53 +69,7 @@ def __init__(self, parameters=None): self.readdy_simulation = actin_simulation.simulation def ports_schema(self): - return { - "monomers": { - "box_center": { - "_default": np.array([1000.0, 0.0, 0.0]), - "_updater": "set", - "_emit": True, - }, - "box_size": { - "_default": 500.0, - "_updater": "set", - "_emit": True, - }, - "topologies": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "particle_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - "particles": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "position": { - "_default": np.zeros(3), - "_updater": "set", - "_emit": True, - }, - "neighbor_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - } - } + return monomer_ports_schema def initial_state(self, config=None): # TODO: make this more general diff --git a/vivarium_models/processes/readdy_simple_process.py b/vivarium_models/processes/readdy_simple_process.py new file mode 100644 index 0000000..9eddffc --- /dev/null +++ b/vivarium_models/processes/readdy_simple_process.py @@ -0,0 +1,231 @@ +import os +import string + +from vivarium.core.process import Process +from vivarium.core.directories import PROCESS_OUT_DIR +from vivarium.core.engine import Engine + +import numpy as np +import readdy +from tqdm import tqdm + +from subcell_analysis.readdy import ReaddyUtil +from vivarium_models.util import create_monomer_update, monomer_ports_schema + +NAME = "ReaDDy_simple" + + +class ReaddySimpleProcess(Process): + """ + This process uses ReaDDy to model 4 diffusing particles. + """ + + name = NAME + + defaults = { + "n_particles" : 4, + "box_size" : 500., # nm + "temperature_C" : 22.0, + "viscosity" : 8.1, # cP + "periodic_boundary" : False, + "force_constant" : 250., + "internal_timestep" : 0.1, # ns + } + + + def __init__(self, parameters=None): + super(ReaddySimpleProcess, self).__init__(parameters) + self.create_system() + + + def calculate_diffusionCoefficient(self, radius): + """ + calculates the theoretical diffusion constant of a spherical particle + in a media with viscosity [cP] at temperature [Kelvin]. + with radius radius[nm] + + returns nm^2/s + """ + return ( + (1.38065 * 10 ** (-23) * self.parameters["temperature_K"]) + / (6 * np.pi * self.parameters["viscosity"] * 10 ** (-3) * radius * 10 ** (-9)) + * 10**18 + / 10**9 + ) + + + def create_system(self): + """ + Create a ReaDDy system with the given number of + types of particles freely diffusing in a box. + """ + # create readdy system + box_size_3d = np.array([self.parameters["box_size"]] * 3) + self.readdy_system = readdy.ReactionDiffusionSystem( + box_size=box_size_3d, + periodic_boundary_conditions=[self.parameters["periodic_boundary"]] * 3, + ) + self.parameters["temperature_K"] = self.parameters["temperature_C"] + 273.15 + self.readdy_system.temperature = self.parameters["temperature_K"] + + # add particle species + letters = list(string.ascii_uppercase) + for p_ix in range(self.parameters["n_particles"]): + diffCoeff = self.calculate_diffusionCoefficient(1. + 0.5 * p_ix) # nm^2/s + self.readdy_system.add_topology_species(letters[p_ix], diffCoeff) + + # add repulsion potentials for excluded volume + for p_ix1 in range(self.parameters["n_particles"]): + for p_ix2 in range(p_ix1 + 1, self.parameters["n_particles"]): + self.readdy_system.potentials.add_harmonic_repulsion( + particle_type1=letters[p_ix1], + particle_type2=letters[p_ix2], + force_constant=self.parameters["force_constant"], + interaction_distance=1.5 + 0.5 * (p_ix1 + p_ix2), + ) + + # add box potential if no periodic boundary + if not self.parameters["periodic_boundary"]: + box_potential_size = box_size_3d - 2.0 + for particle_type in letters[:self.parameters["n_particles"]]: + self.readdy_system.potentials.add_box( + particle_type=particle_type, + force_constant=self.parameters["force_constant"], + origin=-0.5 * box_potential_size, + extent=box_potential_size, + ) + + + def ports_schema(self): + return monomer_ports_schema + + + def create_simulation_with_particles(self, particles): + self.readdy_simulation = self.readdy_system.simulation("CPU") + for p_id in particles: + self.readdy_simulation.add_particle( + type=particles[p_id]["type_name"], + position=particles[p_id]["position"], + ) + + def simulate_readdy(self, timestep): + """ + Simulate in ReaDDy for the given timestep + """ + def loop(): + readdy_actions = self.readdy_simulation._actions + init = readdy_actions.initialize_kernel() + diffuse = readdy_actions.integrator_euler_brownian_dynamics( + self.parameters["internal_timestep"] + ) + calculate_forces = readdy_actions.calculate_forces() + create_nl = readdy_actions.create_neighbor_list( + self.readdy_system.calculate_max_cutoff().magnitude + ) + update_nl = readdy_actions.update_neighbor_list() + react = readdy_actions.reaction_handler_uncontrolled_approximation( + self.parameters["internal_timestep"] + ) + observe = readdy_actions.evaluate_observables() + init() + create_nl() + calculate_forces() + update_nl() + observe(0) + n_steps = int(timestep / self.parameters["internal_timestep"]) + for t in tqdm(range(1, n_steps + 1)): + diffuse() + update_nl() + react() + update_nl() + calculate_forces() + observe(t) + + self.readdy_simulation._run_custom_loop(loop, show_summary=False) + + + def get_particles_from_simulation(self): + result = { + "box_center" : np.zeros(3), + "box_size" : self.parameters["box_size"], + "topologies": {}, + "particles": {}, + } + for particle in self.readdy_simulation.current_particles: + result["particles"][str(particle.id)] = { + "type_name" : particle.type, + "position" : particle.pos, + } + return result + + + def next_update(self, timestep, states): + print(f"ReaDDy simple process updating by {timestep} ns") + + self.create_simulation_with_particles(states["monomers"]["particles"]) + self.simulate_readdy(timestep) + new_monomers = self.get_particles_from_simulation() + + return create_monomer_update(states["monomers"], new_monomers) + + +def random_initial_state(parameters): + result = { + "box_center" : np.zeros(3), + "box_size" : parameters["box_size"], + "topologies": {}, + "particles": {}, + } + letters = list(string.ascii_uppercase) + random_positions = ( + parameters["box_size"] * + (np.random.uniform(size=(parameters["n_particles"], 3)) - 0.5) + ) + for p_ix in range(parameters["n_particles"]): + result["particles"][str(p_ix)] = { + "type_name" : letters[p_ix], + "position" : random_positions[p_ix], + } + return {"monomers" : result} + + +# functions to configure and run the process +def run_readdy_simple_process(): + """ + Run a simulation of the process. + + Returns: + The simulation output. + """ + readdy_process = ReaddySimpleProcess() + composite = readdy_process.generate() + initial_state = random_initial_state(readdy_process.parameters) + + sim = Engine(composite=composite, initial_state=initial_state) + + sim.update(10) # ns + + output = sim.emitter.get_timeseries() + return output + + +def test_readdy_simple_process(): + ''' + Test that the process runs correctly. + This will be executed by pytest. + ''' + output = run_readdy_simple_process() + # TODO: Add assert statements to ensure correct performance. + + +def main(): + '''Simulate the process and plot results.''' + # make an output directory to save results + out_dir = os.path.join(PROCESS_OUT_DIR, NAME) + os.makedirs(out_dir, exist_ok=True) + + output = run_readdy_simple_process() + + +if __name__ == '__main__': + main() diff --git a/vivarium_models/util.py b/vivarium_models/util.py index 9be67e5..62c04bd 100644 --- a/vivarium_models/util.py +++ b/vivarium_models/util.py @@ -1,6 +1,55 @@ import numpy as np +monomer_ports_schema = { + "monomers": { + "box_center": { + "_default": np.array([1000.0, 0.0, 0.0]), + "_updater": "set", + "_emit": True, + }, + "box_size": { + "_default": 500.0, + "_updater": "set", + "_emit": True, + }, + "topologies": { + "*": { + "type_name": { + "_default": "", + "_updater": "set", + "_emit": True, + }, + "particle_ids": { + "_default": [], + "_updater": "set", + "_emit": True, + }, + } + }, + "particles": { + "*": { + "type_name": { + "_default": "", + "_updater": "set", + "_emit": True, + }, + "position": { + "_default": np.zeros(3), + "_updater": "set", + "_emit": True, + }, + "neighbor_ids": { + "_default": [], + "_updater": "set", + "_emit": True, + }, + } + }, + } +} + + def agents_update(existing, projected): update = {"_add": [], "_delete": []} From 17477e89aa076973fc922a09d5e35afa7d1907ab Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 4 May 2023 17:51:13 -0700 Subject: [PATCH 2/5] added radius to monomer ports schema --- vivarium_models/util.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/vivarium_models/util.py b/vivarium_models/util.py index 62c04bd..c928614 100644 --- a/vivarium_models/util.py +++ b/vivarium_models/util.py @@ -44,6 +44,11 @@ "_updater": "set", "_emit": True, }, + "radius": { + "_default": 1.0, + "_updater": "set", + "_emit": True, + }, } }, } From ac2223652d777a82cc51a49994eab1d64af92266 Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 4 May 2023 17:51:51 -0700 Subject: [PATCH 3/5] viz simple readdy output and fixed issues --- .../processes/readdy_simple_process.py | 40 +++++-- .../processes/simularium_emitter.py | 102 ++++-------------- 2 files changed, 53 insertions(+), 89 deletions(-) diff --git a/vivarium_models/processes/readdy_simple_process.py b/vivarium_models/processes/readdy_simple_process.py index 9eddffc..70cf2eb 100644 --- a/vivarium_models/processes/readdy_simple_process.py +++ b/vivarium_models/processes/readdy_simple_process.py @@ -9,7 +9,6 @@ import readdy from tqdm import tqdm -from subcell_analysis.readdy import ReaddyUtil from vivarium_models.util import create_monomer_update, monomer_ports_schema NAME = "ReaDDy_simple" @@ -29,7 +28,8 @@ class ReaddySimpleProcess(Process): "viscosity" : 8.1, # cP "periodic_boundary" : False, "force_constant" : 250., - "internal_timestep" : 0.1, # ns + "internal_timestep" : 1000., # ns + "time_step" : 10000., # ns } @@ -38,6 +38,11 @@ def __init__(self, parameters=None): self.create_system() + @staticmethod + def radius_for_ix(ix: int): + return 10. + 5. * ix + + def calculate_diffusionCoefficient(self, radius): """ calculates the theoretical diffusion constant of a spherical particle @@ -71,17 +76,24 @@ def create_system(self): # add particle species letters = list(string.ascii_uppercase) for p_ix in range(self.parameters["n_particles"]): - diffCoeff = self.calculate_diffusionCoefficient(1. + 0.5 * p_ix) # nm^2/s - self.readdy_system.add_topology_species(letters[p_ix], diffCoeff) + particle_type = letters[p_ix] + radius = ReaddySimpleProcess.radius_for_ix(p_ix) + diffCoeff = self.calculate_diffusionCoefficient(radius) # nm^2/s + self.readdy_system.add_topology_species(particle_type, diffCoeff) # add repulsion potentials for excluded volume for p_ix1 in range(self.parameters["n_particles"]): for p_ix2 in range(p_ix1 + 1, self.parameters["n_particles"]): + particle_type1 = letters[p_ix1] + particle_type2 = letters[p_ix2] self.readdy_system.potentials.add_harmonic_repulsion( - particle_type1=letters[p_ix1], - particle_type2=letters[p_ix2], + particle_type1=particle_type1, + particle_type2=particle_type2, force_constant=self.parameters["force_constant"], - interaction_distance=1.5 + 0.5 * (p_ix1 + p_ix2), + interaction_distance=( + ReaddySimpleProcess.radius_for_ix(p_ix1) + + ReaddySimpleProcess.radius_for_ix(p_ix2) - 0.5 + ), ) # add box potential if no periodic boundary @@ -151,10 +163,13 @@ def get_particles_from_simulation(self): "topologies": {}, "particles": {}, } + letters = list(string.ascii_uppercase) for particle in self.readdy_simulation.current_particles: + p_ix = letters.index(particle.type) result["particles"][str(particle.id)] = { "type_name" : particle.type, "position" : particle.pos, + "radius" : ReaddySimpleProcess.radius_for_ix(p_ix), } return result @@ -185,6 +200,7 @@ def random_initial_state(parameters): result["particles"][str(p_ix)] = { "type_name" : letters[p_ix], "position" : random_positions[p_ix], + "radius" : ReaddySimpleProcess.radius_for_ix(p_ix), } return {"monomers" : result} @@ -201,11 +217,15 @@ def run_readdy_simple_process(): composite = readdy_process.generate() initial_state = random_initial_state(readdy_process.parameters) - sim = Engine(composite=composite, initial_state=initial_state) + sim = Engine( + composite=composite, + initial_state=initial_state, + emitter="simularium", + ) - sim.update(10) # ns + sim.update(200000.) # ns - output = sim.emitter.get_timeseries() + output = sim.emitter.get_data() return output diff --git a/vivarium_models/processes/simularium_emitter.py b/vivarium_models/processes/simularium_emitter.py index 88d6246..90fdbb7 100644 --- a/vivarium_models/processes/simularium_emitter.py +++ b/vivarium_models/processes/simularium_emitter.py @@ -10,6 +10,8 @@ AgentData, MetaData, UnitData, + DisplayData, + DISPLAY_TYPE, ) @@ -35,11 +37,12 @@ def emit(self, data: Dict[str, Any]) -> None: key: value for key, value in emit_data.items() if key not in ["time"] } - def get_simularium_fibers(self, time, fibers, actin_radius, trajectory): + def get_simularium_fibers(self, time, fibers, trajectory): """ Shape fiber state data into Simularium fiber agents """ n_agents = 0 + actin_radius = 3.0 time_index = len(trajectory["times"]) trajectory["times"].append(time) trajectory["unique_ids"].append([]) @@ -59,7 +62,7 @@ def get_simularium_fibers(self, time, fibers, actin_radius, trajectory): trajectory["radii"].append(n_agents * [actin_radius]) return trajectory - def get_simularium_monomers(self, time, monomers, actin_radius, trajectory): + def get_simularium_monomers(self, time, monomers, trajectory): """ Shape monomer state data into Simularium agents """ @@ -68,13 +71,21 @@ def get_simularium_monomers(self, time, monomers, actin_radius, trajectory): trajectory["unique_ids"].append([]) trajectory["type_names"].append([]) trajectory["positions"].append([]) + trajectory["radii"].append([]) edge_ids = [] edge_positions = [] for particle_id in monomers["particles"]: particle = monomers["particles"][particle_id] trajectory["unique_ids"][time_index].append(int(particle_id)) trajectory["type_names"][time_index].append(particle["type_name"]) + # HACK needed until simulariumio default display data fix + if particle["type_name"] not in trajectory["display_data"]: + trajectory["display_data"][particle["type_name"]] = DisplayData( + name=particle["type_name"], + display_type=DISPLAY_TYPE.SPHERE, + ) trajectory["positions"][time_index].append(np.array(particle["position"])) + trajectory["radii"][time_index].append(particle["radius"]) # visualize edges between particles for neighbor_id in particle["neighbor_ids"]: neighbor_id_str = str(neighbor_id) @@ -101,7 +112,6 @@ def get_simularium_monomers(self, time, monomers, actin_radius, trajectory): trajectory["unique_ids"][time_index] += [1000 + i for i in range(n_edges)] trajectory["type_names"][time_index] += ["edge" for edge in range(n_edges)] trajectory["positions"][time_index] += n_edges * [[0.0, 0.0, 0.0]] - trajectory["radii"].append(n_agents * [actin_radius]) trajectory["radii"][time_index] += n_edges * [1.0] trajectory["n_subpoints"].append(n_agents * [0]) trajectory["n_subpoints"][time_index] += n_edges * [2] @@ -187,6 +197,7 @@ def get_agent_data_from_jagged_lists(trajectory, scale_factor) -> AgentData: ).to_numpy(dtype=int), subpoints=scale_factor * SimulariumEmitter.get_subpoints_numpy_array(trajectory), + display_data=trajectory["display_data"], ) @staticmethod @@ -212,30 +223,10 @@ def get_simularium_converter( ) ) - @staticmethod - def get_active_simulator(choices) -> str: - """ - Use choices from state to determine which simulator ran - """ - medyan_active = "medyan_active" in choices and choices["medyan_active"] - readdy_active = "readdy_active" in choices and choices["readdy_active"] - cytosim_active = "cytosim_active" in choices and choices["cytosim_active"] - - if medyan_active and not readdy_active and not cytosim_active: - return "medyan" - elif readdy_active and not medyan_active and not cytosim_active: - return "readdy" - elif cytosim_active and not medyan_active and not readdy_active: - return "cytosim" - def get_data(self) -> dict: """ Save the accumulated timeseries history of "emitted" data to file """ - if "readdy_actin" in self.configuration_data: - actin_radius = self.configuration_data["readdy_actin"]["actin_radius"] - else: - actin_radius = 3.0 # TODO add to MEDYAN config box_dimensions = None trajectory = { "times": [], @@ -247,66 +238,19 @@ def get_data(self) -> dict: "radii": [], "n_subpoints": [], "subpoints": [], + "display_data": {}, } times = list(self.saved_data.keys()) times.sort() - vizualize_time_index = 0 for time, state in self.saved_data.items(): - print(f"time = {time}") - index = times.index(time) - prev_simulator = "none" - if index > 0: - prev_simulator = SimulariumEmitter.get_active_simulator( - self.saved_data[times[index - 1]]["choices"] - ) - current_simulator = SimulariumEmitter.get_active_simulator(state["choices"]) - print(f"prev_simulator = {prev_simulator}") - print(f"current_simulator = {current_simulator}") - """visualize the first frame in the simulator that started first - then subsequent frames according the the simulator that ran - right before that time point""" - if prev_simulator == "none": - if current_simulator == "medyan" or current_simulator == "cytosim": - print(f" visualize current {current_simulator}") - if box_dimensions is None: - box_dimensions = np.array(state["fibers_box_extent"]) - trajectory = self.get_simularium_fibers( - vizualize_time_index, - state["fibers"], - actin_radius, - trajectory, - ) - vizualize_time_index += 1 - if current_simulator == "readdy": - print(" visualize current readdy") - trajectory = self.get_simularium_monomers( - vizualize_time_index, - state["monomers"], - actin_radius, - trajectory, - ) - vizualize_time_index += 1 - elif prev_simulator == "medyan" or prev_simulator == "cytosim": - print(f" visualize previous {prev_simulator}") - if box_dimensions is None: - box_dimensions = np.array(state["fibers_box_extent"]) - trajectory = self.get_simularium_fibers( - vizualize_time_index, - state["fibers"], - actin_radius, - trajectory, - ) - vizualize_time_index += 1 - elif prev_simulator == "readdy": - print(" visualize readdy") - trajectory = self.get_simularium_monomers( - vizualize_time_index, - state["monomers"], - actin_radius, - trajectory, - ) - vizualize_time_index += 1 + if box_dimensions is None: + box_dimensions = np.array([state["monomers"]["box_size"]] * 3) + trajectory = self.get_simularium_monomers( + time, + state["monomers"], + trajectory, + ) simularium_converter = SimulariumEmitter.get_simularium_converter( trajectory, box_dimensions, 0.1 ) - simularium_converter.save("out/actin_test") + simularium_converter.save("out/test") From 63bf43934afd295db1596ade2bfcbd3b097466ac Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Tue, 16 May 2023 18:16:45 -0700 Subject: [PATCH 4/5] turn ReaddySimpleProcess into a test that uses Vivarium-ReaDDy's ReaddyProcess instead --- .../processes/readdy_simple_process.py | 251 ------------------ vivarium_models/tests/test_readdy_simple.py | 101 +++++++ 2 files changed, 101 insertions(+), 251 deletions(-) delete mode 100644 vivarium_models/processes/readdy_simple_process.py create mode 100644 vivarium_models/tests/test_readdy_simple.py diff --git a/vivarium_models/processes/readdy_simple_process.py b/vivarium_models/processes/readdy_simple_process.py deleted file mode 100644 index 70cf2eb..0000000 --- a/vivarium_models/processes/readdy_simple_process.py +++ /dev/null @@ -1,251 +0,0 @@ -import os -import string - -from vivarium.core.process import Process -from vivarium.core.directories import PROCESS_OUT_DIR -from vivarium.core.engine import Engine - -import numpy as np -import readdy -from tqdm import tqdm - -from vivarium_models.util import create_monomer_update, monomer_ports_schema - -NAME = "ReaDDy_simple" - - -class ReaddySimpleProcess(Process): - """ - This process uses ReaDDy to model 4 diffusing particles. - """ - - name = NAME - - defaults = { - "n_particles" : 4, - "box_size" : 500., # nm - "temperature_C" : 22.0, - "viscosity" : 8.1, # cP - "periodic_boundary" : False, - "force_constant" : 250., - "internal_timestep" : 1000., # ns - "time_step" : 10000., # ns - } - - - def __init__(self, parameters=None): - super(ReaddySimpleProcess, self).__init__(parameters) - self.create_system() - - - @staticmethod - def radius_for_ix(ix: int): - return 10. + 5. * ix - - - def calculate_diffusionCoefficient(self, radius): - """ - calculates the theoretical diffusion constant of a spherical particle - in a media with viscosity [cP] at temperature [Kelvin]. - with radius radius[nm] - - returns nm^2/s - """ - return ( - (1.38065 * 10 ** (-23) * self.parameters["temperature_K"]) - / (6 * np.pi * self.parameters["viscosity"] * 10 ** (-3) * radius * 10 ** (-9)) - * 10**18 - / 10**9 - ) - - - def create_system(self): - """ - Create a ReaDDy system with the given number of - types of particles freely diffusing in a box. - """ - # create readdy system - box_size_3d = np.array([self.parameters["box_size"]] * 3) - self.readdy_system = readdy.ReactionDiffusionSystem( - box_size=box_size_3d, - periodic_boundary_conditions=[self.parameters["periodic_boundary"]] * 3, - ) - self.parameters["temperature_K"] = self.parameters["temperature_C"] + 273.15 - self.readdy_system.temperature = self.parameters["temperature_K"] - - # add particle species - letters = list(string.ascii_uppercase) - for p_ix in range(self.parameters["n_particles"]): - particle_type = letters[p_ix] - radius = ReaddySimpleProcess.radius_for_ix(p_ix) - diffCoeff = self.calculate_diffusionCoefficient(radius) # nm^2/s - self.readdy_system.add_topology_species(particle_type, diffCoeff) - - # add repulsion potentials for excluded volume - for p_ix1 in range(self.parameters["n_particles"]): - for p_ix2 in range(p_ix1 + 1, self.parameters["n_particles"]): - particle_type1 = letters[p_ix1] - particle_type2 = letters[p_ix2] - self.readdy_system.potentials.add_harmonic_repulsion( - particle_type1=particle_type1, - particle_type2=particle_type2, - force_constant=self.parameters["force_constant"], - interaction_distance=( - ReaddySimpleProcess.radius_for_ix(p_ix1) - + ReaddySimpleProcess.radius_for_ix(p_ix2) - 0.5 - ), - ) - - # add box potential if no periodic boundary - if not self.parameters["periodic_boundary"]: - box_potential_size = box_size_3d - 2.0 - for particle_type in letters[:self.parameters["n_particles"]]: - self.readdy_system.potentials.add_box( - particle_type=particle_type, - force_constant=self.parameters["force_constant"], - origin=-0.5 * box_potential_size, - extent=box_potential_size, - ) - - - def ports_schema(self): - return monomer_ports_schema - - - def create_simulation_with_particles(self, particles): - self.readdy_simulation = self.readdy_system.simulation("CPU") - for p_id in particles: - self.readdy_simulation.add_particle( - type=particles[p_id]["type_name"], - position=particles[p_id]["position"], - ) - - def simulate_readdy(self, timestep): - """ - Simulate in ReaDDy for the given timestep - """ - def loop(): - readdy_actions = self.readdy_simulation._actions - init = readdy_actions.initialize_kernel() - diffuse = readdy_actions.integrator_euler_brownian_dynamics( - self.parameters["internal_timestep"] - ) - calculate_forces = readdy_actions.calculate_forces() - create_nl = readdy_actions.create_neighbor_list( - self.readdy_system.calculate_max_cutoff().magnitude - ) - update_nl = readdy_actions.update_neighbor_list() - react = readdy_actions.reaction_handler_uncontrolled_approximation( - self.parameters["internal_timestep"] - ) - observe = readdy_actions.evaluate_observables() - init() - create_nl() - calculate_forces() - update_nl() - observe(0) - n_steps = int(timestep / self.parameters["internal_timestep"]) - for t in tqdm(range(1, n_steps + 1)): - diffuse() - update_nl() - react() - update_nl() - calculate_forces() - observe(t) - - self.readdy_simulation._run_custom_loop(loop, show_summary=False) - - - def get_particles_from_simulation(self): - result = { - "box_center" : np.zeros(3), - "box_size" : self.parameters["box_size"], - "topologies": {}, - "particles": {}, - } - letters = list(string.ascii_uppercase) - for particle in self.readdy_simulation.current_particles: - p_ix = letters.index(particle.type) - result["particles"][str(particle.id)] = { - "type_name" : particle.type, - "position" : particle.pos, - "radius" : ReaddySimpleProcess.radius_for_ix(p_ix), - } - return result - - - def next_update(self, timestep, states): - print(f"ReaDDy simple process updating by {timestep} ns") - - self.create_simulation_with_particles(states["monomers"]["particles"]) - self.simulate_readdy(timestep) - new_monomers = self.get_particles_from_simulation() - - return create_monomer_update(states["monomers"], new_monomers) - - -def random_initial_state(parameters): - result = { - "box_center" : np.zeros(3), - "box_size" : parameters["box_size"], - "topologies": {}, - "particles": {}, - } - letters = list(string.ascii_uppercase) - random_positions = ( - parameters["box_size"] * - (np.random.uniform(size=(parameters["n_particles"], 3)) - 0.5) - ) - for p_ix in range(parameters["n_particles"]): - result["particles"][str(p_ix)] = { - "type_name" : letters[p_ix], - "position" : random_positions[p_ix], - "radius" : ReaddySimpleProcess.radius_for_ix(p_ix), - } - return {"monomers" : result} - - -# functions to configure and run the process -def run_readdy_simple_process(): - """ - Run a simulation of the process. - - Returns: - The simulation output. - """ - readdy_process = ReaddySimpleProcess() - composite = readdy_process.generate() - initial_state = random_initial_state(readdy_process.parameters) - - sim = Engine( - composite=composite, - initial_state=initial_state, - emitter="simularium", - ) - - sim.update(200000.) # ns - - output = sim.emitter.get_data() - return output - - -def test_readdy_simple_process(): - ''' - Test that the process runs correctly. - This will be executed by pytest. - ''' - output = run_readdy_simple_process() - # TODO: Add assert statements to ensure correct performance. - - -def main(): - '''Simulate the process and plot results.''' - # make an output directory to save results - out_dir = os.path.join(PROCESS_OUT_DIR, NAME) - os.makedirs(out_dir, exist_ok=True) - - output = run_readdy_simple_process() - - -if __name__ == '__main__': - main() diff --git a/vivarium_models/tests/test_readdy_simple.py b/vivarium_models/tests/test_readdy_simple.py new file mode 100644 index 0000000..39270f0 --- /dev/null +++ b/vivarium_models/tests/test_readdy_simple.py @@ -0,0 +1,101 @@ +import string + +import numpy as np +from vivarium.core.engine import Engine + +from vivarium_readdy import ReaddyProcess + + +def radius_for_ix(ix: int): + return 10.0 + 5.0 * ix + + +def random_free_particles_state(n_particles, box_size, time_units, spatial_units): + box_size_3d = np.array([box_size] * 3) + particles = {} + letters = list(string.ascii_uppercase) + random_positions = box_size * (np.random.uniform(size=(n_particles, 3)) - 0.5) + for p_ix in range(n_particles): + particles[str(p_ix)] = { + "type_name": letters[p_ix], + "position": random_positions[p_ix], + "radius": radius_for_ix(p_ix), + } + return { + "time_units": time_units, + "spatial_units": spatial_units, + "monomers": { + "box_size": box_size_3d, + "topologies": {}, + "particles": particles, + }, + } + + +def run_readdy_simple_process(): + """ + Run a simulation of particles freely diffusing in ReaDDy. + + Returns: + The simulation output. + """ + n_particles = 4 + time_units = "ns" + spatial_units = "nm" + box_size = 500.0 + letters = list(string.ascii_uppercase) + particle_radii = {} + for p_ix in range(n_particles): + particle_radii[letters[p_ix]] = radius_for_ix(p_ix) + readdy_process = ReaddyProcess( + { + "time_step": 10000, + "internal_timestep": 1000, + "box_size": box_size, + "periodic_boundary": False, + "temperature_C": 22.0, + "viscosity": 8.1, # cP, cytoplasm + "force_constant": 250.0, + "n_cpu": 4, + "particle_radii": particle_radii, + "topology_particles": [], + "reactions": [], + "time_units": time_units, + "spatial_units": spatial_units, + } + ) + composite = readdy_process.generate() + engine = Engine( + composite=composite, + initial_state=random_free_particles_state( + n_particles, box_size, time_units, spatial_units + ), + emitter="simularium", + ) + engine.update(200000) + return engine.emitter.get_data() + + +def test_readdy_simple_process(): + """ + Test that the process runs correctly. + This will be executed by pytest. + """ + output = run_readdy_simple_process()["simularium_converter"]._data + np.testing.assert_allclose(output.meta_data.box_size, np.array([50.0, 50.0, 50.0])) + assert output.agent_data.display_data["C"].name == "C" + np.testing.assert_allclose( + output.agent_data.times, np.arange(0.0, 200001.0, 10000.0) + ) + np.testing.assert_allclose(np.unique(output.agent_data.n_agents), np.array([4])) + assert list(np.unique(output.agent_data.types)) == ["A", "B", "C", "D"] + np.testing.assert_allclose( + np.unique(output.agent_data.radii), np.array([1.0, 1.5, 2.0, 2.5]) + ) + assert output.time_units.name == "ns" + assert output.spatial_units.name == "nm" + assert output.spatial_units.magnitude == 10.0 + + +if __name__ == "__main__": + test_readdy_simple_process() From ac65c767da2a6dd03b2f71f4b7bb49128fa4ab79 Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Tue, 16 May 2023 18:17:52 -0700 Subject: [PATCH 5/5] clean up and lint --- setup.py | 3 +- vivarium_models/__init__.py | 3 +- vivarium_models/composites/actin_fiber.py | 6 +- .../composites/filament_alternatives.py | 6 +- vivarium_models/composites/moving_anchor.py | 40 ++--- vivarium_models/data/fibers.py | 22 +-- vivarium_models/processes/fiber_to_monomer.py | 10 +- vivarium_models/processes/monomer_to_fiber.py | 10 +- .../processes/readdy_actin_process.py | 17 +- .../processes/simularium_emitter.py | 40 ++--- .../processes/visualize_filament.py | 126 -------------- .../processes/visualize_monomer.py | 162 ------------------ vivarium_models/tests/test_readdy_actin.py | 39 ----- vivarium_models/util.py | 84 --------- 14 files changed, 75 insertions(+), 493 deletions(-) delete mode 100644 vivarium_models/processes/visualize_filament.py delete mode 100644 vivarium_models/processes/visualize_monomer.py delete mode 100644 vivarium_models/tests/test_readdy_actin.py diff --git a/setup.py b/setup.py index c4be0b4..d50ebb8 100644 --- a/setup.py +++ b/setup.py @@ -41,8 +41,9 @@ "vivarium-core", "vivarium_medyan @ git+https://github.com/vivarium-collective/vivarium-MEDYAN.git", "vivarium_cytosim @ git+https://github.com/vivarium-collective/vivarium-cytosim.git", + "vivarium_readdy @ git+https://github.com/vivarium-collective/vivarium-ReaDDy.git@monomers-ports", "simularium_readdy_models @ git+https://github.com/simularium/readdy-models.git", - "simulariumio>=1.5.0", + "simulariumio", ] extra_requirements = { diff --git a/vivarium_models/__init__.py b/vivarium_models/__init__.py index 44b0426..d3dbf16 100644 --- a/vivarium_models/__init__.py +++ b/vivarium_models/__init__.py @@ -14,7 +14,8 @@ def get_module_version(): from .processes import ReaddyActinProcess # noqa: F401 -from vivarium.core.registry import emitter_registry # noqa: F401 from .processes.simularium_emitter import SimulariumEmitter # noqa: F401 +from vivarium.core.registry import emitter_registry + emitter_registry.register("simularium", SimulariumEmitter) diff --git a/vivarium_models/composites/actin_fiber.py b/vivarium_models/composites/actin_fiber.py index e10e89b..7e96789 100644 --- a/vivarium_models/composites/actin_fiber.py +++ b/vivarium_models/composites/actin_fiber.py @@ -81,7 +81,7 @@ def generate_topology(self, config): } -def test_actin_fiber(): +def run_actin_fiber(): initial_state = centered_initial_fibers() initial_state["choices"] = {"medyan_active": True, "readdy_active": False} medyan_config = { @@ -100,8 +100,8 @@ def test_actin_fiber(): emit_processes=True, ) engine.update(5) - engine.emitter.get_data() + return engine.emitter.get_data() if __name__ == "__main__": - test_actin_fiber() + run_actin_fiber() diff --git a/vivarium_models/composites/filament_alternatives.py b/vivarium_models/composites/filament_alternatives.py index ccfc1ed..92cfdf8 100644 --- a/vivarium_models/composites/filament_alternatives.py +++ b/vivarium_models/composites/filament_alternatives.py @@ -61,7 +61,7 @@ def generate_topology(self, config): } -def test_filament_alternatives(): +def run_filament_alternatives(): initial_state = centered_initial_fibers() initial_state["choices"] = {"medyan_active": False, "cytosim_active": True} medyan_config = { @@ -87,8 +87,8 @@ def test_filament_alternatives(): emit_processes=True, ) engine.update(6) - engine.emitter.get_data() + return engine.emitter.get_data() if __name__ == "__main__": - test_filament_alternatives() + run_filament_alternatives() diff --git a/vivarium_models/composites/moving_anchor.py b/vivarium_models/composites/moving_anchor.py index 998bc6e..7cfa2dc 100644 --- a/vivarium_models/composites/moving_anchor.py +++ b/vivarium_models/composites/moving_anchor.py @@ -48,13 +48,8 @@ def next_update(self, timestep, state): anchor_moved = anchor + np.array(self.parameters["movement_vector"]) points[self.parameters["point_index"]] = anchor_moved - return { - "fibers": { - self.parameters["fiber_id"]: { - "points": points - } - } - } + return {"fibers": {self.parameters["fiber_id"]: {"points": points}}} + class AnchorMoverCytosim(CytosimProcess): new_defaults = { @@ -76,13 +71,12 @@ def next_update(self, timestep, state): update = super().next_update(timestep, state) return update + class BucklingSqueeze(Composer): defaults = { - "periodic_event": { - "periods": [1.0] - }, + "periodic_event": {"periods": [1.0]}, "cytosim_anchor_mover_squeeze": { - 'model_name': 'buckling_squeeze', + "model_name": "buckling_squeeze", "movement_vector": [5, 0, 0], "fiber_id": "1", }, @@ -102,13 +96,15 @@ def generate_processes(self, config): periodic_event = PeriodicEvent(config["periodic_event"]) # cytosim_squeeze = CytosimProcess(config['cytosim_squeeze']) # anchor_mover = AnchorMover(config["anchor_mover"]) - cytosim_anchor_mover_squeeze = AnchorMoverCytosim(config["cytosim_anchor_mover_squeeze"]) + cytosim_anchor_mover_squeeze = AnchorMoverCytosim( + config["cytosim_anchor_mover_squeeze"] + ) return { "periodic_event": periodic_event, # 'cytosim_squeeze': cytosim_squeeze, # "anchor_mover": anchor_mover, - "cytosim_anchor_mover_squeeze": cytosim_anchor_mover_squeeze + "cytosim_anchor_mover_squeeze": cytosim_anchor_mover_squeeze, } def generate_topology(self, config): @@ -131,16 +127,17 @@ def generate_topology(self, config): } -def test_buckling_squeeze(): +def run_buckling_squeeze(): initial_state = single_fiber() - initial_state['choices'] = 'N/A' + initial_state["choices"] = "N/A" cytosim_config = { - 'actin_segmentation': 0.005, + "actin_segmentation": 0.005, "timestep": 5, - "template_directory": "vivarium_models/templates/"} + "template_directory": "vivarium_models/templates/", + } buckling_squeeze_config = { # "cytosim_squeeze": cytosim_config - "cytosim_anchor_mover_squeeze": cytosim_config + "cytosim_anchor_mover_squeeze": cytosim_config } buckling_squeeze = BucklingSqueeze(buckling_squeeze_config) composite = buckling_squeeze.generate() @@ -153,11 +150,8 @@ def test_buckling_squeeze(): emit_processes=True, ) engine.update(100) - - output = engine.emitter.get_data() - # import ipdb; ipdb.set_trace() - + return engine.emitter.get_data() if __name__ == "__main__": - test_buckling_squeeze() + run_buckling_squeeze() diff --git a/vivarium_models/data/fibers.py b/vivarium_models/data/fibers.py index 43a814f..1ed5b72 100644 --- a/vivarium_models/data/fibers.py +++ b/vivarium_models/data/fibers.py @@ -228,13 +228,14 @@ def centered_initial_fibers(): ) return result + def map_fibers(f, initial_fibers): result = copy.deepcopy(initial_fibers) for fiber_id in initial_fibers["fibers"]: fiber_points = initial_fibers["fibers"][fiber_id]["points"] for point_index in range(len(fiber_points)): - result["fibers"][fiber_id]["points"][point_index] = ( - f(fiber_points[point_index], initial_fibers["fibers_box_extent"]) + result["fibers"][fiber_id]["points"][point_index] = f( + fiber_points[point_index], initial_fibers["fibers_box_extent"] ) return result @@ -243,7 +244,7 @@ def single_fiber(): fiber = { "fibers_box_extent": np.array([4000.0, 2000.0, 2000.0]), "fibers": { - '1': { + "1": { "type_name": "Actin-Polymer", "points": [ np.array([-0.250, 0.000, 0.000]), @@ -297,14 +298,15 @@ def single_fiber(): np.array([0.230, 0.000, 0.000]), np.array([0.240, 0.000, 0.000]), np.array([0.250, 0.000, 0.000]), - ]}}} + ], + } + }, + } - scaled = map_fibers( - lambda fiber_points, box_extent: fiber_points * 1000, - fiber) + scaled = map_fibers(lambda fiber_points, box_extent: fiber_points * 1000, fiber) return scaled -if __name__ == '__main__': - fiber = single_fiber() - import ipdb; ipdb.set_trace() + +if __name__ == "__main__": + single_fiber() diff --git a/vivarium_models/processes/fiber_to_monomer.py b/vivarium_models/processes/fiber_to_monomer.py index 14e7acd..72737f3 100644 --- a/vivarium_models/processes/fiber_to_monomer.py +++ b/vivarium_models/processes/fiber_to_monomer.py @@ -1,7 +1,7 @@ import numpy as np from vivarium.core.process import Deriver -from vivarium.core.engine import Engine, pf +from vivarium.core.engine import Engine from simularium_readdy_models.actin import ActinGenerator, ActinTestData, FiberData from ..util import create_monomer_update @@ -114,7 +114,7 @@ def get_initial_fiber_data(): return fibers_dict -def test_fiber_to_monomer(): +def run_fiber_to_monomer(): fiber_data = get_initial_fiber_data() fiber_to_monomer = FiberToMonomer() @@ -132,10 +132,8 @@ def test_fiber_to_monomer(): ) engine.update(1.0) - - output = engine.emitter.get_data() - print(pf(output)) + return engine.emitter.get_data() if __name__ == "__main__": - test_fiber_to_monomer() + run_fiber_to_monomer() diff --git a/vivarium_models/processes/monomer_to_fiber.py b/vivarium_models/processes/monomer_to_fiber.py index 7f09478..84c6c25 100644 --- a/vivarium_models/processes/monomer_to_fiber.py +++ b/vivarium_models/processes/monomer_to_fiber.py @@ -1,7 +1,7 @@ import numpy as np from vivarium.core.process import Deriver -from vivarium.core.engine import Engine, pf +from vivarium.core.engine import Engine from simularium_readdy_models.actin import ActinUtil, ActinTestData from ..util import agents_update @@ -184,7 +184,7 @@ def generate_fibers_from_monomers(monomers, box_size=500.0): return result -def test_fiber_to_monomer(): +def run_fiber_to_monomer(): monomer_data = ActinTestData.linear_actin_monomers() monomer_to_fiber = MonomerToFiber() @@ -202,10 +202,8 @@ def test_fiber_to_monomer(): ) engine.update(1.0) - - output = engine.emitter.get_data() - print(pf(output)) + return engine.emitter.get_data() if __name__ == "__main__": - test_fiber_to_monomer() + run_fiber_to_monomer() diff --git a/vivarium_models/processes/readdy_actin_process.py b/vivarium_models/processes/readdy_actin_process.py index a7ceadd..d638fcc 100644 --- a/vivarium_models/processes/readdy_actin_process.py +++ b/vivarium_models/processes/readdy_actin_process.py @@ -1,7 +1,7 @@ import numpy as np from vivarium.core.process import Process -from vivarium.core.engine import Engine, pf +from vivarium.core.engine import Engine from vivarium.core.control import run_library_cli from tqdm import tqdm @@ -12,7 +12,8 @@ ActinAnalyzer, ) from simularium_readdy_models import ReaddyUtil -from vivarium_models.util import create_monomer_update, format_monomer_results, monomer_ports_schema +from vivarium_readdy import create_monomer_update, monomer_ports_schema +from ..util import format_monomer_results from vivarium_models.library.scan import Scan NAME = "ReaDDy_actin" @@ -162,7 +163,7 @@ def get_monomer_data(): return {"monomers": monomer_data} -def test_readdy_actin_process(): +def run_readdy_actin_process(): monomer_data = get_monomer_data() readdy_actin_process = ReaddyActinProcess() @@ -179,12 +180,10 @@ def test_readdy_actin_process(): ) engine.update(0.0000001) # 1e3 steps + return engine.emitter.get_data() - output = engine.emitter.get_data() - print(pf(output)) - -def test_scan_readdy(): +def run_scan_readdy(): monomer_data = get_monomer_data() parameters = { @@ -308,10 +307,10 @@ def filament_straightness(results): scan = Scan(parameters, ReaddyActinProcess, 0.0000001, metrics=metrics) results = scan.run_scan() - print(results) + return results -library = {"0": test_readdy_actin_process, "1": test_scan_readdy} +library = {"0": run_readdy_actin_process, "1": run_scan_readdy} if __name__ == "__main__": diff --git a/vivarium_models/processes/simularium_emitter.py b/vivarium_models/processes/simularium_emitter.py index 90fdbb7..1cc3196 100644 --- a/vivarium_models/processes/simularium_emitter.py +++ b/vivarium_models/processes/simularium_emitter.py @@ -1,4 +1,4 @@ -from typing import Any, Dict +from typing import Any, Dict, List from vivarium.core.emitter import Emitter @@ -19,7 +19,7 @@ class SimulariumEmitter(Emitter): def __init__(self, config: Dict[str, str]) -> None: super().__init__(config) self.configuration_data = None - self.saved_data: Dict[float, Dict[str, Any]] = {} + self.saved_data: List[Dict[str, Any]] = [] def emit(self, data: Dict[str, Any]) -> None: """ @@ -31,11 +31,7 @@ def emit(self, data: Dict[str, Any]) -> None: self.configuration_data = data["data"] assert "processes" in self.configuration_data, "please emit processes" if data["table"] == "history": - emit_data = data["data"] - time = emit_data["time"] - self.saved_data[time] = { - key: value for key, value in emit_data.items() if key not in ["time"] - } + self.saved_data.append(data["data"]) def get_simularium_fibers(self, time, fibers, trajectory): """ @@ -177,7 +173,7 @@ def get_agent_data_from_jagged_lists(trajectory, scale_factor) -> AgentData: Shape a dictionary of jagged lists into a Simularium AgentData object """ return AgentData( - times=np.arange(len(trajectory["times"])), + times=trajectory["times"], n_agents=np.array(trajectory["n_agents"]), viz_types=SimulariumEmitter.fill_df( pd.DataFrame(trajectory["viz_types"]), 1000.0 @@ -202,13 +198,12 @@ def get_agent_data_from_jagged_lists(trajectory, scale_factor) -> AgentData: @staticmethod def get_simularium_converter( - trajectory, box_dimensions, scale_factor + trajectory, box_dimensions, scale_factor, time_units, spatial_units ) -> TrajectoryConverter: """ Shape a dictionary of jagged lists into a Simularium TrajectoryData object and provide it to a TrajectoryConverter for conversion """ - spatial_units = UnitData("nm") spatial_units.multiply(1 / scale_factor) return TrajectoryConverter( TrajectoryData( @@ -218,7 +213,7 @@ def get_simularium_converter( agent_data=SimulariumEmitter.get_agent_data_from_jagged_lists( trajectory, scale_factor ), - time_units=UnitData("count"), + time_units=time_units, spatial_units=spatial_units, ) ) @@ -227,7 +222,9 @@ def get_data(self) -> dict: """ Save the accumulated timeseries history of "emitted" data to file """ - box_dimensions = None + box_size = None + time_units = None + spatial_units = None trajectory = { "times": [], "n_agents": [], @@ -240,17 +237,20 @@ def get_data(self) -> dict: "subpoints": [], "display_data": {}, } - times = list(self.saved_data.keys()) - times.sort() - for time, state in self.saved_data.items(): - if box_dimensions is None: - box_dimensions = np.array([state["monomers"]["box_size"]] * 3) + for state in self.saved_data: + if box_size is None: + box_size = state["monomers"]["box_size"] + if time_units is None: + time_units = state["time_units"] + if spatial_units is None: + spatial_units = state["spatial_units"] trajectory = self.get_simularium_monomers( - time, + state["time"], state["monomers"], trajectory, ) simularium_converter = SimulariumEmitter.get_simularium_converter( - trajectory, box_dimensions, 0.1 + trajectory, box_size, 0.1, UnitData(time_units), UnitData(spatial_units) ) - simularium_converter.save("out/test") + simularium_converter.save("out/test") + return {"simularium_converter": simularium_converter} diff --git a/vivarium_models/processes/visualize_filament.py b/vivarium_models/processes/visualize_filament.py deleted file mode 100644 index 07765ce..0000000 --- a/vivarium_models/processes/visualize_filament.py +++ /dev/null @@ -1,126 +0,0 @@ -import numpy as np - -from vivarium.core.process import Deriver -from vivarium.core.engine import Engine, pf -from simularium_readdy_models.actin import ActinTestData -from simulariumio import ( - TrajectoryConverter, - TrajectoryData, - AgentData, - MetaData, - UnitData, -) - - -class VisualizeFilament(Deriver): - defaults = {} - - def __init__(self, parameters=None): - super().__init__(parameters) - - def ports_schema(self): - return { - "filaments": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "points": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - "simularium_json": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - } - - def next_update(self, timestep, states): - filaments = states["filaments"] - - box_size = 150.0 - actin_radius = 3.0 - n_agents = 0 - unique_ids = [] - type_names = [] - n_subpoints = [] - subpoints = [] - for filament_id in filaments: - filament = filaments[filament_id] - unique_ids.append(filament_id) - type_names.append(filament["type_name"]) - n_subpoints.append(len(filament["points"])) - subpoints.append(filament["points"]) - n_agents += 1 - subpoints = np.array([subpoints], dtype=float) - - # HACK around a Simularium Viewer bug - subpoints[0][0][1][2] += 0.001 - - simularium_json = TrajectoryConverter( - TrajectoryData( - meta_data=MetaData( - box_size=np.array([box_size, box_size, box_size]), - ), - agent_data=AgentData( - times=np.array([0]), - n_agents=np.array([n_agents]), - viz_types=1001.0 * np.ones((1, n_agents)), - unique_ids=np.array([unique_ids]), - types=[type_names], - positions=np.zeros((1, n_agents, 3)), - radii=actin_radius * np.ones((1, n_agents)), - n_subpoints=np.array([n_subpoints]), - subpoints=subpoints, - ), - time_units=UnitData("ns"), # nanoseconds - spatial_units=UnitData("nm"), # nanometers - ) - ).to_JSON() - - return {"simularium_json": simularium_json} - - -def get_initial_filament_data(): - fibers = ActinTestData.linear_actin_fiber() - filaments_dict = {} - for fiber in fibers: - filaments_dict[fiber.fiber_id] = dict(fiber) - return filaments_dict - - -def test_visualize_filament(): - filament_data = get_initial_filament_data() - visualize_filament = VisualizeFilament() - - engine = Engine( - { - "processes": {"visualize_filament": visualize_filament}, - "topology": { - "visualize_filament": { - "filaments": ("filaments",), - "simularium_json": ("simularium_json",), - } - }, - "initial_state": {"filaments": filament_data, "simularium_json": ""}, - } - ) - - engine.update(1.0) - - output = engine.emitter.get_data() - print(pf(output)) - simularium_json = output[0.0]["simularium_json"] - with open("filaments.simularium", "w+") as outfile: - outfile.write(simularium_json) - print("Saved Simularium output to filaments.simularium") - - -if __name__ == "__main__": - test_visualize_filament() diff --git a/vivarium_models/processes/visualize_monomer.py b/vivarium_models/processes/visualize_monomer.py deleted file mode 100644 index 08f1de9..0000000 --- a/vivarium_models/processes/visualize_monomer.py +++ /dev/null @@ -1,162 +0,0 @@ -import numpy as np - -from vivarium.core.process import Deriver -from vivarium.core.engine import Engine, pf - -from simularium_readdy_models.actin import ActinTestData -from simulariumio import ( - TrajectoryConverter, - TrajectoryData, - AgentData, - MetaData, - UnitData, -) - - -class VisualizeMonomer(Deriver): - defaults = {} - - def __init__(self, parameters=None): - super().__init__(parameters) - - def ports_schema(self): - return { - "monomers": { - "topologies": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "particle_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - "particles": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "position": { - "_default": np.zeros(3), - "_updater": "set", - "_emit": True, - }, - "neighbor_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - }, - "simularium_json": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - } - - def next_update(self, timestep, states): - monomers = states["monomers"] - - box_size = 150.0 - actin_radius = 3.0 - unique_ids = [] - type_names = [] - positions = [] - edge_ids = [] - edge_positions = [] - for particle_id in monomers["particles"]: - particle = monomers["particles"][particle_id] - unique_ids.append(particle_id) - type_names.append(particle["type_name"]) - positions.append(particle["position"]) - for neighbor_id in particle["neighbor_ids"]: - edge = (particle_id, neighbor_id) - reverse_edge = ( - neighbor_id, - particle_id, - ) - if edge not in edge_ids and reverse_edge not in edge_ids: - edge_ids.append(edge) - edge_positions.append( - [ - particle["position"], - monomers["particles"][neighbor_id]["position"], - ] - ) - # add fiber agents for edges - n_agents = len(unique_ids) - n_edges = len(edge_ids) - viz_types = 1000.0 * np.ones((1, n_agents + n_edges)) - viz_types[:, n_agents:] += 1 - unique_ids = np.array([unique_ids + [1000 + i for i, e in enumerate(edge_ids)]]) - type_names += ["edge" for edge in edge_ids] - positions += n_edges * [np.zeros(3)] - radii = actin_radius * np.ones((1, n_agents + n_edges)) - radii[:, n_agents:] = 1.0 - n_subpoints = np.zeros(n_agents + n_edges) - n_subpoints[n_agents:] += 2 - subpoints = np.zeros((n_agents + n_edges, 2, 3)) - subpoints[n_agents:] = np.array(edge_positions) - - simularium_json = TrajectoryConverter( - TrajectoryData( - meta_data=MetaData( - box_size=np.array([box_size, box_size, box_size]), - ), - agent_data=AgentData( - times=np.array([0]), - n_agents=np.array([n_agents + n_edges]), - viz_types=viz_types, - unique_ids=unique_ids, - types=[type_names], - positions=np.array([positions]), - radii=radii, - n_subpoints=np.array([n_subpoints]), - subpoints=np.array([subpoints]), - ), - time_units=UnitData("ns"), # nanoseconds - spatial_units=UnitData("nm"), # nanometers - ) - ).to_JSON() - - return {"simularium_json": simularium_json} - - -def test_visualize_monomer(): - monomer_data = ActinTestData.linear_actin_monomers() - visualize_monomer = VisualizeMonomer() - - engine = Engine( - { - "processes": {"visualize_monomer": visualize_monomer}, - "topology": { - "visualize_monomer": { - "monomers": ("monomers",), - "simularium_json": ("simularium_json",), - } - }, - "initial_state": {"monomers": monomer_data, "simularium_json": ""}, - } - ) - - engine.update(1.0) - - output = engine.emitter.get_data() - print(pf(output)) - simularium_json = output[0.0]["simularium_json"] - with open("monomers.simularium", "w+") as outfile: - outfile.write(simularium_json) - print("Saved Simularium output to monomers.simularium") - - -if __name__ == "__main__": - test_visualize_monomer() diff --git a/vivarium_models/tests/test_readdy_actin.py b/vivarium_models/tests/test_readdy_actin.py deleted file mode 100644 index b1ccb1d..0000000 --- a/vivarium_models/tests/test_readdy_actin.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -Tests for actin ReaDDy models -""" - -# from vivarium_models import ReaddyActinProcess - - -def test_readdy_actin_process(): - """ - Test the initial ReaDDy actin process - """ - assert True - # output = ReaddyActinProcess.run_readdy_actin_process()[5e-09]["monomers"] - # found_monomer = False - # found_dimer = False - # assert len(output["topologies"]) == 2 - # for t in output["topologies"]: - # top = output["topologies"][t] - # if top["type_name"] == "Actin-Monomer": - # found_monomer = True - # particle_ids = top["particle_ids"] - # assert len(particle_ids) == 1 - # assert ( - # "actin#free" in output["particles"][str(particle_ids[0])]["type_name"] - # ) - # assert len(output["particles"][str(particle_ids[0])]["neighbor_ids"]) == 0 - # if top["type_name"] == "Arp23-Dimer": - # found_dimer = True - # particle_ids = top["particle_ids"] - # assert len(particle_ids) == 2 - # assert "arp2" in output["particles"][str(particle_ids[0])]["type_name"] - # assert output["particles"][str(particle_ids[0])]["neighbor_ids"] == [1] - # assert "arp3" in output["particles"][str(particle_ids[1])]["type_name"] - # assert output["particles"][str(particle_ids[1])]["neighbor_ids"] == [0] - # assert found_monomer - # assert found_dimer diff --git a/vivarium_models/util.py b/vivarium_models/util.py index c928614..e3d97b2 100644 --- a/vivarium_models/util.py +++ b/vivarium_models/util.py @@ -1,90 +1,6 @@ import numpy as np -monomer_ports_schema = { - "monomers": { - "box_center": { - "_default": np.array([1000.0, 0.0, 0.0]), - "_updater": "set", - "_emit": True, - }, - "box_size": { - "_default": 500.0, - "_updater": "set", - "_emit": True, - }, - "topologies": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "particle_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - } - }, - "particles": { - "*": { - "type_name": { - "_default": "", - "_updater": "set", - "_emit": True, - }, - "position": { - "_default": np.zeros(3), - "_updater": "set", - "_emit": True, - }, - "neighbor_ids": { - "_default": [], - "_updater": "set", - "_emit": True, - }, - "radius": { - "_default": 1.0, - "_updater": "set", - "_emit": True, - }, - } - }, - } -} - - -def agents_update(existing, projected): - update = {"_add": [], "_delete": []} - - for id, state in projected.items(): - if id in existing: - update[id] = state - else: - update["_add"].append({"key": id, "state": state}) - - for existing_id in existing.keys(): - if existing_id not in projected: - update["_delete"].append(existing_id) - - return update - - -def create_monomer_update(previous_monomers, new_monomers): - topologies_update = agents_update( - previous_monomers["topologies"], new_monomers["topologies"] - ) - - particles_update = agents_update( - previous_monomers["particles"], new_monomers["particles"] - ) - - return { - "monomers": {"topologies": topologies_update, "particles": particles_update} - } - - def format_monomer_results(results): """ Workaround since numpy arrays are not preserved in Vivarium?