Merge pull request #18 from virtualcell/sysbio-course

prepare pyvcell for broader use

Author: jcschaff

examples/scripts/fielddata_from_image_workflow.py

@@ -4,16 +4,19 @@
 
 import numpy as np
 
-from pyvcell.sim_results.var_types import NDArray3D
+from pyvcell._internal.simdata.mesh import CartesianMesh
+from pyvcell.sim_results.var_types import NDArray3D, NDArray4D
 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)
+def create_sinusoid(
+    coords: NDArray4D,
+    freq: float,
+) -> NDArray3D:
     sinusoid: NDArray3D = (
-        np.cos(freq * indices[0, :, :, :]) * np.sin(freq * indices[1, :, :, :]) * np.sin(freq * indices[2, :, :, :])
+        np.cos(freq * coords[:, :, :, 0]) * np.sin(freq * coords[:, :, :, 1]) * np.sin(freq * coords[:, :, :, 2])
     )
     return sinusoid.astype(dtype=np.float64)
 
@@ -47,8 +50,11 @@ def create_sinusoid(shape: tuple[int, ...], freq: float) -> NDArray3D:
     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)
+    coords_array = CartesianMesh.compute_coordinates(
+        mesh_shape=(shape[0], shape[1], shape[2]), origin=app.geometry.origin, extent=app.geometry.extent
+    )
+    fields[0].data_nD = np.multiply(create_sinusoid(coords=coords_array, freq=0.5), 8.0)
+    fields[1].data_nD = np.multiply(create_sinusoid(coords=coords_array, freq=0.3), 4.0)
 
     # ---- add field data to the simulation
 

examples/scripts/sysbio_class.py

@@ -0,0 +1,36 @@
+from pathlib import Path
+
+import pyvcell.vcml as vcml
+from pyvcell.vcml.vcml_simulation import VcmlSpatialSimulation as VCSolver
+
+WORKSPACE_DIR = Path(__file__).parent.resolve()
+
+antimony_str = """
+    // This is a simple example of a spatial model in VCell
+    compartment extracellular;
+    compartment cell;
+    species A in cell;
+    species B in cell;
+    A -> B; k1*A - k2*B
+
+    k1 = 0.1; k2 = 0.2
+    A = 10
+"""
+
+biomodel = vcml.load_antimony_str(antimony_str)
+
+# create a 3D geometry
+geo = vcml.Geometry(name="geo", origin=(0, 0, 0), extent=(10, 10, 10), dim=3)
+# cell = geo.add_sphere(name="cell", radius=4, center=(5,5,5))
+medium = geo.add_background(name="background")
+
+# add an application to the biomodel
+app = biomodel.add_application("app1", geometry=geo)
+app.map_compartment(compartment=biomodel.model.get_compartment("cell"), domain=medium)  # type: ignore[union-attr]
+app.species_mappings = [
+    vcml.SpeciesMapping(species_name="A", init_conc="sin(x)"),
+    vcml.SpeciesMapping(species_name="B", init_conc="cos(x+y+z)"),
+]
+sim = app.add_sim(name="sim1", duration=10.0, output_time_step=0.5, mesh_size=(50, 50, 50))
+results = VCSolver(bio_model=biomodel, out_dir=WORKSPACE_DIR / "sim1").run(simulation_name="sim1")
+results.plotter.plot_concentrations()

poetry.lock

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pyvcell"
-version = "0.1.10"
+version = "0.1.11"
 description = "This is the python wrapper for vcell modeling and simulation"
 authors = ["Jim Schaff <[email protected]>"]
 repository = "https://github.com/virtualcell/pyvcell"
@@ -32,13 +32,14 @@ matplotlib = "^3.10.0"
 lxml = "^5.3.1"
 imageio = "^2.37.0"
 tensorstore = "^0.1.72"
-libvcell = "^0.0.12"
+libvcell = "^0.0.13"
 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.1"
+antimony = "^2.15.0"
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.2.0"

pyproject.toml

@@ -5,6 +5,7 @@
 from numpy._typing import NDArray
 
 from pyvcell._internal.simdata.vtk.vismesh import Box3D
+from pyvcell.sim_results.var_types import NDArray4D
 
 
 class CartesianMesh:
@@ -103,6 +104,35 @@ def dimension(self) -> int:
         else:
             return 3
 
