Fix invalid and lossy output of petab_math_str

CHANGELOG.md

@@ -29,6 +29,11 @@
   https://github.com/PEtab-dev/libpetab-python/pull/533,
   https://github.com/PEtab-dev/libpetab-python/pull/535, and
   https://github.com/PEtab-dev/libpetab-python/pull/536)
+* `petab_math_str` now produces valid PEtab math expressions for boolean
+  operators, (in)equality, `min`/`max`, powers with negative bases, and
+  `pi`/`e`, and prints floats without loss of precision. PEtab v2 condition
+  and measurement tables now write expressions in PEtab syntax
+  (by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/545)
 
 **Features**
 * `petab1to2` and `petab_files_1to2` can skip validation via `validate=False`

petab/v1/math/printer.py

@@ -1,9 +1,11 @@
 """A PEtab-compatible sympy string-printer."""
 
 import typing
+from functools import reduce
 from itertools import chain, islice
 
 import sympy as sp
+from sympy.logic.boolalg import Boolean
 from sympy.printing.str import StrPrinter
 
 __all__ = ["PetabStrPrinter", "petab_math_str"]
@@ -35,12 +37,51 @@ def _print_BooleanTrue(self, expr):
     def _print_BooleanFalse(self, expr):
         return "false"
 
+    def _print_Float(self, expr: sp.Float):
+        """Print double-precision floats such that they round-trip exactly."""
+        if expr._prec != 53:
+            return super()._print_Float(expr)
+        value = float(expr)
+        if value.is_integer() and abs(value) < 1e16:
+            return str(int(value))
+        return repr(value)
+
+    def _print_bool_operand(self, expr) -> str:
+        """Print an operand of a boolean or relational operator."""
+        if isinstance(expr, Boolean) and not expr.is_Atom:
+            return f"({self._print(expr)})"
+        return self._print(expr)
+
+    def _print_Relational(self, expr):
+        """Custom printing for comparisons (sympy prints ``Eq(a, b)``)"""
+        lhs = self._print_bool_operand(expr.lhs)
+        rhs = self._print_bool_operand(expr.rhs)
+        return f"{lhs} {expr.rel_op} {rhs}"
+
+    def _print_And(self, expr):
+        """Custom printing for logical AND"""
+        return " && ".join(map(self._print_bool_operand, expr.args))
+
+    def _print_Or(self, expr):
+        """Custom printing for logical OR"""
+        return " || ".join(map(self._print_bool_operand, expr.args))
+
+    def _print_Not(self, expr):
+        """Custom printing for logical NOT"""
+        return f"!{self._print_bool_operand(expr.args[0])}"
+
     def _print_Pow(self, expr: sp.Pow):
         """Custom printing for the power operator"""
         base, exp = expr.as_base_exp()
         str_base = self._print(base)
         str_exp = self._print(exp)
-        if not base.is_Atom:
+        # A negative or non-integer rational base is an Atom, but without
+        # parentheses "-2 ^ x" or "1/2 ^ x" would re-parse as -(2^x) or 1/(2^x)
+        if (
+            not base.is_Atom
+            or (base.is_number and base.is_negative)
+            or (base.is_Rational and not base.is_Integer)
+        ):
             str_base = f"({str_base})"
         # A non-integer Rational exponent (e.g. sqrt -> 1/2) is an Atom but
         # prints as the multi-token "1/2", so without parentheses "x ^ 1/2"
@@ -49,6 +90,14 @@ def _print_Pow(self, expr: sp.Pow):
             str_exp = f"({str_exp})"
         return f"{str_base} ^ {str_exp}"
 
+    def _print_Pi(self, expr):
+        """PEtab has no pi constant; print its double-precision value"""
+        return repr(float(expr))
+
+    def _print_Exp1(self, expr):
+        """Custom printing for Euler's number (no constant in PEtab)"""
+        return "exp(1)"
+
     def _print_Infinity(self, expr):
         """Custom printing for infinity"""
         return "inf"
