Extend support for more features of VCell expressions

.github/actions/setup-poetry-env/action.yml

@@ -17,7 +17,7 @@ runs:
 
     - name: Install Poetry
       env:
-        POETRY_VERSION: "1.7.1"
+        POETRY_VERSION: "1.8.5"
       run: curl -sSL https://install.python-poetry.org | python - -y
       shell: bash
 

.gitignore

@@ -182,3 +182,4 @@ examples/test_output/SimID_*
 examples/notebooks/workspace/
 
 examples/scripts/workspace/
+workspace

examples/scripts/_internal_data_demo.py

@@ -94,7 +94,7 @@
 # Set labels for axes
 ax.set_xlabel("X")
 ax.set_ylabel("Y")
-ax.set_zlabel("Z")  # type: ignore[attr-defined]
+ax.set_zlabel("Z")
 
 # Show the plot
 plt.show()
@@ -138,7 +138,7 @@
 # plot each series on a different plot arranged in a 2x2 grid with a legend from series_legends
 fig, ax = plt.subplots(2, 2, figsize=(10, 10))
 for i, series_array in enumerate(series_arrays):
-    axis = ax[int(i / 2), i % 2]  # type: ignore[index]
+    axis = ax[int(i / 2), i % 2]
     axis.plot(times, series_array[:, 0], label="mean")
     axis.fill_between(times, series_array[:, 1], series_array[:, 2], alpha=0.2)
     axis.set_title(series_legend[i])

examples/scripts/_internal_n5_download_demo.py

@@ -1,12 +1,12 @@
 import matplotlib.pyplot as plt
-from tensorstore._tensorstore import TensorStore  # type: ignore[import-untyped]
+import tensorstore._tensorstore as ts  # type: ignore[import-not-found]
 
 from pyvcell._internal.simdata.n5_data import vcell_n5_datastore
 from pyvcell.sim_results.var_types import NDArray2D
 
 url = "https://vcell-dev.cam.uchc.edu/n5Data/ACowan/4b5ac930c40d5ba.n5"
 dataset_name = "4248805214"
-data_store: TensorStore = vcell_n5_datastore(base_url=url, dataset_name=dataset_name)
+data_store: ts.TensorStore = vcell_n5_datastore(base_url=url, dataset_name=dataset_name)
 print(f"shape = {data_store.shape}")
 
 np_data: NDArray2D = data_store[:, :, 0, 0, 0].read().result()

pyproject.toml

@@ -32,7 +32,7 @@ matplotlib = "^3.10.0"
 lxml = "^5.3.1"
 imageio = "^2.37.0"
 tensorstore = "^0.1.72"
-libvcell = "^0.0.13"
+libvcell = "^0.0.15"
 trame = "^3.8.1"
 trame-vtk = "^2.8.15"
 trame-vuetify = "^2.8.1"

pyvcell/_internal/simdata/n5_data.py

@@ -1,7 +1,6 @@
-import tensorstore as ts  # type: ignore[import-untyped]
-from tensorstore._tensorstore import TensorStore  # type: ignore[import-untyped]
+import tensorstore as ts
 
 
-def vcell_n5_datastore(base_url: str, dataset_name: str) -> TensorStore:
+def vcell_n5_datastore(base_url: str, dataset_name: str) -> ts._tensorstore.TensorStore:
     spec = {"driver": "n5", "kvstore": {"driver": "http", "base_url": base_url}, "path": dataset_name}
     return ts.open(spec, read=True).result()

pyvcell/_internal/simdata/python_infix.py

@@ -0,0 +1,23 @@
+from libvcell import vcell_infix_to_num_expr_infix
+
+
+def get_numexpr_expression(vcell_expression: str) -> str:
+    # First, get the string from vcell
+    translation_result: bool
+    result_message: str
+    num_expr_expression: str
+    translation_result, result_message, num_expr_expression = vcell_infix_to_num_expr_infix(vcell_expression)
+    if not translation_result:
+        raise ValueError(
+            result_message if result_message is not None else "Unable to convert expression: " + vcell_expression
+        )
+    # NumExpr has a restricted set of what's allowed, in order to protect `eval`
+    sanitized_expression: str = (
+        num_expr_expression.lstrip(" ")
+        .rstrip(" ")
+        .replace(" and ", " & ")
+        .replace(" or ", " | ")
+        .replace(" not ", " ~")
+        .replace("math.", "")
+    )
+    return sanitized_expression

