julian/Calabi-Yau-Period-Geometry
1
0
Fork 0
Code Issues 4 Pull Requests Actions Packages Projects Releases Wiki Activity

Initial setup for period computation #13

Merged
julian merged 11 commits from period-calculations into master 2026-08-17 19:51:33 +02:00
Conversation 0 Commits 11 Files Changed 5
+798 -6
5 changed files with 798 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
README.md
+82 -5
View File
@@ -12,12 +12,21 @@ elliptic_curve_D = ToricPolytopeProjectiveSpace([1, 2, 3], model_name="elliptic_
CY3_quintic = ToricPolytopeProjectiveSpace([1, 1, 1, 1, 1], model_name="quintic")
CY3_bicubic = ToricPolytopeCICY([[3, 3]], model_name="bi-cubic")
CICY3_two_parameter_manual_nef = ToricPolytope([[1,0,0,0,0,0],[0,1,0,0,0,0],[0,0,1,0,0,0],
[0,0,0,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1],
[-1,-1,0,0,0,0],[0,0,-1,-1,-1,-1]],
nef_partition=[[0,1,2,3], [4,5,6,7]])
CICY3_two_parameter_manual_nef = ToricPolytope(
[
[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1],
[-1, -1, 0, 0, 0, 0],
[0, 0, -1, -1, -1, -1],
],
nef_partition=[[0, 1, 2, 3], [4, 5, 6, 7]],
)
CICY5_two_parameter = ToricPolytopeCICY([[6, 1],[0, 2]])
CICY5_two_parameter = ToricPolytopeCICY([[6, 1], [0, 2]])
```
The discriminant factors and topological data, e.g. for the quintic, can then be computed with the methods below.
@@ -72,3 +81,71 @@ DEBUG:__main__:
The result is saved in the folder `data/topdata` as a JSON file — giving a model name helps keeping
track of these outputs.
Note that for Calabi–Yau dimensions larger than four, the additional
## period_computation
`period_computation.sage` provides classes for working with Picard–Fuchs operators and their
period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and
last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given
z- and theta-degree.
A `PFOperator` is parsed from a string in the variables `z0, ..., z<n-1>` and `theta0, ..., theta<n-1>`,
the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard&ndash;Fuchs operator:
```python
L = PFOperator(
"theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0",
no_variables=1,
)
L.simplify().operator_string
```
```term
'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4'
```
An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial
exponents and order.
```python
ideal = PFIdeal([L])
period = ideal.find_power_series_solution(indicials=[0], order=3)[0]
period.period_string
```
```term
'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1'
```
The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for
`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown
coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter
model P_{2,2,2,1,1}[8]:
```python
M1 = PFOperator(
"theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1",
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",
no_variables=2,
)
ideal = PFIdeal([M1, M2])
period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0]
recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
period.period_string
recovered[0].simplify().operator_string
```
```term
'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)'
'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1'
```
`recovered[0]` is, up to scale, `M1` &mdash; recovered purely from `M1`, `M2`'s shared power series
solution.
pyproject.toml
+3
View File
@@ -15,3 +15,6 @@ ignore = [
# test_smoke.py star-imports sage.all and loads a .sage file, so ruff can't see where its names come from.
"tests/test_topdata_and_disc.py" = ["F403", "F405"]
# Same star-import + load() pattern as above.
"tests/test_period_computation.py" = ["F403", "F405"]
sage/period_computation.sage
+503
View File
@@ -0,0 +1,503 @@
import copy
import logging
from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, prod, 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
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 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
class Period:
"""
A class representing a period as a formal power series in z-variables and their logs.
The coefficients are stored in a dictionary of dictionaries, where the first key is the multi-index
of the logarithmic part and the second key is the multi-index of the z-variables.
So, for example, the coefficient of
log(z0)^2 * log(z1) * z0^3 * z1^2 would be stored as coefficients[(2, 1)][(3, 2)]
"""
def _initialise_ring(self):
z_names = ['z%d' % i for i in range(self.no_variables)]
log_names = ['L%d' % i for i in range(self.no_variables)]
ring = PolynomialRing(QQ, z_names + log_names)
z_gens = ring.gens()[:self.no_variables]
log_gens = ring.gens()[self.no_variables:]
return ring, z_gens, log_gens # type: (PolynomialRing, list, list)
def _period_from_string(self, period_string: str) -> dict:
ring, z_gens, log_gens = self._initialise_ring()
log_of_z = dict(zip(z_gens, log_gens))
def log(zi):
try:
return log_of_z[zi]
except (KeyError, TypeError):
raise ValueError(
"log(...) may only be applied to one of the z-variables z0, ..., z%d"
% (self.no_variables - 1)
)
parse_locals = dict(ring.gens_dict())
parse_locals['log'] = log
try:
expression = ring(sage_eval(period_string, locals=parse_locals))
except Exception as e:
raise ValueError("Invalid period string: %s" % e)
coefficients = {}
for coeff, monomial in expression:
exponents = monomial.exponents()[0]
z_index = tuple(exponents[:self.no_variables])
log_index = tuple(exponents[self.no_variables:])
coefficients.setdefault(log_index, {})[z_index] = coeff
return coefficients
def _period_to_string(self) -> str:
ring, z_gens, log_gens = self._initialise_ring()
expression = ring(0)
for log_index, z_dict in self.coefficients.items():
for z_index, coeff in z_dict.items():
monomial = coeff
for i in range(self.no_variables):
monomial *= (log_gens[i] ** log_index[i]) * (z_gens[i] ** z_index[i])
expression += monomial
log_substitutions = {log_gens[i]: log(SR(z_gens[i])) for i in range(self.no_variables)}
result = SR(expression).subs(log_substitutions)
# A nonzero indicial ρ_i means the coefficients above are for z_i^k, but the
# actual solution is z_i^(ρ_i + k); make that explicit in the printed form.
indicial_prefactor = prod(
(SR(z_gens[i]) ** self.indicials[i] for i in range(self.no_variables) if self.indicials[i] != 0),
SR(1),
)
if indicial_prefactor != 1:
result *= indicial_prefactor
return str(result)
def _max_z_degree(self) -> int:
# Highest total z-degree (sum of the z-multi-index) among all coefficients.
if not self.coefficients:
return 0
return max(
sum(z_index)
for z_dict in self.coefficients.values()
for z_index in z_dict
)
def _truncate_coefficients(self, order: int) -> dict:
# Drop every (log_index, z_index) entry whose total z-degree exceeds order,
# removing it from the dictionary rather than merely zeroing it out.
truncated = {}
for log_index, z_dict in self.coefficients.items():
kept = {z_index: coeff for z_index, coeff in z_dict.items() if sum(z_index) <= order}
if kept:
truncated[log_index] = kept
return truncated
def _apply_theta(self, coefficients: dict, index: int) -> dict:
"""
Applies the logarithmic derivative theta_i = z_i * d/dz_i once to a coefficients dict of the same
shape as self.coefficients. It uses the product rule
theta_i(z^a log(z)^k) = a_i * z^a log(z)^k + k_i * z^a log(z)^(k - e_i).
z_index entries are offsets from self.indicials: a stored z_index of a really means
z^(a + self.indicials[index]), so theta_i's eigenvalue is a_i + self.indicials[index]
rather than the bare a_i (self.indicials is all-zero unless explicitly given, in which
case this reduces to the ordinary power-series rule).
"""
result = {}
for log_index, z_dict in coefficients.items():
for z_index, coeff in z_dict.items():
a_i = z_index[index] + self.indicials[index]
if a_i != 0:
inner = result.setdefault(log_index, {})
inner[z_index] = inner.get(z_index, 0) + a_i * coeff
k_i = log_index[index]
if k_i != 0:
lowered_log_index = log_index[:index] + (k_i - 1,) + log_index[index + 1:]
inner = result.setdefault(lowered_log_index, {})
inner[z_index] = inner.get(z_index, 0) + k_i * coeff
return result
def apply_operator(self, pf_operator: PFOperator) -> "Period":
"""
Applies a PFOperator to this period and returns the result as a new Period,
truncated to self.order. Each monomial of the operator is normalised as
coeff * z^p * theta^q (see PFOperator's docstring), so theta^q is applied
to the period first and the result is then multiplied by coeff * z^p.
"""
if pf_operator.no_variables != self.no_variables:
raise ValueError(
"Variable count mismatch: period has %d variable(s), operator has %d."
% (self.no_variables, pf_operator.no_variables)
)
result_coefficients = {}
for monomial_coeff, monomial in pf_operator.operator:
exponents = monomial.exponents()[0]
z_exponents = exponents[:self.no_variables]
theta_exponents = exponents[self.no_variables:]
term = self.coefficients
for i, power in enumerate(theta_exponents):
for _ in range(power):
term = self._apply_theta(term, i)
for log_index, z_dict in term.items():
inner = result_coefficients.setdefault(log_index, {})
for z_index, coeff in z_dict.items():
shifted_z_index = tuple(z_index[i] + z_exponents[i] for i in range(self.no_variables))
inner[shifted_z_index] = inner.get(shifted_z_index, 0) + monomial_coeff * coeff
result_coefficients = {
log_index: {z_index: c for z_index, c in z_dict.items() if c != 0}
for log_index, z_dict in result_coefficients.items()
}
result_coefficients = {log_index: z_dict for log_index, z_dict in result_coefficients.items() if z_dict}
return Period(
no_variables=self.no_variables,
coefficients=result_coefficients,
order=self.order,
indicials=self.indicials,
)
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,
indicials: list = None,
):
self.no_variables = no_variables
# indicials[i] is the (rational) Frobenius exponent ρ_i of z_i: a stored
# z-index of a really represents z^(a + indicials[i]).
self.indicials = list(indicials) if indicials is not None else [0] * no_variables
# Use coefficients or period_string to initialize the period
if period_string is not None:
if coefficients is not None:
raise ValueError("Provide either coefficients or period_string, not both.")
self.coefficients = self._period_from_string(period_string)
self.period_string = period_string
logger.debug("Using string for initialisation.")
elif coefficients is None:
self.coefficients = {}
else:
self.coefficients = coefficients
self.period_string = self._period_to_string()
logger.debug("Using coefficients for initialisation.")
if order is None:
self.order = self._max_z_degree()
logger.debug("Order not provided, using maximum z-degree: %d", self.order)
else:
self.coefficients = self._truncate_coefficients(order)
self.order = order
self.period_string = self._period_to_string()
logger.debug("Truncated coefficients to order %d.", self.order)
logger.info("Initialised Period in %d variables at order %d.", self.no_variables, self.order)
class PeriodAnsatz(Period):
"""
A class for period Ansätze, characterised by number of variables and their z- and log-multi-degrees.
"""
def __init__(self, no_variables: int, z_degree: int, log_degree: int, indicials: list = None):
self.no_variables = no_variables
self.z_degree = z_degree
self.log_degree = log_degree
log_indices = _multi_indices(self.no_variables, log_degree)
z_indices = _multi_indices(self.no_variables, z_degree)
# Every (log_index, z_index) pair allowed by the degree bounds gets its own
# fresh unknown, to be solved for once a PFOperator is applied to the Ansatz.
coefficients = {
log_index: {
z_index: var("a_" + "_".join(str(x) for x in log_index + z_index))
for z_index in z_indices
}
for log_index in log_indices
}
super().__init__(no_variables=no_variables, coefficients=coefficients, order=z_degree, indicials=indicials)
self.expansion_coefficients = self.coefficients
logger.info(
"Initialised PeriodAnsatz with %d unknown coefficient(s).",
sum(len(z_dict) for z_dict in self.expansion_coefficients.values()),
)
sage/util.py
+15
View File
@@ -0,0 +1,15 @@
"""
Utility functions for the Calabi-Yau period geometry project.
"""
def _multi_indices(no_variables, total_degree: int) -> list:
# All no_variables-tuples of non-negative integers whose entries sum to at most total_degree.
def helper(remaining_vars, remaining_degree):
if remaining_vars == 0:
yield ()
return
for first in range(remaining_degree + 1):
for rest in helper(remaining_vars - 1, remaining_degree - first):
yield (first,) + rest
return list(helper(no_variables, total_degree))
tests/test_period_computation.py
+194
View File
@@ -0,0 +1,194 @@
import os
from sage.all import * # noqa: F401
load(os.path.join(os.path.dirname(__file__), "..", "sage", "period_computation.sage"))
def test_apply_operator_univariate():
# theta^2 - z applied to log(z) should give -z*log(z), since theta(log z) = 1
# and theta^2(log z) = theta(1) = 0.
period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5)
op = PFOperator("theta0^2 - z0", no_variables=1)
result = period.apply_operator(op)
assert result.coefficients == {(1,): {(1,): -1}}
def test_apply_operator_mixes_variables():
# z0*theta1 applied to z0*log(z1): theta1 strips log(z1) down to a bare 1
# (leaving z0 untouched), then multiplying by z0 gives z0^2.
period = Period(no_variables=2, coefficients={(0, 1): {(1, 0): 1}}, order=5)
op = PFOperator("z0*theta1", no_variables=2)
result = period.apply_operator(op)
assert result.coefficients == {(0, 0): {(2, 0): 1}}
def test_apply_operator_truncates_to_order():
# Multiplying by z0^2 pushes some terms above the period's order, so they
# should be dropped rather than kept with a nonzero coefficient.
period = Period(no_variables=1, coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1}}, order=2)
op = PFOperator("z0^2", no_variables=1)
result = period.apply_operator(op)
assert result.coefficients == {(0,): {(2,): 1}}
assert result.order == 2
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"
def test_find_power_series_solution_recovers_quintic_period():
# The quintic's Picard-Fuchs operator has indicial equation theta0^4 = 0 at z0 = 0
# (a quadruple root at 0), so its holomorphic (non-logarithmic) power series solution
# is found at indicial 0. Up to normalisation this is the classic quintic period
# 1 + 120*z0 + 113400*z0^2 + ...
op = PFOperator(
"theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)",
no_variables=1,
)
ideal = PFIdeal([op])
solutions = ideal.find_power_series_solution(indicials=[0], order=3)
assert len(solutions) == 1
assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000}}
def test_find_power_series_solution_handles_rational_indicial():
# theta0 - 1/2 kills z0^(1/2 + k) only when k = 0, since its eigenvalue is 1/2 + k;
# so at indicial 1/2 there is exactly one solution (a constant multiple of sqrt(z0)),
# while at indicial 0 no power series solution exists at all.
op = PFOperator("theta0 - 1/2", no_variables=1)
ideal = PFIdeal([op])
solutions = ideal.find_power_series_solution(indicials=[1/2], order=3)
assert len(solutions) == 1
assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 0, (2,): 0, (3,): 0}}
assert solutions[0].period_string == "sqrt(z0)"
assert ideal.find_power_series_solution(indicials=[0], order=3) == []
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",
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",
no_variables=2,
)
ideal = PFIdeal([M1, M2])
solutions = ideal.find_power_series_solution(indicials=[0, QQ(1) / 2], order=6)
assert len(solutions) == 1
period = solutions[0]
assert period.apply_operator(M1).coefficients == {}
assert period.apply_operator(M2).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
# 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
def test_single_operator_recovers_itself_from_its_holomorphic_period():
# Ensuring methods work for operator 4.2.1:
L = PFOperator(
"theta0^4 - 4*z0*(2*theta0 + 1)^2*(7*theta0^2 + 7*theta0 + 2) "
"- 128*z0^2*(2*theta0 + 1)^2*(2*theta0 + 3)^2",
no_variables=1,
)
ideal = PFIdeal([L])
solutions = ideal.find_power_series_solution(indicials=[0], order=15)
assert len(solutions) == 1
period = solutions[0]
assert period.apply_operator(L).coefficients == {}
assert period.coefficients[(0,)][(0,)] == 1
assert period.coefficients[(0,)][(1,)] == 8
assert period.coefficients[(0,)][(2,)] == 360
# L lives entirely within z_degree <= 2, theta_degree <= 4, its stated level.
recovered = period.find_annihilating_operators(z_degree=2, theta_degree=4)
assert len(recovered) == 1
# recovered[0] should be a scalar multiple of L - compare via the coefficient of the
# bare theta0^4 monomial, which is 1 in L.
ratio = recovered[0].operator.monomial_coefficient(L.theta_gens[0] ** 4) / L.operator.monomial_coefficient(
L.theta_gens[0] ** 4
)
assert ratio != 0
assert recovered[0].operator == ratio * L.operator