Repository: matthewwardrop/formulaic

Issue:
Wrong handling of term ordering leads to incorrect differentiation
Hello, differentiation in the latest version sometimes gives incorrect results compared to the earlier version:

```
import formulaic
print(formulaic.__version__)

f = formulaic.Formula("x1 + x2 + x3 + x1:x2 + {x3**2}")
print(f)

for var in ['x1', 'x2', 'x3']:
    print(var, f.differentiate(var, use_sympy=True))
```

Here are the outputs:

```
1.1.0
1 + x1 + x2 + x3 + x3 ** 2 + x1:x2
x1 0 + 1 + 0 + 0 + 0 + x2
x2 0 + 0 + 1 + 0 + 0 + x1
x3 0 + 0 + 0 + 1 + 0 + (2*x3)
```

and

```
1.0.2
1 + x1 + x2 + x3 + x3 ** 2 + x1:x2
x1 root:
    [0, 1, 0, 0, 0, x2]
x2 root:
    [0, 0, 1, 0, 0, x1]
x3 root:
    [0, 0, 0, 1, (2*x3), 0]
```

Seems like the issue is related to the default ordering of the terms, as the correct results can be obtained by explicitly disabling it:

```
f = formulaic.Formula("x1 + x2 + x3 + x1:x2 + {x3**2}")
f.ordering = formulaic.formula.OrderingMethod.NONE
```

Hints from issue comments:
@matthewwardrop: Hi @e-dorigatti ! Thanks for surfacing this! I'll make sure it gets fixed soon.

I didn't realise anyone other than myself had ever shown an interest in the differentiation features. How do you use them in practice?@e-dorigatti: Hi @matthewwardrop thank you for your prompt answer! We are using formulas in our [bofire](https://github.com/experimental-design/bofire) package to express the model for an experimental design, and use the differentiation features to optimize the design with respect to some optimality criteria. In practice we use this to plan experiments :)@matthewwardrop: Nice! Maybe I should upgrade them from experimental status then :).@matthewwardrop: This has been fixed and released as 1.1.1.

Solution patch:
diff --git a/formulaic/formula.py b/formulaic/formula.py
index 3716d4a..b9e02b8 100644
--- a/formulaic/formula.py
+++ b/formulaic/formula.py
@@ -322,20 +322,17 @@ def differentiate(
         use_sympy: bool = False,
     ) -> _SelfType:
         """
-        EXPERIMENTAL: Take the gradient of this formula. When used a linear
-        regression, evaluating a trained model on model matrices generated by
-        this formula is equivalent to estimating the gradient of that fitted
-        form with respect to `wrt`.
+        Take the gradient of this formula with respect to the variables in
+        `wrt`.
+
+        When used a linear regression context, making predictions based on the
+        model matrices generated the differentiated formula is equivalent to
+        estimating the gradient of the fitted model with respect to `wrt`.
 
         Args:
             wrt: The variables with respect to which the gradient should be
                 taken.
             use_sympy: Whether to use sympy to perform symbolic differentiation.
-
-
-        Notes:
-            This method is provisional and may be removed in any future major
-            version.
         """
 
 
@@ -482,27 +479,25 @@ def differentiate(  # pylint: disable=redefined-builtin
         use_sympy: bool = False,
     ) -> SimpleFormula:
         """
-        EXPERIMENTAL: Take the gradient of this formula. When used a linear
-        regression, evaluating a trained model on model matrices generated by
-        this formula is equivalent to estimating the gradient of that fitted
-        form with respect to `wrt`.
+        Take the gradient of this formula with respect to the variables in
+        `wrt`.
+
+        When used a linear regression context, making predictions based on the
+        model matrices generated the differentiated formula is equivalent to
+        estimating the gradient of the fitted model with respect to `wrt`.
 
         Args:
             wrt: The variables with respect to which the gradient should be
                 taken.
             use_sympy: Whether to use sympy to perform symbolic differentiation.
-
-
-        Notes:
-            This method is provisional and may be removed in any future major
-            version.
         """
         return SimpleFormula(
             [
                 differentiate_term(term, wrt, use_sympy=use_sympy)
                 for term in self.__terms
             ],
-            _ordering=self.ordering,
+            # Preserve term ordering even if differentiation modifies degrees/etc.
+            _ordering=OrderingMethod.NONE,
         )
 
     def get_model_matrix(
@@ -784,20 +779,17 @@ def differentiate(  # pylint: disable=redefined-builtin
         use_sympy: bool = False,
     ) -> SimpleFormula:
         """
-        EXPERIMENTAL: Take the gradient of this formula. When used a linear
-        regression, evaluating a trained model on model matrices generated by
-        this formula is equivalent to estimating the gradient of that fitted
-        form with respect to `wrt`.
+        Take the gradient of this formula with respect to the variables in
+        `wrt`.
+
+        When used a linear regression context, making predictions based on the
+        model matrices generated the differentiated formula is equivalent to
+        estimating the gradient of the fitted model with respect to `wrt`.
 
         Args:
             wrt: The variables with respect to which the gradient should be
                 taken.
             use_sympy: Whether to use sympy to perform symbolic differentiation.
-
-
-        Notes:
-            This method is provisional and may be removed in any future major
-            version.
         """
         return cast(
             SimpleFormula,


Test patch:
diff --git a/tests/test_formula.py b/tests/test_formula.py
index 7ce6ad9..a1bf1f5 100644
--- a/tests/test_formula.py
+++ b/tests/test_formula.py
@@ -197,6 +197,9 @@ def test_differentiate(self):
         assert f.differentiate("a") == ["1", "0", "0"]
         assert f.differentiate("c") == ["0", "0", "0"]
 
+        g = Formula("a:b + b:c + c:d - 1")
+        assert g.differentiate("b") == ["a", "c", "0"]  # order preserved
+
     def test_differentiate_with_sympy(self):
         pytest.importorskip("sympy")
         f = Formula("a + b + log(c) - 1")
@@ -208,6 +211,9 @@ def test_differentiate_with_sympy(self):
             "rhs": ["0", "(1/x)"],
         }
 
+        h = Formula("a + {a**2} + b - 1").differentiate("a", use_sympy=True)
+        assert h == ["1", "(2*a)", "0"]  # order preserved
+
     def test_repr(self, formula_expr, formula_exprs):
         assert repr(formula_expr) == "1 + a + b + c + a:b + a:c + b:c + a:b:c"
         assert repr(formula_exprs) == ".lhs:\n    a\n.rhs:\n    1 + b"