Merge pull request #48 from virtualcell/support-velocity

Add advection velocity to species mapping + 'user defined' parameter role

Author: jcschaff

pyvcell/vcml/models.py

@@ -162,8 +162,8 @@ def add_species(self, name: str, compartment: str | Compartment) -> Species:
         self.species.append(species)
         return species
 
-    def add_model_parameter(self, name: str, value: float | str) -> ModelParameter:
-        model_parameter = ModelParameter(name=name, value=value, role="model_parameter", unit="")
+    def add_model_parameter(self, name: str, value: float | str, role: str = "user defined") -> ModelParameter:
+        model_parameter = ModelParameter(name=name, value=value, role=role, unit="")
         self.model_parameters.append(model_parameter)
         return model_parameter
 

pyvcell/vcml/models_app.py

@@ -56,6 +56,9 @@ class SpeciesMapping(VcmlNode):
     init_conc: float | str | None = None
     init_count: float | str | None = None
     diff_coef: float | str | None = None
+    velocity_x: float | str | None = None
+    velocity_y: float | str | None = None
+    velocity_z: float | str | None = None
     boundary_values: list[float | str | None] = Field(default_factory=list)
 
     @property
@@ -67,6 +70,9 @@ def expressions(self) -> list[str]:
             exps.append(self.init_count)
         if isinstance(self.diff_coef, str):
             exps.append(self.diff_coef)
+        for velocity in (self.velocity_x, self.velocity_y, self.velocity_z):
+            if isinstance(velocity, str):
+                exps.append(velocity)
         if self.boundary_values:
             for value in self.boundary_values:
                 if isinstance(value, str):

pyvcell/vcml/vcml_reader.py

@@ -70,10 +70,10 @@ def print_biomodel(cls, xml_string: str) -> None:
 class XMLVisitor:
     def visit(self, element: _Element, node: vc.VcmlNode) -> None:
         method_name = "visit_" + strip_namespace(element.tag)
-        method = getattr(self, method_name, self.generic_visit)
+        method = getattr(self, method_name, self.generic_visit_children)
         method(element=element, node=node)
 
-    def generic_visit(self, element: _Element, node: vc.VcmlNode) -> None:
+    def generic_visit_children(self, element: _Element, node: vc.VcmlNode) -> None:
         for child in element:
             # Best-effort parse: an unmodeled or malformed physiology subtree should
             # cost us that subtree, not abort the whole document (and in particular not
@@ -93,7 +93,7 @@ def __init__(self, document: vc.VCMLDocument) -> None:
     def visit_BioModel(self, element: _Element, node: vc.VCMLDocument) -> None:
         name = element.get("Name", default="unnamed")
         node.biomodel = vc.Biomodel(name=name, version=self._parse_version(element))
-        self.generic_visit(element, node.biomodel)
+        self.generic_visit_children(element, node.biomodel)
 
     def _parse_version(self, element: _Element) -> vc.Version | None:
         """Extract a Version child element, if present."""
@@ -118,21 +118,21 @@ def _parse_version(self, element: _Element) -> vc.Version | None:
     def visit_Model(self, element: _Element, node: vc.Biomodel) -> None:
         name: str = element.get("Name", default="unnamed")
         node.model = vc.Model(name=name)
-        self.generic_visit(element, node.model)
+        self.generic_visit_children(element, node.model)
 
     def visit_SimpleReaction(self, element: _Element, node: vc.Model) -> None:
         name: str = element.get("Name", default="unnamed")
         compartment_name: str = element.get("Structure", default="unknown")
         reaction = vc.Reaction(name=name, is_flux=False, compartment_name=compartment_name)
         node.reactions.append(reaction)
-        self.generic_visit(element, reaction)
+        self.generic_visit_children(element, reaction)
 
     def visit_FluxStep(self, element: _Element, node: vc.Model) -> None:
         name: str = element.get("Name", default="unnamed")
         compartment_name: str = element.get("Structure", default="unknown")
         reaction = vc.Reaction(name=name, is_flux=True, compartment_name=compartment_name)
         node.reactions.append(reaction)
