julianandJulian Piribauer
6600379681
--------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #20
387 lines
18 KiB
Python
387 lines
18 KiB
Python
import copy
|
|
import logging
|
|
from sage.all import sage_eval, PolynomialRing, QQ, SR, function, prod, solve, var
|
|
|
|
load("sage/util.py")
|
|
|
|
# Logger
|
|
logger = logging.getLogger(__name__)
|
|
logging.basicConfig(level=logging.DEBUG)
|
|
|
|
class PFOperator:
|
|
"""
|
|
A class representing a Picard-Fuchs operator in a given number of variables (moduli).
|
|
Independent of the given order, the variables are assumed to be left of
|
|
the derivatives.
|
|
|
|
The variables extra_locals are needed for symbolic coefficients used for
|
|
operator Ansätze.
|
|
"""
|
|
def _operator_from_string(self, operator_string: str, extra_locals: dict = None):
|
|
z_names = ['z%d' % i for i in range(self.no_variables)]
|
|
theta_names = ['theta%d' % i for i in range(self.no_variables)]
|
|
base_ring = SR if extra_locals else QQ
|
|
self.ring = PolynomialRing(base_ring, z_names + theta_names)
|
|
self.z_gens = self.ring.gens()[:self.no_variables]
|
|
self.theta_gens = self.ring.gens()[self.no_variables:]
|
|
|
|
parse_locals = dict(self.ring.gens_dict())
|
|
if extra_locals:
|
|
parse_locals.update(extra_locals)
|
|
|
|
try:
|
|
operator = self.ring(sage_eval(operator_string, locals=parse_locals))
|
|
return operator
|
|
except Exception as e:
|
|
raise ValueError("Invalid operator string: %s" % e)
|
|
|
|
def __init__(self, operator_string: str, no_variables: int = 1, extra_locals: dict = None):
|
|
self.no_variables = no_variables
|
|
self.operator_string = operator_string
|
|
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
|
|
|
|
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.
|
|
"""
|
|
def __init__(self, no_variables: int, theta_degree: int, z_degree: int):
|
|
self.no_variables = no_variables
|
|
self.theta_degree = theta_degree
|
|
self.z_degree = z_degree
|
|
|
|
z_indices = _multi_indices(self.no_variables, z_degree)
|
|
theta_indices = _multi_indices(self.no_variables, theta_degree)
|
|
|
|
# Every (z_index, theta_index) pair allowed by the degree bounds gets its own
|
|
# fresh unknown, to be solved for once the Ansatz is applied to a period.
|
|
operator_terms = [
|
|
(var("b_" + "_".join(str(x) for x in z_index + theta_index)), (z_index, theta_index))
|
|
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]
|
|
|
|
operator_string = " + ".join(
|
|
"*".join([str(coeff), *_monomial_factors(z_index, "z"), *_monomial_factors(theta_index, "theta")])
|
|
for coeff, (z_index, theta_index) 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)
|
|
|
|
|
|
class PFIdeal:
|
|
"""
|
|
A class representing a Picard-Fuchs ideal, which is a collection of PFOperators.
|
|
"""
|
|
def __init__(self, operators: list):
|
|
if not all(op.no_variables == operators[0].no_variables for op in operators):
|
|
raise ValueError("All operators must have the same number of variables.")
|
|
self.no_variables = operators[0].no_variables
|
|
self.operators = operators
|
|
logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators))
|
|
|
|
def add_operator(self, pf_operator: PFOperator):
|
|
self.operators.append(pf_operator)
|
|
logger.info("Added PFOperator to PFIdeal: %s", pf_operator.operator_string)
|
|
|
|
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.
|
|
|
|
A term c*z^p*theta^q sends a_k*z^k (true exponent k+indicials) to
|
|
c*(k+indicials)^q*a_k*z^(k+p): it shifts index k up by p (p>=0), never down or sideways.
|
|
E.g. z0*theta0 sends a_k*z0^k to (k+rho0)*a_k*z0^(k+1).
|
|
|
|
This is the multivariate Frobenius method: canonical series solutions of a regular
|
|
holonomic D-ideal via its indicial ideal. See M. Saito, B. Sturmfels, N. Takayama,
|
|
"Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and
|
|
Computation in Mathematics vol. 6, Springer, 2000, chs. 2-3.
|
|
"""
|
|
if len(indicials) != self.no_variables:
|
|
raise ValueError(
|
|
"Indicials must have length no_variables=%d, got %d." % (self.no_variables, len(indicials))
|
|
)
|
|
indicials = [QQ(rho) for rho in indicials]
|
|
|
|
operator_terms = [
|
|
[
|
|
(monomial.exponents()[0][:self.no_variables], monomial.exponents()[0][self.no_variables:], coeff)
|
|
for coeff, monomial in pf_operator.operator
|
|
]
|
|
for pf_operator in self.operators
|
|
]
|
|
|
|
def eigenvalue(k, q):
|
|
# theta_i^q_i acts on z_i^(k_i + indicials[i]) as multiplication by
|
|
# (k_i + indicials[i])^q_i; theta^q's combined eigenvalue is the product over i.
|
|
value = QQ(1)
|
|
for i in range(self.no_variables):
|
|
if q[i]:
|
|
value *= (k[i] + indicials[i]) ** q[i]
|
|
return value
|
|
|
|
# Process multi-indices in order of increasing total degree: since every operator
|
|
# monomial has p >= 0 (componentwise), the z^m coefficient of L(y) only ever
|
|
# depends on a_k for k <= m, so by this point every k < m has already been solved.
|
|
z_indices = sorted(_multi_indices(self.no_variables, order), key=sum)
|
|
|
|
# solved[k] holds a_k written as a vector of coefficients over the `dimension`
|
|
# independent solutions found so far (a basis of the solution space up to k).
|
|
dimension = 0
|
|
solved = {}
|
|
|
|
for m in z_indices:
|
|
# For each operator, "coefficient of z^m in L(y) = 0" splits into a diagonal
|
|
# part (the p=0, theta-only monomials, whose unknown is a_m itself) plus a
|
|
# known part contributed by already-solved a_k with k = m - p, p > 0.
|
|
diagonals = []
|
|
known_contributions = []
|
|
for terms in operator_terms:
|
|
diagonal = QQ(0)
|
|
contribution = [QQ(0)] * dimension
|
|
for p, q, coeff in terms:
|
|
k = tuple(m[i] - p[i] for i in range(self.no_variables))
|
|
if any(ki < 0 for ki in k):
|
|
continue # this monomial would need a_k for a negative multi-index k: no such term
|
|
ev = eigenvalue(k, q)
|
|
if ev == 0:
|
|
continue # theta^q kills z_i^(k_i + indicials[i]) here, so this monomial contributes nothing
|
|
if k == m:
|
|
diagonal += coeff * ev # p = 0: coefficient multiplying the still-unknown a_m
|
|
else:
|
|
k_vector = solved[k] # p > 0: a_k is already known, add its contribution
|
|
for i in range(dimension):
|
|
contribution[i] += coeff * ev * k_vector[i]
|
|
diagonals.append(diagonal)
|
|
known_contributions.append(contribution)
|
|
|
|
# An operator with a nonzero diagonal lets us solve a_m = -(known part)/diagonal
|
|
# directly; this is exactly the indicial equation being nonzero at m + indicials.
|
|
active = next((r for r, d in enumerate(diagonals) if d != 0), None)
|
|
|
|
if active is not None:
|
|
value = [-known_contributions[active][i] / diagonals[active] for i in range(dimension)]
|
|
# Every operator's equation at m must independently be satisfied by this
|
|
# same a_m; disagreement means the ideal is inconsistent with these indicials.
|
|
for r, d in enumerate(diagonals):
|
|
if any(d * value[i] + known_contributions[r][i] != 0 for i in range(dimension)):
|
|
raise ValueError(
|
|
"Inconsistent Picard-Fuchs ideal or indicial exponents %s at multidegree %s."
|
|
% (indicials, m)
|
|
)
|
|
solved[m] = value
|
|
else:
|
|
# Resonance: every operator's indicial part vanishes at m + indicials, so a_m
|
|
# cannot be pinned down by this equation. If the already-known lower-degree
|
|
# data still forces a nonzero constraint here, satisfying it would require a
|
|
# log(z)-term solution, which this method (deliberately) does not compute.
|
|
if any(x != 0 for contribution in known_contributions for x in contribution):
|
|
raise ValueError(
|
|
"Resonance at multidegree %s for indicials %s would require a logarithmic "
|
|
"solution, which find_power_series_solution does not compute." % (m, indicials)
|
|
)
|
|
# Otherwise a_m is genuinely free: it starts a new independent solution, so
|
|
# extend every previously solved coefficient with a 0 in this new direction.
|
|
dimension += 1
|
|
for v in solved.values():
|
|
v.append(QQ(0))
|
|
solved[m] = [QQ(0)] * (dimension - 1) + [QQ(1)]
|
|
|
|
log_index = tuple([0] * self.no_variables)
|
|
solutions = [
|
|
Period(
|
|
no_variables=self.no_variables,
|
|
coefficients={log_index: {m: solved[m][i] for m in z_indices}},
|
|
order=order,
|
|
indicials=indicials,
|
|
)
|
|
for i in range(dimension)
|
|
]
|
|
|
|
logger.info(
|
|
"Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order
|
|
)
|
|
return solutions
|