Merge pull request #14 from virtualcell/fielddata

add support for field data

Author: jcschaff

examples/scripts/fielddata_workflow.py

@@ -0,0 +1,52 @@
+import os
+import shutil
+from pathlib import Path
+
+from pyvcell.vcml import Simulation, VcmlReader
+from pyvcell.vcml.vcml_simulation import VcmlSpatialSimulation as Solver
+
+# ----- make a workspace
+workspace_dir = Path(os.getcwd()) / "workspace"
+sim1_dir = workspace_dir / "sim1_dir"
+if sim1_dir.exists():
+    shutil.rmtree(sim1_dir)
+sim2_dir = workspace_dir / "sim2_dir"
+if sim2_dir.exists():
+    shutil.rmtree(sim2_dir)
+
+# ---- read in VCML file
+model_fp = Path(os.getcwd()).parent / "models" / "SmallSpatialProject_3D.vcml"
+bio_model1 = VcmlReader.biomodel_from_file(model_fp)
+
+# ---- get the application and the species mappings for species "s0" and "s1"
+app = bio_model1.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")
+
+# ---- 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)
+
+# ---- 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)
+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()
+
+
+# ----- use field data from sim1_dir to set initial concentration of species "s0"
+s0_mapping.init_conc = "vcField('sim1_dir','s0',0.0,'Volume') * vcField('sim1_dir','s1',0.0,'Volume')"
+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.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.4"
+version = "0.1.6"
 description = "This is the python wrapper for vcell modeling and simulation"
 authors = ["Jim Schaff <[email protected]>"]
 repository = "https://github.com/virtualcell/pyvcell"
@@ -33,7 +33,7 @@ matplotlib = "^3.10.0"
 lxml = "^5.3.1"
 imageio = "^2.37.0"
 tensorstore = "^0.1.72"
-libvcell = "^0.0.6"
+libvcell = "^0.0.8"
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.2.0"

pyproject.toml

@@ -13,6 +13,7 @@
 
 PYTHON_ENDIANNESS: Literal["little", "big"] = "big"
 NUMPY_FLOAT_DTYPE = ">f8"
+UNSUPPORTED_VCELL_FUNCTIONS = ["vcField(", "vcRegionVolume(", "vcRegionArea(", "vcConv(", "vcGrad(", "vcProject("]
 
 
 class SpecialLogFileType(Enum):
@@ -290,5 +291,8 @@ def read(self) -> None:
                 _unknown_skipped = parts[2]
                 variable_type = VariableType.from_string(parts[3].strip(" "))
                 _boolean_skipped = parts[4]
+
+                if any(u in expression for u in UNSUPPORTED_VCELL_FUNCTIONS):
+                    continue
                 function = NamedFunction(name=name, vcell_expression=expression, variable_type=variable_type)
                 self.named_functions.append(function)

pyvcell/_internal/simdata/simdata_models.py

@@ -5,6 +5,8 @@
     sbml_spatial_model_3d_path,
 )
 from tests.fixtures.vcell_model_fixtures import (  # noqa: F401
+    vcml_field_data_demo_path,
+    vcml_field_data_tgz_archive_path,
     vcml_spatial_bunny_3d_path,
     vcml_spatial_model_1d_path,
     vcml_spatial_small_3d_path,

tests/conftest.py

@@ -1,3 +1,5 @@
+import tarfile
+import tempfile
 from pathlib import Path
 
 import numpy as np
@@ -271,3 +273,49 @@ def test_vcml_model_parse_3d(vcml_spatial_model_1d_path: Path) -> None:
             dtype=np.float64,
         ),
     )
+
+
+def test_vcml_field_data_demo(vcml_field_data_demo_path: Path, vcml_field_data_tgz_archive_path: Path) -> None:
+    assert vcml_field_data_demo_path.is_file()
+
+    # create a temporary parent directory with a subdirectory for simulation data and a subdirectory to extract the archive
+    parent_dir = Path(tempfile.mkdtemp())
+    sim_dir = parent_dir / "sim"
+    sim_dir.mkdir()
+    assert sim_dir.is_dir()
+
+    with tarfile.open(vcml_field_data_tgz_archive_path, "r:gz") as tar:
+        tar.extractall(parent_dir)
+
+    data_dir = parent_dir / "test2_lsm_DEMO"
+    assert data_dir.is_dir()
+
+    bio_model: Biomodel = VcmlReader.biomodel_from_file(vcml_source=vcml_field_data_demo_path)
+    assert bio_model.model is not None
+    parameters: list[ModelParameter] = bio_model.model.model_parameters
+    assert {p.name: p.value for p in parameters} == {}
+
+    sim_name = bio_model.applications[0].simulations[0].name
+    simulation_orig = VcmlSpatialSimulation(bio_model=bio_model, out_dir=sim_dir)
+    # # create a temporary file to write the VCML content to
+    # vcml_path = Path(tempfile.mktemp(prefix="model_", suffix=".xml"))
+    # vcml_spatial_model.export(vcml_path)
+    results_orig: Result = simulation_orig.run(simulation_name=sim_name)
+
+    channels_orig = results_orig.channel_data
+    assert [channel.label for channel in channels_orig] == [
+        "region_mask",
+        "t",
+        "x",
+        "y",
+        "z",
+        "species0_cyt",
+        "species0_ec",
+    ]
+    # assert np.allclose(
+    #     results_orig.concentrations[0, 0::10],
+    #     np.array(
+    #         [500000.00000000006, 236185.00811512678, 111567.6111715472],
+    #         dtype=np.float64,
+    #     ),
+    # )