Finding operators and simplifying them

sage/period_computation.sage

@@ -1,5 +1,6 @@
+import copy
 import logging
-from sage.all import sage_eval, PolynomialRing, QQ, SR, log, var
+from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, var
 
 load("sage/util.py")
 
@@ -29,7 +30,7 @@ class PFOperator:
             parse_locals.update(extra_locals)
 
         try:
-            operator = sage_eval(operator_string, locals=parse_locals)
+            operator = self.ring(sage_eval(operator_string, locals=parse_locals))
             return operator
         except Exception as e:
             raise ValueError("Invalid operator string: %s" % e)
@@ -40,6 +41,36 @@ class PFOperator:
         self.operator = self._operator_from_string(operator_string, extra_locals=extra_locals)
         logger.info("Initialised PFOperator: %s", self.operator)
 
+    def simplify(self) -> "PFOperator":
+        """
+        Groups the operator's terms by z-monomial and factorises the theta-polynomial
+        multiplying each z-monomial, e.g. turning a computed quintic operator into the
+        well-known theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4).
+        """
+        z_monomial_theta_parts = {}
+        for coeff, monomial in self.operator:
+            exponents = monomial.exponents()[0]
+            z_exponents = exponents[:self.no_variables]
+            theta_exponents = exponents[self.no_variables:]
+
+            theta_monomial = SR(1)
+            for gen, exp in zip(self.theta_gens, theta_exponents):
+                theta_monomial *= SR(gen) ** exp
+
+            z_monomial_theta_parts[z_exponents] = z_monomial_theta_parts.get(z_exponents, SR(0)) + SR(coeff) * theta_monomial
+
+        simplified_expr = SR(0)
+        for z_exponents, theta_part in z_monomial_theta_parts.items():
+            z_monomial = SR(1)
+            for gen, exp in zip(self.z_gens, z_exponents):
+                z_monomial *= SR(gen) ** exp
+            simplified_expr += theta_part.factor() * z_monomial
+
+        simplified = copy.copy(self)
+        simplified.operator_string = str(simplified_expr)
+        logger.debug("Simplified PFOperator to: %s", simplified.operator_string)
+        return simplified
+
 class PFOperatorAnsatz(PFOperator):
     """
     A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees.
@@ -59,6 +90,7 @@ class PFOperatorAnsatz(PFOperator):
             for z_index in z_indices
             for theta_index in theta_indices
         ]
+        self.unknowns = [coeff for coeff, _ in operator_terms]
 
         def _monomial_factors(index, name):
             return [f"{name}{i}^{exp}" for i, exp in enumerate(index) if exp > 0]
@@ -68,7 +100,7 @@ class PFOperatorAnsatz(PFOperator):
             for coeff, (z_index, theta_index) in operator_terms
         )
 
-        extra_locals = {str(coeff): coeff for coeff, _ in operator_terms}
+        extra_locals = {str(coeff): coeff for coeff in self.unknowns}
         super().__init__(operator_string=operator_string, no_variables=no_variables, extra_locals=extra_locals)
 
 
@@ -220,6 +252,31 @@ class Period:
             order=self.order,
         )
 
+    def find_annihilating_operators(self, z_degree: int, theta_degree: int) -> list:
+        """
+        Finds operators annihilating this period among PFOperator Ansätze of the given
+        z- and theta-degree.
+        """
+        ansatz = PFOperatorAnsatz(self.no_variables, theta_degree, z_degree)
+        unknowns = ansatz.unknowns
+
+        equations = [
+            coeff
+            for z_dict in self.apply_operator(ansatz).coefficients.values()
+            for coeff in z_dict.values()
+        ]
+
+        coefficient_rows = [[SR(equation).coefficient(b) for b in unknowns] for equation in equations]
+        kernel_basis = matrix(QQ, coefficient_rows, ncols=len(unknowns)).right_kernel().basis()
+
+        operators = []
+        for basis_vector in kernel_basis:
+            solution = {unknowns[j]: basis_vector[j] for j in range(len(unknowns))}
+            solved_operator = ansatz.operator.map_coefficients(lambda c: SR(c).subs(solution))
+            operators.append(PFOperator(str(solved_operator), no_variables=self.no_variables))
+
+        return operators
+
     def __init__(self, no_variables: int = 1, coefficients: dict = None, period_string: str = None, order: int = None):
         self.no_variables = no_variables
 
@@ -278,10 +335,3 @@ class PeriodAnsatz(Period):
             "Initialised PeriodAnsatz with %d unknown coefficient(s).",
             sum(len(z_dict) for z_dict in self.expansion_coefficients.values()),
         )
-
-
-"""
-To do: 
-    - add method to Period that finds operators that annihilate it (using operator ansatz class)
-    - add class PFideal containing a list of PFOperators
-"""

sage/util.py

@@ -13,4 +13,3 @@ def _multi_indices(no_variables, total_degree: int) -> list:
                 yield (first,) + rest
 
     return list(helper(no_variables, total_degree))
-

tests/test_period_computation.py

@@ -46,3 +46,56 @@ def test_apply_operator_rejects_variable_count_mismatch():
 
     with pytest.raises(ValueError):
         period.apply_operator(op)
+
+
+def test_find_annihilating_operators_recovers_theta_squared():
+    # theta0^2 annihilates log(z0): theta0(log z0) = 1, theta0^2(log z0) = theta0(1) = 0.
+    # Among degree-(0, 2) Ansatze this should be the only solution, up to scaling.
+    period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5)
+
+    ops = period.find_annihilating_operators(z_degree=0, theta_degree=2)
+
+    # A one-dimensional solution space means exactly one basis operator.
+    assert len(ops) == 1
+    assert period.apply_operator(ops[0]).coefficients == {}
+    assert str(ops[0].operator) == "theta0^2"
+
+
+def test_find_annihilating_operators_recovers_geometric_series_operator():
+    # sum_{k=0}^{4} z0^k is annihilated (up to truncation order) by
+    # (1 - z0)*theta0 - z0, i.e. theta0 - z0*theta0 - z0.
+    period = Period(
+        no_variables=1,
+        coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1, (3,): 1, (4,): 1}},
+        order=4,
+    )
+
+    ops = period.find_annihilating_operators(z_degree=1, theta_degree=1)
+
+    assert len(ops) == 1
+    assert period.apply_operator(ops[0]).coefficients == {}
+
+
+def test_find_annihilating_operators_returns_empty_list_when_no_solution_exists():
+    # No degree-0 (constant) operator other than the zero operator can annihilate a
+    # nonzero constant period, so the linear system's only solution is trivial.
+    period = Period(no_variables=1, coefficients={(0,): {(0,): 1}}, order=0)
+
+    ops = period.find_annihilating_operators(z_degree=0, theta_degree=0)
+
+    assert ops == []
+
+
+def test_simplify_factorises_theta_polynomial_per_z_monomial():
+    # theta0^4 - 5*z0*(5*theta0+1)*(5*theta0+2)*(5*theta0+3)*(5*theta0+4), expanded, is the
+    # quintic's Picard-Fuchs operator. simplify() should recover the factorised form: the
+    # z0^0 part (theta0^4) has no theta-factor to pull out, while the z0^1 part factorises
+    # into the four linear pieces.
+    expanded = "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0"
+    op = PFOperator(expanded, no_variables=1)
+
+    simplified = op.simplify()
+
+    # The underlying (expanded) operator is unchanged - only the display string differs.
+    assert simplified.operator == op.operator
+    assert simplified.operator_string == "-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4"