@@ -70,6 +119,11 @@ def _print_Function(self, expr):
 
     def _print_Piecewise(self, expr):
         """Custom printing for Piecewise function"""
+        if expr.args[-1].cond != sp.true:
+            raise ValueError(
+                "Piecewise without a default value cannot be represented "
+                f"in PEtab: {expr}"
+            )
         # merge the tuples and drop the final `True` condition
         str_args = map(
             self._print,
@@ -78,12 +132,18 @@ def _print_Piecewise(self, expr):
         return f"piecewise({', '.join(str_args)})"
 
     def _print_Min(self, expr):
-        """Custom printing for Min function"""
-        return f"min({', '.join(map(self._print, expr.args))})"
+        """Custom printing for Min function (binary in PEtab)"""
+        return reduce(
+            lambda acc, arg: f"min({acc}, {arg})",
+            map(self._print, expr.args),
+        )
 
     def _print_Max(self, expr):
-        """Custom printing for Max function"""
-        return f"max({', '.join(map(self._print, expr.args))})"
+        """Custom printing for Max function (binary in PEtab)"""
+        return reduce(
+            lambda acc, arg: f"max({acc}, {arg})",
+            map(self._print, expr.args),
+        )
 
 
 def petab_math_str(expr: sp.Basic | sp.Expr | None) -> str:

petab/v2/core.py

@@ -561,7 +561,7 @@ def to_df(self) -> pd.DataFrame:
             record[C.TARGET_VALUE] = (
                 float(record[C.TARGET_VALUE])
                 if record[C.TARGET_VALUE].is_number
-                else str(record[C.TARGET_VALUE])
+                else petab_math_str(record[C.TARGET_VALUE])
             )
         return (
             pd.DataFrame(records)
@@ -840,10 +840,10 @@ def to_df(self) -> pd.DataFrame:
         records = self.model_dump(by_alias=True)["elements"]
         for record in records:
             record[C.OBSERVABLE_PARAMETERS] = C.PARAMETER_SEPARATOR.join(
-                map(str, record[C.OBSERVABLE_PARAMETERS])
+                map(petab_math_str, record[C.OBSERVABLE_PARAMETERS])
             )
             record[C.NOISE_PARAMETERS] = C.PARAMETER_SEPARATOR.join(
-                map(str, record[C.NOISE_PARAMETERS])
+                map(petab_math_str, record[C.NOISE_PARAMETERS])
             )
 
         return pd.DataFrame(records)

tests/v1/math/test_math.py

@@ -47,6 +47,65 @@ def test_printer():
     # re-parses as (a^1)/2 (i.e. sqrt(a) would round-trip to a/2).
     assert petab_math_str(sp.sqrt(a)) == "a ^ (1/2)"
     assert petab_math_str(a ** sp.Rational(2, 3)) == "a ^ (2/3)"
+    # negative or non-integer rational bases must be parenthesized
+    assert petab_math_str(sp.Float(-2) ** a) == "(-2) ^ a"
+    assert petab_math_str(sp.Rational(2, 3) ** a) == "(2/3) ^ a"
+    # boolean and relational operators
+    assert petab_math_str(sp.And(a > 0, b < 1)) == "(a > 0) && (b < 1)"
+    assert petab_math_str(sp.Or(a > 0, b < 1)) == "(a > 0) || (b < 1)"
+    assert petab_math_str(sp.Not(sp.And(a > 0, b < 1))) == (
+        "!((a > 0) && (b < 1))"
+    )
+    assert petab_math_str(sp.Eq(a, 1)) == "a == 1"
+    assert petab_math_str(sp.Ne(a, 1)) == "a != 1"
+    # min/max are binary in PEtab
+    assert petab_math_str(sp.Min(a, b, c)) == "min(min(a, b), c)"
+    assert petab_math_str(sp.Max(a, b, c)) == "max(max(a, b), c)"
+    # no pi / e constants in PEtab
+    assert petab_math_str(sp.pi) == "3.141592653589793"
+    assert petab_math_str(sp.E) == "exp(1)"
+    assert petab_math_str(sp.E * a) == "exp(1)*a"
+    # floats round-trip exactly; integer-valued floats are printed as integers
+    assert petab_math_str(sp.Float(2.0) * a) == "2*a"
+    assert petab_math_str(sp.Float(2 / 3) * a) == "0.6666666666666666*a"
+    assert petab_math_str(sp.Float(1e-300)) == "1e-300"
+    assert petab_math_str(sympify_petab("1.2345678901234567 * a")) == (
+        "1.2345678901234567*a"
+    )
+
+    with pytest.raises(ValueError, match="default"):
+        petab_math_str(sp.Piecewise((1, a > 0)))
+
+
+@pytest.mark.parametrize(
+    "expr_str",
+    [
+        "100 * x",
+        "x ^ 2",
+        "2 * x / 3",
+        "log(x, 10)",
+        "1.2345678901234567 * x",
+        "(-2) ^ x",
+        "(1 / 3) ^ x",
+        "x ^ -0.5 * y",
+        "piecewise(1, x > 0 && y < 2 || x == y, 0)",
+        "piecewise(1, !(x > 0 && y < 2), 0)",
+        "piecewise(1, x != y, 0)",
+        "min(min(x, y), 3) + max(x, max(y, 2))",
+        "2 * arccos(0) + exp(1)",
+    ],
+)
+@pytest.mark.parametrize("evaluate", [True, False])
+def test_print_parse_round_trip(expr_str, evaluate):
+    """Parsing, printing, and parsing again gives the same value."""
+    expr = sympify_petab(expr_str, evaluate=evaluate)
+    printed = petab_math_str(expr)
+    reparsed = sympify_petab(printed, evaluate=evaluate)
+    x, y = sp.symbols("x y", real=True)
+    values = {x: 0.7, y: -1.3}
+    expected = complex(expr.subs(values))
+    actual = complex(reparsed.subs(values))
+    assert np.isclose(actual, expected, rtol=1e-15, atol=0), printed
 
 
 def read_cases():

tests/v2/test_core.py

@@ -461,6 +461,50 @@ def test_condition_table():
     ).free_symbols == {x, y}
 
 
+def test_to_df_math_expressions():
+    """Expressions are written in PEtab syntax and without precision loss.
+
+    Previously, floats were written with only 15 significant digits.
+    """
+    conditions = ConditionTable(
+        [
+            Condition(
+                id="c1",
+                changes=[
+                    Change(target_id="k1", target_value="k2 ^ 2"),
+                    Change(
+                        target_id="k3",
+                        target_value="piecewise(1, k2 > 0 && k2 < 2, 0)",
+                    ),
+                ],
+            )
+        ]
+    )
+    df = conditions.to_df()
+    assert df[TARGET_VALUE].tolist() == [
+        "k2 ^ 2",
+        "piecewise(1, (k2 < 2) && (k2 > 0), 0)",
+    ]
+    assert_frame_equal(ConditionTable.from_df(df).to_df(), df)
+
+    measurements = MeasurementTable(
+        [
+            Measurement(
+                observable_id="obs1",
+                time=1,
+                experiment_id="exp1",
+                measurement=1,
+                observable_parameters="x ^ y;p1",
+                noise_parameters=2 / 3,
+            )
+        ]
+    )
+    df = measurements.to_df()
+    assert df[C.OBSERVABLE_PARAMETERS].tolist() == ["x ^ y;p1"]
+    # full double precision, not str(sp.Float(2 / 3)) == "0.666666666666667"
+    assert df[C.NOISE_PARAMETERS].tolist() == ["0.6666666666666666"]
+
+
 def test_load_remote():
     """Test loading remote files"""
     from jsonschema.exceptions import ValidationError