-        self.generic_visit(element, reaction)
+        self.generic_visit_children(element, reaction)
 
     def visit_Reactant(self, element: _Element, node: vc.Reaction) -> None:
         compound_ref: str = element.get("LocalizedCompoundRef", default="unknown")
@@ -141,7 +141,7 @@ def visit_Reactant(self, element: _Element, node: vc.Reaction) -> None:
             name=compound_ref, stoichiometry=stoichiometry, species_ref_type=vc.SpeciesRefType.reactant
         )
         node.reactants.append(reaction)
-        self.generic_visit(element, reaction)
+        self.generic_visit_children(element, reaction)
 
     def visit_Product(self, element: _Element, node: vc.Reaction) -> None:
         compound_ref: str = element.get("LocalizedCompoundRef", default="unknown")
@@ -150,7 +150,7 @@ def visit_Product(self, element: _Element, node: vc.Reaction) -> None:
             name=compound_ref, stoichiometry=stoichiometry, species_ref_type=vc.SpeciesRefType.product
         )
         node.products.append(reaction)
-        self.generic_visit(element, reaction)
+        self.generic_visit_children(element, reaction)
 
     def visit_Kinetics(self, element: _Element, node: vc.VcmlNode) -> None:
         # Only attach kinetics to an actual reaction. A <Kinetics> reached with a
@@ -160,7 +160,7 @@ def visit_Kinetics(self, element: _Element, node: vc.VcmlNode) -> None:
         kinetics_type: str = element.get("KineticsType", default="GeneralKinetics")
         kinetics = vc.Kinetics(kinetics_type=kinetics_type)
         node.kinetics = kinetics
-        self.generic_visit(element, kinetics)
+        self.generic_visit_children(element, kinetics)
 
     def visit_Feature(self, element: _Element, node: vc.Model) -> None:
         name = element.get("Name", default="unnamed")
@@ -209,15 +209,15 @@ def visit_Parameter(self, element: _Element, node: vc.VcmlNode) -> None:
             # A <Parameter> in a context the data model doesn't represent (rate rules,
             # structure/species-context mappings, electrical params, …) — skip it.
             return
-        self.generic_visit(element, parameter)
+        self.generic_visit_children(element, parameter)
 
     def visit_SimulationSpec(self, element: _Element, node: vc.Biomodel) -> None:
         name: str = element.get("Name", default="unnamed")
         stochastic: bool = element.get("Stochastic", default="false").lower() == "true"
         default_geometry = vcg.Geometry(name="default", dim=3)
         application = vc.Application(name=name, stochastic=stochastic, geometry=default_geometry)
         node.applications.append(application)
-        self.generic_visit(element, application)
+        self.generic_visit_children(element, application)
 
     def visit_Simulation(self, element: _Element, node: vc.Application) -> None:
         name: str = element.get("Name", default="unnamed")
@@ -468,7 +468,7 @@ def visit_Geometry(self, element: _Element, node: vc.Application) -> None:
         dim = int(element.get("Dimension", default="0"))
         geometry = vcg.Geometry(name=name, dim=dim)
         node.geometry = geometry
-        self.generic_visit(element, geometry)
+        self.generic_visit_children(element, geometry)
 
     def visit_Extent(self, element: _Element, node: vcg.Geometry) -> None:
         X = float(element.get("X", default="1"))
@@ -529,7 +529,7 @@ def visit_SubVolume(self, element: _Element, node: vcg.Geometry) -> None:
             name=name, handle=handle, subvolume_type=subvolume_type, image_pixel_value=image_pixel_value
         )
         node.subvolumes.append(subvolume)
-        self.generic_visit(element, subvolume)
+        self.generic_visit_children(element, subvolume)
 
     def visit_AnalyticExpression(self, element: _Element, node: vcg.SubVolume) -> None:
         expr: str | None = element.text
@@ -554,7 +554,7 @@ def visit_FeatureMapping(self, element: _Element, node: vc.Application) -> None:
             size_exp=size_exp,
         )
         node.compartment_mappings.append(mapping)
