sage/period_computation.sage

@@ -1,6 +1,6 @@
 import copy
 import logging
-from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, prod, var
+from sage.all import sage_eval, PolynomialRing, QQ, SR, function, log, matrix, prod, solve, var
 
 load("sage/util.py")
 
@@ -71,6 +71,139 @@ class PFOperator:
         logger.debug("Simplified PFOperator to: %s", simplified.operator_string)
         return simplified
 
+    def change_coordinates(self, *new_coordinates: str) -> "PFOperator":
+        """
+        Rewrites this operator in new coordinates w0, ..., w{n-1} := new_coordinates(z0, ...,
+        z{n-1}), where new_coordinates is given as n expressions in the *old* coordinates
+        z0, ..., z{n-1}; the inverse change of coordinates needed to do this is solved for
+        automatically (the two aren't independent - one determines the other). The new
+        origin w = 0 is the point z0, ..., z{n-1} at which new_coordinates vanishes.
+
+        E.g. pf_operator.change_coordinates("1 - 5*z0", "1/z1") moves the origin to
+        z0 = 1/5, z1 = infinity.
+        """
+        n = self.no_variables
+        if len(new_coordinates) != n:
+            raise ValueError(
+                "Expected %d new coordinate expression(s), got %d." % (n, len(new_coordinates))
+            )
+
+        z_gens = [var('z%d' % i) for i in range(n)]
+        parse_locals = {str(z): z for z in z_gens}
+        try:
+            new_coords = [SR(sage_eval(expr, locals=parse_locals)) for expr in new_coordinates]
+        except Exception as e:
+            raise ValueError("Invalid coordinate expression: %s" % e)
+
+        w_gens = [var('w%d' % i) for i in range(n)]
+        equations = [w_gens[i] == new_coords[i] for i in range(n)]
+        solutions = solve(equations, z_gens, solution_dict=True)
+        if not solutions:
+            raise ValueError("Could not invert the given change of coordinates %s." % (new_coordinates,))
+        inverse_coords = [solutions[0][z] for z in z_gens]
+
+        Y = function('Y')(*w_gens)
+        y = Y.subs({w_gens[i]: new_coords[i] for i in range(n)})
+
+        def apply_theta(expr, i):
+            return z_gens[i] * expr.diff(z_gens[i])
+
+        total = SR(0)
+        for coeff, monomial in self.operator:
+            exponents = monomial.exponents()[0]
+            z_exponents = exponents[:n]
+            theta_exponents = exponents[n:]
+
+            term = y
+            for i, power in enumerate(theta_exponents):
+                for _ in range(power):
+                    term = apply_theta(term, i)
+
+            z_monomial = prod(z_gens[i] ** z_exponents[i] for i in range(n))
+            total += SR(coeff) * z_monomial * term
+
+        total = total.subs({z_gens[i]: inverse_coords[i] for i in range(n)})
+        total = total.simplify_full().expand()
+
+        total_theta_degree = max(sum(monomial.exponents()[0][n:]) for _, monomial in self.operator)
+        theta_ring = PolynomialRing(QQ, n, ['theta%d' % i for i in range(n)])
+        theta_gens = theta_ring.gens()
+
+        # Peel off the coefficient of each derivative order k, highest total order first (so
+        # a coefficient can never still contain an as-yet-unextracted higher derivative as a
+        # factor), and rewrite w^k*D^k Y as a falling-factorial polynomial in theta applied
+        # to Y - the Euler-operator identity theta*(theta-1)*...*(theta-k+1) = w^k*d^k/dw^k,
+        # taken variable by variable. Each term's coefficient is left as a general rational
+        # function of w for now; only once every term has been collected is the whole
+        # operator scaled by their common denominator (see below).
+        remainder = total
+        terms = []
+        for k in sorted(_multi_indices(n, total_theta_degree), key=lambda k: -sum(k)):
+            if sum(k) == 0:
+                dterm = Y
+            else:
+                args = [arg for i in range(n) if k[i] for arg in (w_gens[i], k[i])]
+                dterm = Y.diff(*args)
+
+            coeff = remainder.coefficient(dterm)
+            if coeff == 0:
+                continue
+            remainder -= coeff * dterm
+
+            falling_factorial = theta_ring(1)
+            for i in range(n):
+                for j in range(k[i]):
+                    falling_factorial *= (theta_gens[i] - j)
+
+            w_monomial = prod(w_gens[i] ** k[i] for i in range(n))
+            terms.append(((coeff / w_monomial).simplify_rational(), falling_factorial))
+
+        remainder = remainder.simplify_full()
+        if remainder != 0:
+            raise ValueError("Residual nonzero after change of coordinates: %s" % remainder)
+        if not terms:
+            raise ValueError("Change of coordinates produced the zero operator.")
+
+        w_ring = PolynomialRing(QQ, n, ['w%d' % i for i in range(n)])
+        common_denominator = w_ring(1)
+        for rational_coeff, _ in terms:
+            common_denominator = common_denominator.lcm(w_ring(rational_coeff.denominator()))
+
+        result_terms = {}
+        for rational_coeff, falling_factorial in terms:
+            scaled = w_ring((rational_coeff * SR(common_denominator)).simplify_rational())
+            for w_coeff, w_monomial in scaled:
+                p = w_monomial.exponents()[0]
+                for theta_coeff, theta_monomial in falling_factorial:
+                    q = theta_monomial.exponents()[0]
+                    key = (p, q)
+                    result_terms[key] = result_terms.get(key, QQ(0)) + w_coeff * theta_coeff
+
+        content = gcd([QQ(coeff) for coeff in result_terms.values() if coeff != 0])
+        if content not in (0, 1):
+            result_terms = {key: coeff / content for key, coeff in result_terms.items()}
+
+        def _format_monomial(coeff, p, q):
+            factors = [str(QQ(coeff))]
+            for i in range(n):
+                if p[i]:
+                    factors.append("z%d^%d" % (i, p[i]) if p[i] != 1 else "z%d" % i)
+            for i in range(n):
+                if q[i]:
+                    factors.append("theta%d^%d" % (i, q[i]) if q[i] != 1 else "theta%d" % i)
+            return "*".join(factors)
+
+        terms_string = [
+            _format_monomial(coeff, p, q) for (p, q), coeff in result_terms.items() if coeff != 0
+        ]
+        if not terms_string:
+            raise ValueError("Change of coordinates produced the zero operator.")
+        operator_string = " + ".join(terms_string)
+
+        new_operator = PFOperator(operator_string, no_variables=n)
+        logger.debug("Changed coordinates of PFOperator to: %s", new_operator.operator_string)
+        return new_operator
+
 class PFOperatorAnsatz(PFOperator):
     """
     A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees.
@@ -122,7 +255,17 @@ class PFIdeal:
     def remove_operator(self, pf_operator: PFOperator):
         self.operators.remove(pf_operator)
         logger.info("Removed PFOperator from PFIdeal: %s", pf_operator.operator_string)
-    
+
+    def change_coordinates(self, *new_coordinates: str) -> "PFIdeal":
+        """
+        Moves the origin of the whole system of differential equations to a new point, by
+        returning a new PFIdeal whose operators are each rewritten in new_coordinates; see
+        PFOperator.change_coordinates for what new_coordinates means.
+        """
+        new_operators = [pf_operator.change_coordinates(*new_coordinates) for pf_operator in self.operators]
+        logger.info("Changed coordinates of PFIdeal with %d operator(s) to %s.", len(new_operators), new_coordinates)
+        return PFIdeal(new_operators)
+
     def find_power_series_solution(self, indicials: list, order: int) -> list:
         """
         Finds power series solutions (no logs) to this PFIdeal at given indicial exponents/order.