+    @property
+    def coordinates(self) -> NDArray[np.float64]:
+        """
+        Returns the coordinates of the mesh as a 3D array of shape (size[0], size[1], size[2], 3)
+        """
+        return self.compute_coordinates(mesh_shape=self.size, origin=self.origin, extent=self.extent)
+
+    @staticmethod
+    def compute_coordinates(
+        mesh_shape: tuple[int, int, int] | list[int],
+        origin: tuple[float, float, float] | list[float],
+        extent: tuple[float, float, float] | list[float],
+    ) -> NDArray4D:
+        """
+        Computes the coordinates of the mesh based on the mesh shape, origin, and extent.
+        Returns a 3D array of shape (size[0], size[1], size[2], 3)
+        """
+        coords = np.empty((*mesh_shape, 3), dtype=np.float64)
+        coords[..., 0] = origin[0] + np.arange(mesh_shape[0])[:, None, None] * extent[0] / (mesh_shape[0] - 1)
+        if mesh_shape[1] > 1:
+            coords[..., 1] = origin[1] + np.arange(mesh_shape[1])[None, :, None] * extent[1] / (mesh_shape[1] - 1)
+        else:
+            coords[..., 1] = origin[1]
+        if mesh_shape[2] > 1:
+            coords[..., 2] = origin[2] + np.arange(mesh_shape[2])[None, None, :] * extent[2] / (mesh_shape[2] - 1)
+        else:
+            coords[..., 2] = origin[2]
+        return coords
+
     def read(self) -> None:
         # read file as lines and parse
         with self.mesh_file.open("r") as f:

pyvcell/_internal/simdata/mesh.py

@@ -5,6 +5,7 @@
 NDArray1D: TypeAlias = np.ndarray[tuple[int], np.dtype[np.float64]]
 NDArray2D: TypeAlias = np.ndarray[tuple[int, int], np.dtype[np.float64]]
 NDArray3D: TypeAlias = np.ndarray[tuple[int, int, int], np.dtype[np.float64]]
+NDArray4D: TypeAlias = np.ndarray[tuple[int, int, int, int], np.dtype[np.float64]]
 NDArrayND: TypeAlias = np.ndarray[tuple[int, ...], np.dtype[np.float64]]
 
 NDArray1Du8: TypeAlias = np.ndarray[tuple[int], np.dtype[np.uint8]]

pyvcell/sim_results/var_types.py

@@ -1,3 +1,8 @@
+import logging
+import tempfile
+from os import PathLike
+from pathlib import Path
+
 from pyvcell.vcml.models import (
     Application,
     Biomodel,
@@ -50,3 +55,71 @@
     "Application",
     "Simulation",
 ]
+
+
+def load_antimony_str(antimony_str: str) -> Biomodel:
+    import antimony  # type: ignore[import-untyped]
+
+    antimony_success = antimony.loadAntimonyString(antimony_str)
+    if antimony_success != -1:
+        sbml_str = antimony.getSBMLString()
+        sbml_str = sbml_str.replace("sboTerm", "metaid")
+        logging.info(f"Hack - introduced a metaid in place of sboTerm to SBML string:\n{sbml_str}")
+        return load_sbml_str(sbml_str)
+    else:
+        raise ValueError("Error loading model:", antimony.getLastError())
+
+
+def load_antimony_file(antimony_file: PathLike[str]) -> Biomodel:
+    import antimony  # ignore
+
+    antimony_success = antimony.loadAntimonyFile(antimony_file)
+    if antimony_success != -1:
+        sbml_str = antimony.getSBMLString()
+        return load_sbml_str(sbml_str)
+    else:
+        raise ValueError("Error loading model:", antimony.getLastError())
+
+
+def load_sbml_file(sbml_file: PathLike[str]) -> Biomodel:
+    import libvcell
+
+    with tempfile.TemporaryDirectory() as tempdir:
+        with open(sbml_file) as f:
+            sbml_str = f.read()
+        vcml_path = Path(tempdir) / "model.vcml"
+        vc_success, vc_errmsg = libvcell.sbml_to_vcml(sbml_content=sbml_str, vcml_file_path=vcml_path)
+        if vc_success:
+            return VcmlReader.biomodel_from_file(vcml_path=vcml_path)
+        else:
+            raise ValueError("Error loading model:", vc_errmsg)
+
+
+def load_sbml_str(sbml_str: str) -> Biomodel:
+    import libvcell
+
+    with tempfile.TemporaryDirectory() as tempdir:
+        vcml_path = Path(tempdir) / "model.vcml"
+        vc_success, vc_errmsg = libvcell.sbml_to_vcml(sbml_content=sbml_str, vcml_file_path=vcml_path)
+        if vc_success:
+            return VcmlReader.biomodel_from_file(vcml_path=vcml_path)
+        else:
+            raise ValueError("Error loading model:", vc_errmsg)
+
+
+def write_vcml_file(bio_model: Biomodel, vcml_file: PathLike[str]) -> None:
+    vcml_document = VCMLDocument(biomodel=bio_model)
+    VcmlWriter.write_to_file(vcml_document=vcml_document, file_path=vcml_file)
+
+
+def to_vcml_str(bio_model: Biomodel) -> str:
+    vcml_document = VCMLDocument(biomodel=bio_model)
+    vcml_str: str = VcmlWriter().write_vcml(document=vcml_document)
+    return vcml_str
+
+
+def refresh_biomodel(bio_model: Biomodel) -> Biomodel:
+    with tempfile.TemporaryDirectory() as tempdir:
+        vcml_path = Path(tempdir) / "model.vcml"
+        write_vcml_file(bio_model=bio_model, vcml_file=vcml_path)
+        return VcmlReader.biomodel_from_file(vcml_path=vcml_path)