-        self.generic_visit(element, mapping)
+        self.generic_visit_children(element, mapping)
 
     def visit_MembraneMapping(self, element: _Element, node: vc.Application) -> None:
         compartment_name: str = element.get("Membrane", default="unknown")
@@ -568,7 +568,7 @@ def visit_MembraneMapping(self, element: _Element, node: vc.Application) -> None
             size_exp=size,
         )
         node.compartment_mappings.append(mapping)
-        self.generic_visit(element, mapping)
+        self.generic_visit_children(element, mapping)
 
     def visit_BoundariesTypes(self, element: _Element, node: vc.CompartmentMapping) -> None:
         switch = {"Flux": vc.BoundaryType.flux, "Value": vc.BoundaryType.value}
@@ -584,7 +584,7 @@ def visit_LocalizedCompoundSpec(self, element: _Element, node: vc.Application) -
         species_name: str = element.get("LocalizedCompoundRef", default="unnamed")
         species_mapping = vc.SpeciesMapping(species_name=species_name)
         node.species_mappings.append(species_mapping)
-        self.generic_visit(element, species_mapping)
+        self.generic_visit_children(element, species_mapping)
 
     def visit_InitialConcentration(self, element: _Element, node: vc.SpeciesMapping) -> None:
         text: str = element.text or "0"
@@ -620,12 +620,19 @@ def visit_Diffusion(self, element: _Element, node: vc.SpeciesMapping) -> None:
             value: str | float = float_or_formula(text)
             node.diff_coef = value
 
+    def visit_Velocity(self, element: _Element, node: vc.SpeciesMapping) -> None:
+        parent = element.getparent()
+        if parent is None or strip_namespace(parent.tag) != "LocalizedCompoundSpec":
+            return
+        values: list[float | str] = [float_or_formula(element.get(dim, default="0.0")) for dim in ["X", "Y", "Z"]]
+        node.velocity_x, node.velocity_y, node.velocity_z = tuple(values)
+
 
 class PrintVisitor(XMLVisitor):
     def visit_root(self, element: _Element, node: vc.VcmlNode) -> None:
         print(f"Visiting root: {element.tag}")
-        self.generic_visit(element, node)
+        self.generic_visit_children(element, node)
 
     def visit_child(self, element: _Element, node: vc.VcmlNode) -> None:
         print(f"Visiting child: {element.tag}")
-        self.generic_visit(element, node)
+        self.generic_visit_children(element, node)

pyvcell/vcml/vcml_writer.py

@@ -381,6 +381,21 @@ def write_species_mapping(self, mapping: SpeciesMapping, parent: _Element) -> No
             parent.append(boundaries_element)
         elif boundary_value_count > 0:
             raise ValueError(f"SpeciesMapping {mapping.species_name} has {boundary_value_count} boundary values")
+        if mapping.velocity_x is not None or mapping.velocity_y is not None or mapping.velocity_z is not None:
+            velocity_element = Element("Velocity")
+            if mapping.velocity_x is not None:
+                str_val_x = str(mapping.velocity_x)
+                if str_val_x != "0.0":
+                    velocity_element.set("X", str_val_x)
+            if mapping.velocity_y is not None:
+                str_val_y = str(mapping.velocity_y)
+                if str_val_y != "0.0":
+                    velocity_element.set("Y", str_val_y)
+            if mapping.velocity_z is not None:
+                str_val_z = str(mapping.velocity_z)
+                if str_val_z != "0.0":
+                    velocity_element.set("Z", str_val_z)
+            parent.append(velocity_element)
 
     @staticmethod
     def _append_text_element(parent: _Element, tag: str, text: str | None) -> None:

tests/conftest.py

@@ -15,6 +15,8 @@
     sbml_spatial_bunny_3d_path,
     sbml_spatial_model_1d_path,
     sbml_spatial_model_3d_path,
+    vcml_sasco_model_path,
+    vcml_sasco_model_with_velocity_path,
 )
 from tests.fixtures.vcell_model_fixtures import (  # noqa: F401
     vcml_field_data_demo_arrays,