tests/test_period_computation.py

@@ -125,42 +125,39 @@ def test_find_power_series_solution_handles_rational_indicial():
 
 
 def test_ideal_finds_power_series_solution_and_recovers_first_operator():
-    # Picard--Fuchs ideal for P_{22211}[8]:
-    #   M1 = theta2*(-2*theta1+2*theta2-1) + 2*(theta1-2*theta2-1)*(theta1-2*theta2)*z2
-    #   M2 = theta1^2*(2*(theta1-2*theta2)*z2-theta2) - 16*(2*theta1+1)*(4*theta1+1)*(4*theta1+3)*z1*z2
-    M1 = PFOperator(
-        "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1",
+    # Picard--Fuchs ideal for P_{22211}[8] at the MUM point:
+    #   D1 = Theta_x^2*(Theta_x - 2*Theta_y) - 4*x*(4*Theta_x + 3)*(4*Theta_x + 2)*(4*Theta_x + 1)
+    #   D2 = Theta_y^2 - y*(2*Theta_y - Theta_x + 1)*(2*Theta_y - Theta_x)
+    #
+    # We test moving the ideal to the intersection with the locus of Strong Coupling.
+    D1 = PFOperator(
+        "theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
         no_variables=2,
     )
-    M2 = PFOperator(
-        "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1",
+    D2 = PFOperator(
+        "theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
         no_variables=2,
     )
-    ideal = PFIdeal([M1, M2])
+    moved = PFIdeal([D1, D2]).change_coordinates("z0", "z1 - 1/4")
+    op1, op2 = moved.operators
 
-    solutions = ideal.find_power_series_solution(indicials=[0, QQ(1) / 2], order=6)
+    solutions = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=14)
     assert len(solutions) == 1
 
     period = solutions[0]
-    assert period.apply_operator(M1).coefficients == {}
-    assert period.apply_operator(M2).coefficients == {}
+    assert period.apply_operator(op1).coefficients == {}
+    assert period.apply_operator(op2).coefficients == {}
     # Normalisation convention: the free parameter at the leading (0, 0) coefficient is 1.
     assert period.coefficients[(0, 0)][(0, 0)] == 1
-    assert period.coefficients[(0, 0)][(1, 1)] == -32
+    assert period.coefficients[(0, 0)][(1, 0)] == 0
 
-    # M1 lives entirely within z_degree <= 1, theta_degree <= 2 (a single bare factor of
-    # z2, quadratic in theta), so this is the natural Ansatz level to look for it at.
     recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
     assert len(recovered) == 1
-
-    # recovered[0] should be a scalar multiple of M1 - compare via the coefficient of the
-    # bare theta2 (theta1 in code) monomial, which is nonzero in M1.
-    ratio = (
-        recovered[0].operator.monomial_coefficient(M1.theta_gens[1])
-        / M1.operator.monomial_coefficient(M1.theta_gens[1])
-    )
-    assert ratio != 0
-    assert recovered[0].operator == ratio * M1.operator
+    assert recovered[0].operator == PFOperator(
+        "-2*z1*theta0^2 + 8*z1*theta0*theta1 - 8*z1*theta1^2 + 2*z1*theta0 - 4*z1*theta1 "
+        "+ 2*theta0*theta1 - 2*theta1^2 + theta1",
+        no_variables=2,
+    ).operator
 
 
 def test_single_operator_recovers_itself_from_its_holomorphic_period():
@@ -192,3 +189,15 @@ def test_single_operator_recovers_itself_from_its_holomorphic_period():
     )
     assert ratio != 0
     assert recovered[0].operator == ratio * L.operator
+
+
+def test_change_coordinates_handles_non_monomial_denominator():
+    op = PFOperator("theta0^2 - z0", no_variables=1)
+
+    result = op.change_coordinates("z0/(1-z0)")
+
+    expected = PFOperator(
+        "(1+z0)^3*theta0^2 + z0*(1+z0)^2*theta0 - z0",
+        no_variables=1,
+    )
+    assert result.operator == expected.operator