pyvcell/_internal/simdata/simdata_models.py

@@ -9,6 +9,7 @@
 import numpy as np
 from numpy._typing import NDArray
 
+from pyvcell._internal.simdata.python_infix import get_numexpr_expression
 from pyvcell.sim_results.var_types import NDArray1D
 
 PYTHON_ENDIANNESS: Literal["little", "big"] = "big"
@@ -350,23 +351,23 @@ def get_data(self, variable: VariableInfo | str, time: float) -> NDArray1D:
 class NamedFunction:
     name: str
     vcell_expression: str
-    python_expression: str
+    num_expr_expression: str
     variables: list[str]
     variable_type: VariableType
 
     def __init__(self, name: str, vcell_expression: str, variable_type: VariableType) -> None:
         self.name = name
         self.vcell_expression = vcell_expression
-        self.python_expression = vcell_expression.replace("^", "**").lstrip(" ").rstrip(" ")
+        self.num_expr_expression = get_numexpr_expression(self.vcell_expression)
         self.variable_type = variable_type
 
         # Parse the python expression into an AST and extract all Name nodes (which represent variables)
-        tree = ast.parse(self.python_expression)
+        tree = ast.parse(self.num_expr_expression)
         self.variables = [node.id for node in ast.walk(tree) if isinstance(node, ast.Name)]
 
     def evaluate(self, variable_bindings: dict[str, NDArray[np.float64]]) -> NDArray[np.float64]:
         ne.set_num_threads(1)
-        expression = self.python_expression
+        expression = self.num_expr_expression
         result = ne.evaluate(expression, local_dict=variable_bindings)
         if not isinstance(result, np.ndarray):
             raise TypeError(f"Expression {expression} did not evaluate to a numpy array")

pyvcell/_internal/simdata/zarr_writer.py

@@ -144,9 +144,9 @@ def write_zarr(pde_dataset: PdeDataSet, data_functions: DataFunctions, mesh: Car
                     "max_values": [],
                     "mean_values": [],
                 })