pyvcell/vcml/__init__.py

@@ -8,6 +8,6 @@ cd ${ROOT_DIR} || (echo "Failed to cd into ${ROOT_DIR}" && exit 1)
 
 python -m venv .venv_jupyter
 source .venv_jupyter/bin/activate
-pip install jupyter ipython matplotlib tensorstore
+pip install jupyter ipython
 poetry install
 jupyter notebook

scripts/run_notebook.sh

@@ -1,5 +1,7 @@
 from pathlib import Path
 
+import numpy as np
+
 from pyvcell._internal.simdata.mesh import CartesianMesh
 from pyvcell._internal.simdata.vtk.fv_mesh_mapping import from_mesh3d_membrane, from_mesh3d_volume
 from pyvcell._internal.simdata.vtk.vtkmesh_fv import write_finite_volume_smoothed_vtk_grid_and_index_data
@@ -46,3 +48,45 @@ def test_mesh_parse(temp_sim_946368938_path: Path) -> None:
         "cytosol.json",
         "SimID_946368938_0_.functions",
     }
+
+
+def test_mesh_coords(temp_sim_946368938_path: Path) -> None:
+    input_filenames = {"SimID_946368938_0_.mesh", "SimID_946368938_0_.functions"}
+    # remove all files not in input_filenames
+    for p in temp_sim_946368938_path.iterdir():
+        if p.name not in input_filenames:
+            p.unlink()
+
+    filenames = {p.name for p in temp_sim_946368938_path.iterdir()}
+    assert filenames == input_filenames
+
+    mesh = CartesianMesh(mesh_file=temp_sim_946368938_path / "SimID_946368938_0_.mesh")
+    mesh.read()
+
+    assert mesh.size == [71, 71, 25]
+    assert mesh.dimension == 3
+
+    # method 1 - using mesh property
+    coords_1 = mesh.coordinates
+    assert coords_1.shape == (71, 71, 25, 3)
+    assert np.all(coords_1[0, 0, 0, :] == tuple(mesh.origin))
+    assert np.all(
+        coords_1[-1, -1, -1, :]
+        == (mesh.origin[0] + mesh.extent[0], mesh.origin[1] + mesh.extent[1], mesh.origin[2] + mesh.extent[2])
+    )
+
+    # method 2 - using static method
+    coords = CartesianMesh.compute_coordinates(
+        (mesh.size[0], mesh.size[1], mesh.size[2]),
+        (mesh.origin[0], mesh.origin[1], mesh.origin[2]),
+        (mesh.extent[0], mesh.extent[1], mesh.extent[2]),
+    )
+    assert coords.shape == (71, 71, 25, 3)
+    assert np.all(coords[0, 0, 0, :] == tuple(mesh.origin))
+    assert np.all(
+        coords[-1, -1, -1, :]
+        == (mesh.origin[0] + mesh.extent[0], mesh.origin[1] + mesh.extent[1], mesh.origin[2] + mesh.extent[2])
+    )
+
+    # check that the two methods give the same result
+    assert np.all(coords == coords_1)