Merge pull request #17 from virtualcell/image-fielddata

add image-based field data support

Author: jcschaff

examples/scripts/fielddata_from_image_workflow.py

@@ -0,0 +1,60 @@
+import os
+import tempfile
+from pathlib import Path
+
+import numpy as np
+
+from pyvcell.sim_results.var_types import NDArray3D
+from pyvcell.vcml import VcmlReader
+from pyvcell.vcml.field import Field
+from pyvcell.vcml.vcml_simulation import VcmlSpatialSimulation as Solver
+
+
+def create_sinusoid(shape: tuple[int, ...], freq: float) -> NDArray3D:
+    indices = np.indices(shape)
+    sinusoid: NDArray3D = (
+        np.cos(freq * indices[0, :, :, :]) * np.sin(freq * indices[1, :, :, :]) * np.sin(freq * indices[2, :, :, :])
+    )
+    return sinusoid.astype(dtype=np.float64)
+
+
+with tempfile.TemporaryDirectory() as temp_dir_name, Path(temp_dir_name) as temp_dir:
+    print(f"temp_dir: {temp_dir}, exists={temp_dir.exists()}")
+    # ----- make a workspace
+    workspace_dir = temp_dir / "workspace"
+    workspace_dir.mkdir(parents=True, exist_ok=True)
+    print(f"workspace_dir: {workspace_dir}, exists={workspace_dir.exists()}")
+    sim_dir = workspace_dir / "sim1_dir"
+    sim_dir.mkdir(parents=True, exist_ok=True)
+
+    print(f"sim_dir: {sim_dir}, exists={sim_dir.exists()}")
+
+    # ---- read in VCML file
+    model_path = Path(os.getcwd()).parent / "models" / "SmallSpatialProject_3D.vcml"
+    bio_model = VcmlReader.biomodel_from_file(vcml_path=model_path)
+
+    # ---- get the species mappings for species "s0" and "s1"
+    app = bio_model.applications[0]
+    s1_mapping = next(s for s in app.species_mappings if s.species_name == "s1")
+    s0_mapping = next(s for s in app.species_mappings if s.species_name == "s0")
+
+    # ---- set initial concentrations to reference external field data
+    s0_mapping.init_conc = "vcField('test2_lsm_DEMO', 'species0_cyt', 0.5, 'Volume')"
+    s1_mapping.init_conc = "vcField('checkerboard', 'v', 0.0, 'Volume')"
+
+    sim = app.add_sim(name="new_sim", duration=10.0, output_time_step=0.1, mesh_size=(20, 20, 20))
+    fields = Field.create_fields(bio_model=bio_model, sim=sim)
+    print(fields)
+
+    shape = fields[0].data_nD.shape
+    fields[0].data_nD = np.multiply(create_sinusoid(shape=shape, freq=0.5), 8.0)
+    fields[1].data_nD = np.multiply(create_sinusoid(shape=shape, freq=0.3), 4.0)
+
+    # ---- add field data to the simulation
+
+    sim1_result = Solver(bio_model=bio_model, out_dir=sim_dir, fields=fields).run(sim.name)
+    print([c.label for c in sim1_result.channel_data])
+    print(sim1_result.time_points[::11])
+    sim1_result.plotter.plot_slice_3d(time_index=0, channel_id="s0")
+    sim1_result.plotter.plot_slice_3d(time_index=0, channel_id="s1")
+    sim1_result.plotter.plot_concentrations()

examples/scripts/fielddata_workflow.py …s/scripts/fielddata_from_sim_workflow.pyexamples/scripts/fielddata_workflow.py renamed to examples/scripts/fielddata_from_sim_workflow.py examples/scripts/fielddata_workflow.py renamed to examples/scripts/fielddata_from_sim_workflow.py

@@ -3,7 +3,7 @@
 import tempfile
 from pathlib import Path
 
-from pyvcell.vcml import Simulation, VcmlReader
+from pyvcell.vcml import VcmlReader
 from pyvcell.vcml.vcml_simulation import VcmlSpatialSimulation as Solver
 
 with tempfile.TemporaryDirectory() as temp_dir_name, Path(temp_dir_name) as temp_dir:
@@ -26,16 +26,15 @@
     s0_mapping = next(s for s in app.species_mappings if s.species_name == "s0")
 
     # ---- add a simulation to the first application in the biomodel (didn't already have a simulation in the VCML file)
-    new_sim = Simulation(name="new_sim", duration=10.0, output_time_step=0.1, mesh_size=(20, 20, 20))
-    app.simulations.append(new_sim)
+    sim = app.add_sim(name="new_sim", duration=10.0, output_time_step=0.1, mesh_size=(20, 20, 20))
 
     # ---- set the initial concentration of species "s0" and "s1" in the first application
     s0_mapping.init_conc = "3+sin(x)+cos(y)+sin(z)"
     s1_mapping.init_conc = "3+sin(x+y+z)"
 
     # ---- run simulation, store in sim1_dir, and plot results
     # >>>>> This forms the data for the "Field Data" identified by 'sim1_dir' <<<<<<
-    sim1_result = Solver(bio_model=bio_model1, out_dir=sim1_dir).run(new_sim.name)
+    sim1_result = Solver(bio_model=bio_model1, out_dir=sim1_dir).run(sim.name)
     print([c.label for c in sim1_result.channel_data])
     print(sim1_result.time_points[::11])
     sim1_result.plotter.plot_slice_3d(time_index=0, channel_id="s0")
@@ -47,7 +46,7 @@
     s1_mapping.init_conc = "5.0"
     # ---- re-run simulation and store in sim2_dir
     # note that the solution of s0 draws from the data from sim1_dir
-    sim2_result = Solver(bio_model=bio_model1, out_dir=sim2_dir).run(new_sim.name)
+    sim2_result = Solver(bio_model=bio_model1, out_dir=sim2_dir).run(sim.name)
     sim2_result.plotter.plot_slice_3d(time_index=0, channel_id="s0")
     sim2_result.plotter.plot_slice_3d(time_index=0, channel_id="s1")
     sim2_result.plotter.plot_concentrations()

poetry.lock

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pyvcell"
-version = "0.1.7"
+version = "0.1.8"
 description = "This is the python wrapper for vcell modeling and simulation"
 authors = ["Jim Schaff <[email protected]>"]
 repository = "https://github.com/virtualcell/pyvcell"
@@ -32,12 +32,13 @@ matplotlib = "^3.10.0"
 lxml = "^5.3.1"
 imageio = "^2.37.0"
 tensorstore = "^0.1.72"
-libvcell = "^0.0.8"
+libvcell = "^0.0.12"
 trame = "^3.8.1"
 trame-vtk = "^2.8.15"
 trame-vuetify = "^2.8.1"
 pyvista = "^0.44.2"
 trame-server = "^3.4.0"
+sympy = "^1.13.3"
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.2.0"