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/setup.py b/setup.py index 6eee88f..d50ebb8 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,7 @@ "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", ] diff --git a/vivarium_models/processes/readdy_actin_process.py b/vivarium_models/processes/readdy_actin_process.py index d1d1733..d638fcc 100644 --- a/vivarium_models/processes/readdy_actin_process.py +++ b/vivarium_models/processes/readdy_actin_process.py @@ -12,7 +12,8 @@ ActinAnalyzer, ) from simularium_readdy_models import ReaddyUtil -from vivarium_models.util import create_monomer_update, format_monomer_results +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" @@ -69,53 +70,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/simularium_emitter.py b/vivarium_models/processes/simularium_emitter.py index 88d6246..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 @@ -10,6 +10,8 @@ AgentData, MetaData, UnitData, + DisplayData, + DISPLAY_TYPE, ) @@ -17,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: """ @@ -29,17 +31,14 @@ 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, 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 +58,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 +67,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 +108,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] @@ -167,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 @@ -187,17 +193,17 @@ 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 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( @@ -207,36 +213,18 @@ 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, ) ) - @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 + box_size = None + time_units = None + spatial_units = None trajectory = { "times": [], "n_agents": [], @@ -247,66 +235,22 @@ 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 + 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( + 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/actin_test") + simularium_converter.save("out/test") + return {"simularium_converter": simularium_converter} 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() diff --git a/vivarium_models/util.py b/vivarium_models/util.py index 9be67e5..e3d97b2 100644 --- a/vivarium_models/util.py +++ b/vivarium_models/util.py @@ -1,36 +1,6 @@ import numpy as np -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?