-            channel_metadata[c]["min_values"].append(np.min(func_data))
-            channel_metadata[c]["max_values"].append(np.max(func_data))
-            channel_metadata[c]["mean_values"].append(np.mean(func_data))
+            channel_metadata[c]["min_values"].append(float(np.min(func_data)))
+            channel_metadata[c]["max_values"].append(float(np.max(func_data)))
+            channel_metadata[c]["mean_values"].append(float(np.mean(func_data)))
             c = c + 1
 
     z1.attrs["metadata"] = {

pyvcell/sim_results/plotter.py

@@ -123,7 +123,7 @@ def plot_slice_3d(self, time_index: int, channel_id: str) -> None:
         # Set labels for axes
         ax.set_xlabel("X")
         ax.set_ylabel("Y")
-        ax.set_zlabel("Z")  # type: ignore[attr-defined]
+        ax.set_zlabel("Z")
         t = self.times[time_index]
         title = f"{channel.label} (in {channel.domain_name}) at t={t}"
         plt.title(title)
@@ -159,7 +159,7 @@ def get_3d_slice_animation(self, channel_id: str, interval: int = 200) -> animat
         # Set labels for axes
         ax.set_xlabel("X")
         ax.set_ylabel("Y")
-        ax.set_zlabel("Z")  # type: ignore[attr-defined]
+        ax.set_zlabel("Z")
         sc = None
 
         def update(frame: int) -> tuple[PathCollection]:

tests/_internal/simdata/vcell_parsing_test.py

@@ -143,7 +143,7 @@ def test_function_eval(temp_sim_946368938_path: Path) -> None:
     assert function_J_r0.name.split("::")[1] == "J_r0"
 
     assert function_J_r0.variables == ["RanC_cyt", "Ran_cyt", "C_cyt"]
-    assert function_J_r0.python_expression == "(RanC_cyt - (1000.0 * C_cyt * Ran_cyt))"
+    assert function_J_r0.num_expr_expression == "(RanC_cyt - (1000.0 * C_cyt * Ran_cyt))"
     min_values = []
     max_values = []
     for t in pde_dataset.times():

tests/conftest.py

@@ -1,3 +1,7 @@
+from pathlib import Path
+
+import pytest
+
 from tests.fixtures.data_fixtures import (  # noqa: F401
     fielddata_file_path,
     solver_output_path,
@@ -21,3 +25,13 @@
     vcml_spatial_small_3d_path,
     vcml_tutorial_multiapp_pde_path,
 )
+
+
+@pytest.fixture
+def test_root_dir() -> Path:
+    return Path(__file__).parent
+
+
+@pytest.fixture
+def fixture_root_dir(test_root_dir: Path) -> Path:
+    return test_root_dir / "fixtures"

tests/expressions/test_expressions.py

@@ -0,0 +1,27 @@
+import numexpr  # type: ignore[import-untyped]
+from numpy import ndarray
+
+from pyvcell._internal.simdata.python_infix import get_numexpr_expression
+
+
+def test_vcell_infix_to_num_expr_infix() -> None:
+    vcell_infix = "id_1 * csc(id_0 ^ 2.2)"
+    value = get_numexpr_expression(vcell_infix)
+    assert value == "(id_1 * (1.0/sin(((id_0)**(2.2)))))"
+
+
+def test_bad_vcell_infix() -> None:
+    vcellInfix = "id_1 / + / /-  cos(/ / /) id_2"
+    try:
+        get_numexpr_expression(vcellInfix)
+    except ValueError:
+        return
+    raise AssertionError("Should have raised ValueError")
+
+
+def test_if_else_in_expression() -> None:
+    vcell_infix = "(id_2 && id_3) * 1.2"
+    value = get_numexpr_expression(vcell_infix)
+    assert value == "(where(((0.0!=id_2) & (0.0!=id_3)), 1.2, 0.0))"
+    result: ndarray = ndarray(1)
+    numexpr.evaluate(value, {"id_2": 2.2, "id_3": 3.3, "float": float, bool: bool}, out=result)

tests/fixtures/data/TinySpatialLogicProject.vcml

@@ -0,0 +1,152 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--This biomodel was generated in VCML Version Alpha_Version_7.7.0_build_49-->
+<vcml xmlns="http://sourceforge.net/projects/vcell/vcml" Version="Alpha_Version_7.7.0_build_49">
+  <BioModel Name="TinySpatialProject_Application0">
+    <Model Name="unnamed">
+      <ModelParameters>
+        <Parameter Name="Kf_r0" Role="user defined" Unit="s-1">(ceil(sin(t)) &amp;&amp; s0)</Parameter>
+        <Parameter Name="Kr_r0" Role="user defined" Unit="s-1">0.5</Parameter>
+      </ModelParameters>
+      <Compound Name="s0">
+        <Annotation>s0</Annotation>
+      </Compound>
+      <Compound Name="s1">
+        <Annotation>s1</Annotation>
+      </Compound>
+      <Feature Name="c0" />
+      <LocalizedCompound Name="s0" SbmlName="s0" CompoundRef="s0" Structure="c0" OverrideName="true" />
+      <LocalizedCompound Name="s1" SbmlName="s1" CompoundRef="s1" Structure="c0" OverrideName="true" />
+      <SimpleReaction Structure="c0" Name="r0" Reversible="true" FluxOption="MolecularOnly" SbmlName="r0">
+        <Reactant LocalizedCompoundRef="s0" Stoichiometry="1" />
+        <Product LocalizedCompoundRef="s1" Stoichiometry="1" />
+        <Kinetics KineticsType="GeneralKinetics">
+          <Parameter Name="J" Role="reaction rate" Unit="umol.l-1.s-1">((Kf_r0 * s0) - (Kr_r0 * s1))</Parameter>
+        </Kinetics>
+      </SimpleReaction>
+      <Diagram Name="c0" Structure="c0">
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="470" LocationY="220" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="244" LocationY="129" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="14" LocationY="34" />
+      </Diagram>
+      <ModelUnitSystem VolumeSubstanceUnit="umol" MembraneSubstanceUnit="umol" LumpedReactionSubstanceUnit="umol" VolumeUnit="l" AreaUnit="dm2" LengthUnit="dm" TimeUnit="s" />
+    </Model>
+    <SimulationSpec Name="unnamed_spatialGeom" Stochastic="false" UseConcentration="true" SpringSaLaD="false" RuleBased="false" MassConservationModelReduction="false" InsufficientIterations="false" InsufficientMaxMolecules="false" CharacteristicSize="0.021090187684898395">
+      <NetworkConstraints RbmMaxIteration="1" RbmMaxMoleculesPerSpecies="10" RbmSpeciesLimit="800" RbmReactionsLimit="2500" />
+      <Annotation />
+      <Geometry Name="spatialGeom" Dimension="1">
+        <Extent X="10.0" Y="1.0" Z="1.0" />
+        <Origin X="0.0" Y="0.0" Z="0.0" />
+        <SubVolume Name="subdomain0" Handle="0" Type="Analytical">
+          <AnalyticExpression>1.0</AnalyticExpression>
+        </SubVolume>
+        <SurfaceDescription NumSamplesX="50" NumSamplesY="1" NumSamplesZ="1" CutoffFrequency="0.3">
+          <VolumeRegion Name="subdomain00" RegionID="0" SubVolume="subdomain0" Size="10.0" Unit="um" />
+        </SurfaceDescription>
+      </Geometry>
+      <GeometryContext>
+        <FeatureMapping Feature="c0" GeometryClass="subdomain0" SubVolume="subdomain0" Size="50000.0" VolumePerUnitVolume="1.0">
+          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
+        </FeatureMapping>
+      </GeometryContext>
+      <ReactionContext>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s0" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(100000.0 * x)</InitialConcentration>
+          <Diffusion>1.0E-9</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" />
+        </LocalizedCompoundSpec>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s1" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(10.0 - (1.0 * 100000.0 * x))</InitialConcentration>
+          <Diffusion>1.0E-9</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" />
+        </LocalizedCompoundSpec>
+        <ReactionSpec ReactionStepRef="r0" ReactionMapping="included" />
+      </ReactionContext>
+      <MathDescription Name="unnamed_spatialGeom_generated">
+        <Constant Name="_F_">96485.3321</Constant>
+        <Constant Name="_F_nmol_">9.64853321E-5</Constant>
+        <Constant Name="_K_GHK_">1.0E-9</Constant>
+        <Constant Name="_N_pmol_">6.02214179E11</Constant>
+        <Constant Name="_PI_">3.141592653589793</Constant>
+        <Constant Name="_R_">8314.46261815</Constant>
+        <Constant Name="_T_">300.0</Constant>
+        <Constant Name="K_millivolts_per_volt">1000.0</Constant>
+        <Constant Name="KMOLE">0.001660538783162726</Constant>
+        <Constant Name="Kr_r0">0.5</Constant>
+        <Constant Name="s0_boundaryXm">0.0</Constant>
+        <Constant Name="s0_boundaryXp">0.0</Constant>
+        <Constant Name="s0_diffusionRate">1.0E-9</Constant>
+        <Constant Name="s1_boundaryXm">0.0</Constant>
+        <Constant Name="s1_boundaryXp">0.0</Constant>
+        <Constant Name="s1_diffusionRate">1.0E-9</Constant>
+        <Constant Name="VolumePerUnitVolume_c0">1.0</Constant>
+        <VolumeVariable Name="s0" Domain="subdomain0" />
+        <VolumeVariable Name="s1" Domain="subdomain0" />
+        <Function Name="J_r0" Domain="subdomain0">((Kf_r0 * s0) - (Kr_r0 * s1))</Function>
+        <Function Name="Kf_r0" Domain="subdomain0">(ceil(sin(t)) &amp;&amp; s0)</Function>
+        <Function Name="s0_init_umol_l_1" Domain="subdomain0">(100000.0 * x)</Function>
+        <Function Name="s1_init_umol_l_1" Domain="subdomain0">(10.0 - (1.0 * 100000.0 * x))</Function>
+        <Function Name="Size_c0" Domain="subdomain0">(VolumePerUnitVolume_c0 * vcRegionVolume('subdomain0'))</Function>
+        <Function Name="vobj_subdomain00_size" Domain="subdomain0">vcRegionVolume('subdomain0')</Function>
+        <CompartmentSubDomain Name="subdomain0">
+          <BoundaryType Boundary="Xm" Type="Flux" />
+          <BoundaryType Boundary="Xp" Type="Flux" />
+          <BoundaryType Boundary="Ym" Type="Value" />
+          <BoundaryType Boundary="Yp" Type="Value" />
+          <BoundaryType Boundary="Zm" Type="Value" />
+          <BoundaryType Boundary="Zp" Type="Value" />
+          <PdeEquation Name="s0" SolutionType="Unknown">
+            <Boundaries Xm="s0_boundaryXm" Xp="s0_boundaryXp" />
+            <Rate>- J_r0</Rate>
+            <Diffusion>s0_diffusionRate</Diffusion>
+            <Initial>s0_init_umol_l_1</Initial>
+          </PdeEquation>
+          <PdeEquation Name="s1" SolutionType="Unknown">
+            <Boundaries Xm="s1_boundaryXm" Xp="s1_boundaryXp" />
+            <Rate>J_r0</Rate>
+            <Diffusion>s1_diffusionRate</Diffusion>
+            <Initial>s1_init_umol_l_1</Initial>
+          </PdeEquation>
+        </CompartmentSubDomain>
+      </MathDescription>
+      <Simulation Name="Simulation0">
+        <SolverTaskDescription TaskType="Unsteady" UseSymbolicJacobian="false" Solver="Sundials Stiff PDE Solver (Variable Time Step)">
+          <TimeBound StartTime="0.0" EndTime="1.0" />
+          <TimeStep DefaultTime="0.05" MinTime="0.0" MaxTime="0.1" />
+          <ErrorTolerance Absolut="1.0E-9" Relative="1.0E-7" />
+          <OutputOptions OutputTimeStep="0.05" />
+          <SundialsSolverOptions>
+            <maxOrderAdvection>2</maxOrderAdvection>
+          </SundialsSolverOptions>
+          <NumberProcessors>1</NumberProcessors>
+        </SolverTaskDescription>
+        <MathOverrides />
+        <MeshSpecification>
+          <Size X="200" Y="1" Z="1" />
+        </MeshSpecification>
+      </Simulation>
+      <SpatialObjects>
+        <SpatialObject Name="vobj_subdomain00" Type="Volume" subVolume="subdomain0" regionId="0">
+          <QuantityCategoryList>
+            <QuantityCategory Name="VolumeCentroid" Enabled="false" />
+            <QuantityCategory Name="InteriorVelocity" Enabled="false" />
+            <QuantityCategory Name="VolumeRegionSize" Enabled="true" />
+          </QuantityCategoryList>
+        </SpatialObject>
+      </SpatialObjects>
+      <MicroscopeMeasurement Name="fluor">
+        <ConvolutionKernel Type="ProjectionZKernel" />
+      </MicroscopeMeasurement>
+    </SimulationSpec>
+    <pathwayModel>
+      <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:bp="http://www.biopax.org/release/biopax-level3.owl#" version="3.0" />
+    </pathwayModel>
+    <relationshipModel>
+      <RMNS version="3.0" />
+    </relationshipModel>
+    <vcmetadata>
+      <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" />
+      <nonrdfAnnotationList />
+      <uriBindingList />
+    </vcmetadata>
+  </BioModel>
+</vcml>

tests/sim_results/test_model_with_logic_in_expressions.py

@@ -0,0 +1,17 @@
+from pathlib import Path
+
+import pyvcell.vcml as vcml
+
+
+def test_model_with_logic_in_expressions(fixture_root_dir: Path) -> None:
+    model_fp = fixture_root_dir / "data" / "TinySpatialLogicProject.vcml"
+    assert model_fp.exists()
+    bio_model = vcml.load_vcml_file(model_fp)
+    for sim in [sim for app in bio_model.applications for sim in app.simulations]:
+        result = vcml.simulate(biomodel=bio_model, simulation=sim.name)
+
+        print([c.label for c in result.channel_data])
+
+        result.plotter.plot_slice_3d(time_index=3, channel_id="s1")
+        result.plotter.plot_concentrations()
+        result.cleanup()