Compare commits
10
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Julian Piribauer
|
d7ed05771b | ||
|
Julian Piribauer
|
aadc1f37ea | ||
|
Julian Piribauer
|
d17bdc8949 | ||
|
Julian Piribauer
|
3cc7bc0dda | ||
|
Julian Piribauer
|
53a1eb6ca7 | ||
|
Julian Piribauer
|
80072b291f | ||
|
julian
|
d649d5047f | ||
|
julian
|
3e6660a414 | ||
|
Julian Piribauer
|
1f930b01fc | ||
|
Julian Piribauer
|
ece770ab48 |
@@ -1,191 +1,2 @@
|
||||
# Calabi-Yau-Period-Geometry
|
||||
|
||||
This project collects code used for analysing Calabi–Yau families.
|
||||
It allows for computation of discriminant loci and topological data of manifolds defined as hypersurfaces or complete intersections in toric ambient spaces.
|
||||
|
||||
## toric_topdata
|
||||
|
||||
Supported initialisations are represented by the following examples.
|
||||
|
||||
```python
|
||||
elliptic_curve_D = ToricPolytopeProjectiveSpace([3, 2, 1], model_name="elliptic_curve_D")
|
||||
|
||||
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]],
|
||||
)
|
||||
|
||||
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.
|
||||
```python
|
||||
CY3_quintic.disc()
|
||||
|
||||
CY3_quintic.topdata()
|
||||
```
|
||||
|
||||
We list the output of above two lines.
|
||||
|
||||
```term
|
||||
INFO:__main__:Discriminant factors: [[z1 + 1/3125, 0]]
|
||||
INFO:__main__:
|
||||
|
||||
--- Polytope and GLSM table ------------------
|
||||
1|1|1|1|1|| 1
|
||||
|
||||
-1|1|0|0|0|| 0
|
||||
-1|0|1|0|0|| 0
|
||||
-1|0|0|1|0|| 0
|
||||
-1|0|0|0|1|| 0
|
||||
--------------
|
||||
1,1,1,1,1; -5
|
||||
|
||||
--- L-vectors (GLSM charges) -----------------
|
||||
[[1, 1, 1, 1, 1, -5]]
|
||||
|
||||
--- Intersection numbers CY ------------------
|
||||
{(0, 0, 0): 5}
|
||||
|
||||
--- Intersection ring ------------------------
|
||||
5*t0^3
|
||||
|
||||
--- Intersection ring no multiplicities ------
|
||||
5*t0^3
|
||||
|
||||
--- Chern polynomials ------------------------
|
||||
[[0], [10*t0^3], [-40*t0^3]]
|
||||
|
||||
--- Integrated Chern classes -----------------
|
||||
[[0], [50], [-200]]
|
||||
DEBUG:__main__:
|
||||
|
||||
--- Kähler cone generators (ambient space) ---
|
||||
['[z4]']
|
||||
|
||||
--- Intersection numbers ambient space -------
|
||||
{(0, 0, 0, 0): 1}
|
||||
```
|
||||
|
||||
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
|
||||
|
||||
## Operators and periods
|
||||
|
||||
`operators.sage` and `periods.sage` provide 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. The operator classes live in `operators.sage`, the period classes in `periods.sage`;
|
||||
since the two depend on each other, load `periods.sage`, which also loads `operators.sage`.
|
||||
|
||||
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–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` — recovered purely from `M1`, `M2`'s shared power series
|
||||
solution.
|
||||
|
||||
### Frobenius bases
|
||||
|
||||
A `PFIdeal` also computes its complete Frobenius basis — all solutions, logarithmic ones included —
|
||||
at the origin of its current coordinates; for another point, move it to the origin with `change_coordinates`
|
||||
first. The exponents are the `indicials`, the roots of the system's indicial ideal, which is completed by linear
|
||||
algebra where the operators' z-free parts alone don't determine it. The optional `hodge_numbers` — the
|
||||
numbers of solutions of top log-degree 0, 1, ..., n at a MUM point, e.g. `[1, 2, 2, 1]` for a two-parameter
|
||||
threefold — certify that completion through their sum, the rank, and bound the powers of logarithms. For
|
||||
P_{2,2,2,1,1}[8] at its MUM point and at the point (z0, z1) = (0, 1/4) on its strong coupling divisor:
|
||||
|
||||
```python
|
||||
D1 = PFOperator(
|
||||
"theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
|
||||
no_variables=2,
|
||||
)
|
||||
D2 = PFOperator("theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)", no_variables=2)
|
||||
ideal = PFIdeal([D1, D2], hodge_numbers=[1, 2, 2, 1])
|
||||
strong_coupling = ideal.change_coordinates("z0", "z1 - 1/4")
|
||||
|
||||
ideal.indicials()
|
||||
strong_coupling.indicials()
|
||||
basis = strong_coupling.frobenius_basis(order=2)
|
||||
basis[5].period_string
|
||||
```
|
||||
|
||||
```term
|
||||
[((0, 0), 6)]
|
||||
[((0, 0), 4), ((0, 1/2), 2)]
|
||||
INFO:__main__:Completed the indicial ideal at z = 0 at multiplier degree 2 with 1 extra operator(s): colength 6.
|
||||
INFO:__main__:Exponent class of (0, 0): 4 solution(s) with leading terms ['1', 'log(z0)', 'log(z0)^2', 'log(z0)^3'].
|
||||
INFO:__main__:Exponent class of (0, 1/2): 2 solution(s) with leading terms ['sqrt(z1)', 'sqrt(z1)*(log(z0) + log(z1))'].
|
||||
INFO:__main__:Found a Frobenius basis of 6 solution(s) at order 2.
|
||||
'-1/225*(7200*z0*z1*log(z0) - 720*z1^2*log(z0) + 7200*z0*z1*log(z1) - 720*z1^2*log(z1) + 33600*z0*z1 - 732*z1^2 + 300*z1*log(z0) + 300*z1*log(z1) + 250*z1 - 225*log(z0) - 225*log(z1))*sqrt(z1)'
|
||||
```
|
||||
|
||||
Within each exponent class the basis is canonical: every solution has its own leading term with coefficient 1,
|
||||
at which the class' other solutions vanish. Coefficients are those of plain powers of log(z_i), without factors
|
||||
of 2πi.
|
||||
+1
-4
@@ -14,7 +14,4 @@ ignore = [
|
||||
"*.sage.py" = ["F403", "F405", "F821", "E741", "E402", "E702", "I001", "W291", "W293", "W292"]
|
||||
|
||||
# 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"]
|
||||
"tests/test_smoke.py" = ["F403", "F405"]
|
||||
|
||||
@@ -1,811 +0,0 @@
|
||||
"""
|
||||
Picard-Fuchs operators and ideals in the logarithmic derivatives theta_i = z_i * d/dz_i.
|
||||
|
||||
Contents
|
||||
--------
|
||||
PFOperator
|
||||
__init__ Construct from a string in z- and theta-variables
|
||||
simplify Group by z-monomials and factorise theta-polynomials
|
||||
change_coordinates Rewrite the operator in new coordinates
|
||||
_operator_from_string Parse a string into the operator polynomial
|
||||
_theta_parts Theta-polynomial multiplying each z-monomial
|
||||
|
||||
PFOperatorAnsatz(PFOperator)
|
||||
__init__ Ansatz of given theta- and z-multi-degrees
|
||||
|
||||
PFIdeal
|
||||
__init__ Construct from a list of PFOperators and optional Hodge numbers
|
||||
add_operator Append a PFOperator
|
||||
remove_operator Remove a PFOperator
|
||||
change_coordinates Rewrite all operators in new coordinates
|
||||
indicial_ideal Indicial ideal at z = 0, completed by linear algebra
|
||||
indicials Indicials (exponents) at z = 0 with multiplicities
|
||||
frobenius_basis Basis of all solutions at z = 0, logarithmic ones included
|
||||
find_power_series_solution Power series solutions at given indicials/order
|
||||
_operator_theta_parts Theta-parts of all operators in a common theta-ring
|
||||
_complete_indicial_ideal Indicial ideal and the extra operators needed to see it
|
||||
_solve_exponent_class Solutions with given indicials, up to a degree in log(z)
|
||||
|
||||
Helpers
|
||||
_z_order_key Graded order on z-monomials (initial forms, solving order)
|
||||
_left_multiply Normal-ordered product z^a * theta^c * operator
|
||||
_indicial_roots Roots of an indicial ideal with their multiplicities
|
||||
_leading_term Lowest term of a solution, for logging
|
||||
"""
|
||||
import copy
|
||||
import logging
|
||||
from sage.all import (
|
||||
sage_eval, PolynomialRing, QQ, SR, ZZ, binomial, block_matrix, floor, function, log, matrix, 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
|
||||
|
||||
def _theta_parts(self, theta_ring) -> dict:
|
||||
"""
|
||||
Splits this operator as sum_p z^p * P_p(theta) and returns {p: P_p}, each P_p an element of
|
||||
theta_ring = QQ[theta0, ..., theta{n-1}]; z-monomials whose theta-polynomial vanishes are dropped.
|
||||
theta_ring has to be built as PolynomialRing(QQ, n, names): without the explicit n, a single
|
||||
variable gives a univariate ring, whose exponents() and monomial() work differently.
|
||||
"""
|
||||
if self.ring.base_ring() != QQ:
|
||||
raise ValueError("Only operators with rational coefficients can be split into theta-parts.")
|
||||
theta_gens = theta_ring.gens()
|
||||
parts = {}
|
||||
for coeff, monomial in self.operator:
|
||||
exponents = monomial.exponents()[0]
|
||||
z_exponents = tuple(exponents[:self.no_variables])
|
||||
theta_monomial = prod(theta_gens[i] ** exponents[self.no_variables + i] for i in range(self.no_variables))
|
||||
parts[z_exponents] = parts.get(z_exponents, theta_ring(0)) + coeff * theta_monomial
|
||||
return {p: P for p, P in parts.items() if P != 0}
|
||||
|
||||
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.
|
||||
|
||||
hodge_numbers optionally lists the numbers of solutions of top log-degree 0, 1, ..., n at a point of
|
||||
maximal unipotent monodromy, i.e. h^{n,0}, ..., h^{0,n} of the (horizontal part of the) middle
|
||||
cohomology - e.g. [1, 2, 2, 1] for a two-parameter Calabi-Yau threefold. Their sum is the holonomic
|
||||
rank, the number of solutions at every point, and n bounds the degree in log(z) of every solution at
|
||||
every point (monodromy theorem: N^(n+1) = 0 for a weight-n variation of Hodge structure). Without them,
|
||||
indicials and frobenius_basis still work, just without the certificates and checks these allow.
|
||||
"""
|
||||
def __init__(self, operators: list, hodge_numbers: list = None):
|
||||
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.")
|
||||
if hodge_numbers is not None and (not hodge_numbers or any(h not in ZZ or h < 1 for h in hodge_numbers)):
|
||||
raise ValueError("hodge_numbers must be a non-empty list of positive integers, got %s." % (hodge_numbers,))
|
||||
self.no_variables = operators[0].no_variables
|
||||
self.operators = operators
|
||||
self.hodge_numbers = list(hodge_numbers) if hodge_numbers is not None else None
|
||||
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, hodge_numbers=self.hodge_numbers)
|
||||
|
||||
def indicial_ideal(self, max_multiplier_degree: int = 4):
|
||||
"""
|
||||
Returns the indicial ideal of this PFIdeal at z = 0, an ideal of QQ[theta0, ..., theta{n-1}] whose
|
||||
roots are the indicials (exponents) of the solutions there; see _complete_indicial_ideal.
|
||||
"""
|
||||
indicial_ideal, _ = self._complete_indicial_ideal(max_multiplier_degree)
|
||||
return indicial_ideal
|
||||
|
||||
def indicials(self, max_multiplier_degree: int = 4) -> list:
|
||||
"""
|
||||
Returns the indicials (exponents) of this PFIdeal at z = 0, i.e. the roots of its indicial ideal, as
|
||||
a sorted list of (rho, multiplicity) pairs, the multiplicities adding up to the holonomic rank. E.g.
|
||||
the quintic has [((0,), 4)] at its MUM point and [((0,), 1), ((1,), 2), ((2,), 1)] at its conifold point.
|
||||
"""
|
||||
return _indicial_roots(self.indicial_ideal(max_multiplier_degree))
|
||||
|
||||
def frobenius_basis(self, order: int, max_log_degree: int = None, max_multiplier_degree: int = 4) -> list:
|
||||
"""
|
||||
Computes a basis of all solutions of this PFIdeal at z = 0, logarithmic ones included, as Periods
|
||||
to the given order. For the basis at another point, move that point to the origin first with
|
||||
change_coordinates.
|
||||
|
||||
The solutions are found one exponent class at a time: the indicials are grouped by their residue
|
||||
mod Z^n, and each class is solved once, from its componentwise smallest exponent rho_min (which
|
||||
becomes the Periods' indicials), so that solutions starting at rho_min + k with k > 0 - e.g. at the
|
||||
indicials 0, 1, 1, 2 of a conifold point - come out of the same solve. Their degree in log(z) is at
|
||||
most max_log_degree if given, and otherwise at most the class' total multiplicity minus 1 (the
|
||||
logarithms of the monodromies' unipotent parts act on the class' solutions as commuting nilpotent
|
||||
maps, so any product of that many of them vanishes) and at most n = len(hodge_numbers) - 1.
|
||||
|
||||
Within each class the basis is canonical (see _solve_exponent_class): every solution has its own
|
||||
leading term with coefficient 1, at which the class' other solutions vanish. Coefficients are those
|
||||
of plain log(z_i)^b, without factors of 2*pi*i or factorials.
|
||||
|
||||
This assumes the origin is a normal-crossing point of the discriminant, where every solution is a
|
||||
sum of terms z^rho * (power series) * (polynomial in log z). Elsewhere fewer solutions than the
|
||||
rank turn up, which is logged as a warning - as is a log structure at a MUM point that differs from
|
||||
hodge_numbers.
|
||||
"""
|
||||
n = self.no_variables
|
||||
indicial_ideal, operator_parts = self._complete_indicial_ideal(max_multiplier_degree)
|
||||
theta_ring = indicial_ideal.ring()
|
||||
|
||||
roots = _indicial_roots(indicial_ideal)
|
||||
classes = {}
|
||||
for rho, multiplicity in roots:
|
||||
classes.setdefault(tuple(r - floor(r) for r in rho), []).append((rho, multiplicity))
|
||||
exponent_classes = sorted(
|
||||
((tuple(min(rho[i] for rho, _ in members) for i in range(n)), members) for members in classes.values()),
|
||||
key=lambda exponent_class: exponent_class[0],
|
||||
)
|
||||
|
||||
cy_dimension = len(self.hodge_numbers) - 1 if self.hodge_numbers is not None else None
|
||||
basis = []
|
||||
for rho_min, members in exponent_classes:
|
||||
class_size = sum(multiplicity for _, multiplicity in members)
|
||||
# A solution starting at the indicial rho only turns up once the solve reaches z^(rho - rho_min).
|
||||
required_order = max(sum(rho[i] - rho_min[i] for i in range(n)) for rho, _ in members)
|
||||
if order < required_order:
|
||||
raise ValueError(
|
||||
"Order %d is too low to reach the indicials %s; it must be at least %d."
|
||||
% (order, [rho for rho, _ in members], required_order)
|
||||
)
|
||||
if max_log_degree is not None:
|
||||
log_degree = max_log_degree
|
||||
elif cy_dimension is not None:
|
||||
log_degree = min(class_size - 1, cy_dimension)
|
||||
else:
|
||||
log_degree = class_size - 1
|
||||
|
||||
solutions = self._solve_exponent_class(operator_parts, theta_ring, rho_min, order, log_degree)
|
||||
logger.info(
|
||||
"Exponent class of %s: %d solution(s) with leading terms %s.",
|
||||
rho_min, len(solutions), [_leading_term(coefficients, rho_min) for coefficients in solutions],
|
||||
)
|
||||
if len(solutions) != class_size:
|
||||
logger.warning(
|
||||
"Found %d solution(s) in the exponent class of %s, whose indicials have total multiplicity "
|
||||
"%d: z = 0 may not be a normal-crossing point of the discriminant, or the degree in log(z), "
|
||||
"capped at %d, is too low.",
|
||||
len(solutions), rho_min, class_size, log_degree,
|
||||
)
|
||||
basis.extend(
|
||||
Period(no_variables=n, coefficients=coefficients, order=order, indicials=list(rho_min))
|
||||
for coefficients in solutions
|
||||
)
|
||||
|
||||
if self.hodge_numbers is not None:
|
||||
if len(basis) != sum(self.hodge_numbers):
|
||||
logger.warning(
|
||||
"Found %d solution(s), but the rank sum(hodge_numbers) is %d.", len(basis), sum(self.hodge_numbers)
|
||||
)
|
||||
# At a point of maximal unipotent monodromy - all indicials 0, with a solution of top log-degree
|
||||
# n - the numbers of solutions of each top log-degree are the Hodge numbers. A single exponent
|
||||
# class at 0 isn't enough: at (0, -1/27) of P_{333111}[12], on the MUM divisor z0 = 0, the
|
||||
# indicials are (0, 0) and (0, 1) and there is a log(z0)^4, but no MUM point.
|
||||
top_log_degrees = [max(sum(log_index) for log_index in period.coefficients) for period in basis]
|
||||
if (
|
||||
len(roots) == 1
|
||||
and all(r == 0 for r in roots[0][0])
|
||||
and max(top_log_degrees, default=-1) == cy_dimension
|
||||
):
|
||||
log_structure = [top_log_degrees.count(degree) for degree in range(cy_dimension + 1)]
|
||||
if log_structure != self.hodge_numbers:
|
||||
logger.warning(
|
||||
"The numbers of solutions per top log-degree %s at this MUM point differ from "
|
||||
"hodge_numbers %s.", log_structure, self.hodge_numbers,
|
||||
)
|
||||
|
||||
logger.info("Found a Frobenius basis of %d solution(s) at order %d.", len(basis), order)
|
||||
return basis
|
||||
|
||||
def find_power_series_solution(self, indicials: list, order: int) -> list:
|
||||
"""
|
||||
Finds the power series solutions (no logs) to this PFIdeal at given indicial exponents/order, i.e.
|
||||
all solutions of the form z^indicials * (power series in z); frobenius_basis finds all solutions,
|
||||
logarithmic ones included.
|
||||
|
||||
This is _solve_exponent_class with no room for logarithms, on the operators as given (without
|
||||
completing their indicial ideal). A combination of solutions that could only be continued with a
|
||||
log(z)-term, at a resonance, is not a power series solution and is dropped.
|
||||
"""
|
||||
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]
|
||||
theta_ring, operator_parts = self._operator_theta_parts()
|
||||
|
||||
solutions = [
|
||||
Period(no_variables=self.no_variables, coefficients=coefficients, order=order, indicials=indicials)
|
||||
for coefficients in self._solve_exponent_class(operator_parts, theta_ring, indicials, order, 0)
|
||||
]
|
||||
|
||||
logger.info(
|
||||
"Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order
|
||||
)
|
||||
return solutions
|
||||
|
||||
def _operator_theta_parts(self) -> tuple:
|
||||
# The theta-ring QQ[theta0, ..., theta{n-1}] and, in it, the theta-parts (see
|
||||
# PFOperator._theta_parts) of every nonzero operator of this PFIdeal.
|
||||
theta_ring = PolynomialRing(QQ, self.no_variables, ['theta%d' % i for i in range(self.no_variables)])
|
||||
operator_parts = [parts for parts in (op._theta_parts(theta_ring) for op in self.operators) if parts]
|
||||
if not operator_parts:
|
||||
raise ValueError("The PFIdeal contains no nonzero operator.")
|
||||
return theta_ring, operator_parts
|
||||
|
||||
def _complete_indicial_ideal(self, max_multiplier_degree: int = 4) -> tuple:
|
||||
"""
|
||||
Computes the indicial ideal of this PFIdeal at z = 0, together with the operators needed to see it.
|
||||
|
||||
The indicial ideal consists of the theta-polynomials Q for which some operator in the (left) ideal
|
||||
has initial part z^e * Q(theta), i.e. this as its part at the lowest z-monomial w.r.t. _z_order_key.
|
||||
Its roots are the indicials (exponents) of the solutions at z = 0, and for a regular holonomic system,
|
||||
as every Picard-Fuchs system is, its colength - the number of roots counted with multiplicity - is
|
||||
the holonomic rank. See M. Saito, B. Sturmfels, N. Takayama, "Gröbner Deformations of
|
||||
Hypergeometric Differential Equations", Algorithms and Computation in Mathematics vol. 6,
|
||||
Springer, 2000, ch. 2.
|
||||
|
||||
The initial parts of the generators alone can generate a strictly smaller ideal. E.g. for
|
||||
P_{22211}[8] moved to (z0, z1) = (0, 1/4) they are theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1),
|
||||
which vanish on the whole line theta1 = 0. The missing elements are initial parts of combinations of
|
||||
the generators in which the generators' own initial parts cancel (S-pairs), and are found here by
|
||||
linear algebra: the left multiples z^a * theta^c * L of the generators L with |a| + |c| <= d span a
|
||||
finite-dimensional space (a Macaulay matrix), whose echelon form w.r.t. the order on z-monomials
|
||||
exposes every initial part it contains. d grows until the result is zero-dimensional and complete:
|
||||
of colength rank = sum(hodge_numbers) if those are known, and otherwise unchanged from d - 1.
|
||||
|
||||
Returns (indicial ideal, operator parts): the theta-parts (see PFOperator._theta_parts) of the
|
||||
generators, followed by those of the few extra operators whose initial parts were needed - solving
|
||||
with all of them sees exactly the right starting terms.
|
||||
"""
|
||||
n = self.no_variables
|
||||
theta_ring, generator_parts = self._operator_theta_parts()
|
||||
rank = sum(self.hodge_numbers) if self.hodge_numbers is not None else None
|
||||
generator_initials = [parts[min(parts, key=_z_order_key)] for parts in generator_parts]
|
||||
# A single ordinary differential operator needs no completion: initial parts are multiplicative, so
|
||||
# its indicial ideal is generated by its own initial part (Fuchs' indicial polynomial).
|
||||
single_ode = n == 1 and len(generator_parts) == 1
|
||||
|
||||
previous = None
|
||||
for degree in range(max_multiplier_degree + 1):
|
||||
# One row per left multiple, one column per (z-monomial, theta-monomial), the z-monomials in
|
||||
# z-order, so that each row of the echelon form starts at its initial z-monomial.
|
||||
rows = [
|
||||
_left_multiply(parts, multiplier[:n], multiplier[n:], theta_ring)
|
||||
for parts in generator_parts
|
||||
for multiplier in _multi_indices(2 * n, degree)
|
||||
]
|
||||
columns = sorted(
|
||||
{(p, tuple(q)) for row in rows for p, P in row.items() for q in P.dict()},
|
||||
key=lambda column: (_z_order_key(column[0]), column[1]),
|
||||
)
|
||||
column_position = {column: j for j, column in enumerate(columns)}
|
||||
# Dense on purpose, here and in _solve_exponent_class: Sage's echelon form of *sparse* rational
|
||||
# matrices can be wrong (seen with Sage 9.7 for a 1 x 2 matrix met when solving the quintic at
|
||||
# its orbifold point).
|
||||
macaulay = matrix(QQ, len(rows), len(columns))
|
||||
for i, row in enumerate(rows):
|
||||
for p, P in row.items():
|
||||
for q, coeff in P.dict().items():
|
||||
macaulay[i, column_position[(p, tuple(q))]] = coeff
|
||||
echelon = macaulay.echelon_form()
|
||||
|
||||
# Keep, greedily, every echelon row whose initial part isn't in the ideal yet, then drop again
|
||||
# those made redundant by rows kept later: every extra operator costs a full pass of the
|
||||
# solver, and very few are needed (one for P_{22211}[8] at (0, 1/4)).
|
||||
extra_parts, extra_initials = [], []
|
||||
current = theta_ring.ideal(generator_initials)
|
||||
for i, pivot in enumerate(echelon.pivots()):
|
||||
row = echelon.row(i).dict()
|
||||
initial_z_index = columns[pivot][0]
|
||||
initial = sum(
|
||||
coeff * theta_ring.monomial(*columns[j][1])
|
||||
for j, coeff in row.items()
|
||||
if columns[j][0] == initial_z_index
|
||||
)
|
||||
if current.reduce(initial) == 0:
|
||||
continue
|
||||
parts = {}
|
||||
for j, coeff in row.items():
|
||||
p, q = columns[j]
|
||||
parts[p] = parts.get(p, theta_ring(0)) + coeff * theta_ring.monomial(*q)
|
||||
extra_parts.append(parts)
|
||||
extra_initials.append(initial)
|
||||
current = theta_ring.ideal(generator_initials + extra_initials)
|
||||
for i in reversed(range(len(extra_parts))):
|
||||
others = generator_initials + extra_initials[:i] + extra_initials[i + 1:]
|
||||
if theta_ring.ideal(others).reduce(extra_initials[i]) == 0:
|
||||
del extra_parts[i], extra_initials[i]
|
||||
|
||||
indicial_ideal = theta_ring.ideal(generator_initials + extra_initials)
|
||||
finite = indicial_ideal.dimension() <= 0
|
||||
colength = indicial_ideal.vector_space_dimension() if finite else None
|
||||
logger.debug(
|
||||
"Indicial ideal at multiplier degree %d with %d extra operator(s): %s.", degree, len(extra_parts),
|
||||
"colength %d" % colength if finite else "not zero-dimensional",
|
||||
)
|
||||
if finite and rank is not None:
|
||||
if colength == rank:
|
||||
break
|
||||
if colength < rank:
|
||||
raise ValueError(
|
||||
"The indicial ideal at z = 0 has colength %d, less than the rank sum(hodge_numbers) = %d: "
|
||||
"an operator is not in the Picard-Fuchs ideal, or hodge_numbers is wrong." % (colength, rank)
|
||||
)
|
||||
elif finite and (single_ode or indicial_ideal == previous):
|
||||
if not single_ode:
|
||||
logger.warning(
|
||||
"Taking the indicial ideal at z = 0 as complete since it didn't change from multiplier "
|
||||
"degree %d to %d; give the PFIdeal hodge_numbers to certify this by its rank.",
|
||||
degree - 1, degree,
|
||||
)
|
||||
break
|
||||
previous = indicial_ideal
|
||||
else:
|
||||
raise ValueError(
|
||||
"Could not complete the indicial ideal at z = 0 up to multiplier degree %d: the operators may not "
|
||||
"generate the whole Picard-Fuchs ideal (whose rank then exceeds sum(hodge_numbers)), or a higher "
|
||||
"max_multiplier_degree is needed." % max_multiplier_degree
|
||||
)
|
||||
|
||||
logger.info(
|
||||
"Completed the indicial ideal at z = 0 at multiplier degree %d with %d extra operator(s): colength %d.",
|
||||
degree, len(extra_parts), colength,
|
||||
)
|
||||
return indicial_ideal, generator_parts + extra_parts
|
||||
|
||||
def _solve_exponent_class(
|
||||
self, operator_parts: list, theta_ring, indicials, order: int, max_log_degree: int
|
||||
) -> list:
|
||||
"""
|
||||
Finds all solutions y = sum_k z^(indicials + k) * g_k(log z) of the operators given by their
|
||||
theta-parts (see PFOperator._theta_parts), to the given order: k runs over the multi-indices >= 0,
|
||||
and every g_k is a polynomial in L_i = log(z_i) of total degree at most max_log_degree. Returns one
|
||||
coefficients dictionary {log_index: {z_index: coeff}} (see Period) per solution of a canonical basis.
|
||||
|
||||
Since theta_i acts on z^(indicials + k) * g(L) as indicials_i + k_i + d/dL_i, a term c*z^p*theta^q
|
||||
sends it to c*z^(indicials + k + p) * (indicials + k + d/dL)^q g(L): it shifts k up by p (p >= 0),
|
||||
never down or sideways. So, as in Frobenius' method, the solutions can be found z-monomial by
|
||||
z-monomial - only with every unknown a polynomial g_k rather than a number, on whose coefficients a
|
||||
theta-polynomial P acts by the matrix of P(indicials + k + d/dL). Logarithms, indicials differing by
|
||||
integers (resonance) and several operators constraining each other are then all handled by the
|
||||
kernels and consistency conditions of these small matrices, with no derivatives with respect to the
|
||||
exponent and no case analysis.
|
||||
|
||||
The basis is canonical: in the reduced echelon form of all solutions, each has its own leading
|
||||
coefficient 1 at which all others vanish, so in particular no solution's tail contains another one's
|
||||
leading term - as for the canonical series solutions of M. Saito, B. Sturmfels, N. Takayama,
|
||||
"Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and Computation in
|
||||
Mathematics vol. 6, Springer, 2000, ch. 2.
|
||||
"""
|
||||
n = self.no_variables
|
||||
theta_gens = theta_ring.gens()
|
||||
indicials = [QQ(rho) for rho in indicials]
|
||||
|
||||
# Monomials L^b spanning the polynomials of total degree <= max_log_degree in L, lowest degree first
|
||||
# (1, L0, L1, L0^2, L0*L1, ...); this order decides which coefficient of a solution is normalised.
|
||||
log_indices = sorted(_multi_indices(n, max_log_degree), key=lambda b: (sum(b), tuple(-x for x in b)))
|
||||
log_position = {log_index: i for i, log_index in enumerate(log_indices)}
|
||||
size = len(log_indices)
|
||||
|
||||
# By Taylor's formula, P(s + d/dL) sends L^b to the sum over b' <= b of
|
||||
# binom(b, b') * (d^(b - b')P/dtheta^(b - b'))(s) * L^b'. With this pattern of entries and every
|
||||
# theta-polynomial's derivatives computed once, each matrix below costs only a few evaluations.
|
||||
action_pattern = [
|
||||
(
|
||||
log_position[lower],
|
||||
log_position[upper],
|
||||
tuple(upper[i] - lower[i] for i in range(n)),
|
||||
prod(binomial(upper[i], lower[i]) for i in range(n)),
|
||||
)
|
||||
for upper in log_indices
|
||||
for lower in log_indices
|
||||
if all(lower[i] <= upper[i] for i in range(n))
|
||||
]
|
||||
|
||||
def taylor_derivatives(P):
|
||||
derivatives = {}
|
||||
for alpha in _multi_indices(n, max_log_degree):
|
||||
derivative = P
|
||||
for i in range(n):
|
||||
if alpha[i]:
|
||||
derivative = derivative.derivative(theta_gens[i], alpha[i])
|
||||
if derivative != 0:
|
||||
derivatives[alpha] = derivative
|
||||
return derivatives
|
||||
|
||||
def action_matrix(derivatives, shift):
|
||||
values = {alpha: derivative(*shift) for alpha, derivative in derivatives.items()}
|
||||
entries = {}
|
||||
for row, column, alpha, binomial_factor in action_pattern:
|
||||
if values.get(alpha, 0) != 0:
|
||||
entries[(row, column)] = binomial_factor * values[alpha]
|
||||
# matrix() would make a sparse matrix from a dictionary; see the comment on dense matrices in
|
||||
# _complete_indicial_ideal.
|
||||
return matrix(QQ, size, size, entries, sparse=False)
|
||||
|
||||
operator_derivatives = [{p: taylor_derivatives(P) for p, P in parts.items()} for parts in operator_parts]
|
||||
supports = [sorted(parts, key=_z_order_key) for parts in operator_parts]
|
||||
z_indices = sorted(_multi_indices(n, order), key=_z_order_key)
|
||||
|
||||
# solved[k] holds g_k as a size x dimension matrix, whose column j holds g_k's coefficients
|
||||
# (w.r.t. log_indices) in the j-th of the `dimension` independent solutions found so far.
|
||||
dimension = 0
|
||||
solved = {}
|
||||
|
||||
for k in z_indices:
|
||||
# Collect the equations whose latest unknown is g_k. The coefficient of z^(indicials + m) in L(y)
|
||||
# involves g_(m - p) for all p <= m (componentwise) in L's support, the latest of them in z-order
|
||||
# being g_(m - p') for the smallest such p'. Taking m = k + p' for those p' that are the smallest
|
||||
# element p <= k + p' of the support thus visits every equation exactly once - also for operators
|
||||
# whose lowest z-monomial isn't 1, like the extra ones from _complete_indicial_ideal.
|
||||
diagonal_blocks, known_blocks = [], []
|
||||
for derivatives, support in zip(operator_derivatives, supports):
|
||||
for position, pivot_exponent in enumerate(support):
|
||||
m = tuple(k[i] + pivot_exponent[i] for i in range(n))
|
||||
if any(all(p[i] <= m[i] for i in range(n)) for p in support[:position]):
|
||||
continue
|
||||
diagonal_blocks.append(
|
||||
action_matrix(derivatives[pivot_exponent], [indicials[i] + k[i] for i in range(n)])
|
||||
)
|
||||
known = matrix(QQ, size, dimension)
|
||||
if dimension:
|
||||
for p in support[position + 1:]:
|
||||
if all(p[i] <= m[i] for i in range(n)):
|
||||
j = tuple(m[i] - p[i] for i in range(n)) # earlier than k in z-order: solved
|
||||
known += action_matrix(derivatives[p], [indicials[i] + j[i] for i in range(n)]) * solved[j]
|
||||
known_blocks.append(known)
|
||||
diagonal = block_matrix(QQ, [[block] for block in diagonal_blocks], subdivide=False)
|
||||
known = block_matrix(QQ, [[block] for block in known_blocks], subdivide=False)
|
||||
|
||||
# Row-reduce [diagonal | known]: rows with a pivot in the diagonal part give g_k in terms of the
|
||||
# solutions so far, rows without are linear conditions on those solutions, and the diagonal
|
||||
# part's free columns start new solutions at z^(indicials + k), where the indicial ideal vanishes.
|
||||
reduced = diagonal.augment(known).rref()
|
||||
conditions = [row[size:] for row in reduced.rows() if row[:size] == 0 and row[size:] != 0]
|
||||
if conditions:
|
||||
# Some combinations of the solutions so far can't be continued past this point - e.g. as they
|
||||
# would need a higher power of log(z) than max_log_degree. Keep those that can, by rewriting
|
||||
# every g_j in a basis of the combinations satisfying all conditions.
|
||||
restriction = matrix(QQ, conditions).right_kernel_matrix().transpose()
|
||||
for j in solved:
|
||||
solved[j] = solved[j] * restriction
|
||||
dimension = restriction.ncols()
|
||||
reduced = diagonal.augment(known * restriction).rref()
|
||||
|
||||
pivots = [column for column in reduced.pivots() if column < size]
|
||||
free = [column for column in range(size) if column not in pivots]
|
||||
# The free coefficients are 0 in the continuations of the solutions so far, so that no tail
|
||||
# contains the leading term of a solution starting here; each new solution has one of them 1.
|
||||
particular = matrix(QQ, size, dimension)
|
||||
new_solutions = matrix(QQ, size, len(free))
|
||||
for row_index, column in enumerate(pivots):
|
||||
particular.set_row(column, -reduced.row(row_index)[size:])
|
||||
for j, free_column in enumerate(free):
|
||||
new_solutions[column, j] = -reduced[row_index, free_column]
|
||||
for j, free_column in enumerate(free):
|
||||
new_solutions[free_column, j] = 1
|
||||
|
||||
if free:
|
||||
for j in solved:
|
||||
solved[j] = solved[j].augment(matrix(QQ, size, len(free)))
|
||||
solved[k] = particular.augment(new_solutions)
|
||||
dimension += len(free)
|
||||
|
||||
if dimension == 0:
|
||||
return []
|
||||
|
||||
# Canonical basis: the reduced echelon form of all solutions' coefficients, ordered by z-monomial (in
|
||||
# z-order) and then log-monomial. Each solution then has its own leading coefficient 1 at which all
|
||||
# others vanish, independently of the order in which the solutions turned up above.
|
||||
coefficient_labels = [(k, log_index) for k in z_indices for log_index in log_indices]
|
||||
echelon = block_matrix(QQ, [[solved[k]] for k in z_indices], subdivide=False).transpose().rref()
|
||||
|
||||
solutions = []
|
||||
for solution_vector in echelon.rows()[:dimension]:
|
||||
by_log_index = {}
|
||||
for (k, log_index), coeff in zip(coefficient_labels, solution_vector):
|
||||
by_log_index.setdefault(log_index, {})[k] = coeff
|
||||
# Every z-index up to order is kept, zeros included, for each log-monomial that occurs.
|
||||
solutions.append({
|
||||
log_index: z_dict for log_index, z_dict in by_log_index.items()
|
||||
if any(coeff != 0 for coeff in z_dict.values())
|
||||
})
|
||||
return solutions
|
||||
|
||||
|
||||
def _z_order_key(z_index) -> tuple:
|
||||
# The order on z-monomials used throughout: by total degree, ties broken lexicographically. Being graded,
|
||||
# it is compatible with truncating at a total degree; being a monomial order (p < q implies p + r < q + r),
|
||||
# it acts like a generic weight vector, so an operator's initial (lowest) part is a single z-monomial
|
||||
# times a theta-polynomial, and solving in this order only ever needs earlier unknowns.
|
||||
return (sum(z_index),) + tuple(z_index)
|
||||
|
||||
|
||||
def _left_multiply(parts: dict, z_shift, theta_power, theta_ring) -> dict:
|
||||
# Normal-ordered product (z^a * theta^c) * (sum_p z^p * P_p(theta)) = sum_p z^(a + p) * (theta + p)^c * P_p(theta)
|
||||
# of theta-parts (see PFOperator._theta_parts), by theta_i * z^p = z^p * (theta_i + p_i).
|
||||
theta_gens = theta_ring.gens()
|
||||
n = len(theta_gens)
|
||||
product = {}
|
||||
for p, P in parts.items():
|
||||
shifted_power = prod((theta_gens[i] + p[i]) ** theta_power[i] for i in range(n))
|
||||
product[tuple(z_shift[i] + p[i] for i in range(n))] = shifted_power * P
|
||||
return product
|
||||
|
||||
|
||||
def _indicial_roots(indicial_ideal) -> list:
|
||||
# The roots of a zero-dimensional indicial ideal with their multiplicities, as a sorted list of
|
||||
# (rho, multiplicity): each primary component belongs to one root, and its colength is that multiplicity.
|
||||
theta_gens = indicial_ideal.ring().gens()
|
||||
roots = []
|
||||
for primary, prime in indicial_ideal.complete_primary_decomposition():
|
||||
points = prime.variety()
|
||||
if len(points) != 1 or prime.vector_space_dimension() != 1:
|
||||
raise NotImplementedError(
|
||||
"The indicial ideal has a non-rational root (prime component %s), but Picard-Fuchs systems of "
|
||||
"geometric families only have rational indicials (quasi-unipotent monodromy)." % (prime,)
|
||||
)
|
||||
roots.append((tuple(points[0][theta] for theta in theta_gens), primary.vector_space_dimension()))
|
||||
return sorted(roots)
|
||||
|
||||
|
||||
def _leading_term(coefficients: dict, indicials) -> str:
|
||||
# The terms of a solution (a coefficients dictionary, see Period) at its z-order-smallest z-monomial,
|
||||
# e.g. 'sqrt(z1)*(log(z0) + log(z1))', to log the structure of a basis.
|
||||
n = len(indicials)
|
||||
z_gens = [SR.var('z%d' % i) for i in range(n)]
|
||||
leading = min(
|
||||
(z_index for z_dict in coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0),
|
||||
key=_z_order_key,
|
||||
)
|
||||
log_polynomial = sum(
|
||||
z_dict[leading] * prod(log(z_gens[i]) ** log_index[i] for i in range(n))
|
||||
for log_index, z_dict in coefficients.items()
|
||||
)
|
||||
return str(prod(z_gens[i] ** (indicials[i] + leading[i]) for i in range(n)) * log_polynomial)
|
||||
@@ -1,349 +0,0 @@
|
||||
"""
|
||||
Periods of Calabi-Yau manifolds as formal power series in z-variables and their logs.
|
||||
|
||||
Contents
|
||||
--------
|
||||
Period
|
||||
__init__ Construct from a string or a coefficients dictionary
|
||||
from_l_vectors (classmethod) Holomorphic period from GLSM/Mori charge vectors
|
||||
apply_operator Apply a PFOperator, truncated to self.order
|
||||
find_annihilating_operators Find PFOperator Ansätze annihilating this period
|
||||
_initialise_ring Polynomial ring in z- and log-variables
|
||||
_period_from_string Parse a string into a coefficients dictionary
|
||||
_period_to_string Render the coefficients as a string
|
||||
_max_z_degree Highest total z-degree among the coefficients
|
||||
_truncate_coefficients Drop terms above a given total z-degree
|
||||
_apply_theta Apply theta_i = z_i * d/dz_i to a coefficients dictionary
|
||||
|
||||
PeriodAnsatz(Period)
|
||||
__init__ Ansatz of given z- and log-multi-degrees
|
||||
"""
|
||||
import logging
|
||||
from sage.all import sage_eval, PolynomialRing, QQ, SR, ceil, factorial, floor, log, matrix, prod, var
|
||||
|
||||
load("sage/operators.sage")
|
||||
|
||||
# Logger
|
||||
logger = logging.getLogger(__name__)
|
||||
logging.basicConfig(level=logging.DEBUG)
|
||||
|
||||
|
||||
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)
|
||||
|
||||
@classmethod
|
||||
def from_l_vectors(cls, l_vectors: list, order: int) -> "Period":
|
||||
"""
|
||||
Computes the holomorphic (fundamental) period directly from the GLSM/Mori charge
|
||||
vectors l^(1), ..., l^(m) of a toric Calabi-Yau (m = number of moduli), via the
|
||||
hypergeometric-type sum
|
||||
|
||||
varpi_0(z) = sum_{a in Z_{>=0}^m}
|
||||
(prod_j Gamma(1 - l_num^(j).a)) / (prod_i Gamma(1 + l_den^(i).a)) * z^a.
|
||||
"""
|
||||
if any(len(l_vector) != 2 for l_vector in l_vectors):
|
||||
raise ValueError("Each charge vector must be given as [numerator part, denominator part].")
|
||||
m = len(l_vectors)
|
||||
numerator_charges = [l_vector[0] for l_vector in l_vectors]
|
||||
denominator_charges = [l_vector[1] for l_vector in l_vectors]
|
||||
len0 = len(numerator_charges[0]) if m else 0
|
||||
len1 = len(denominator_charges[0]) if m else 0
|
||||
|
||||
# Each denominator Gamma-factor is linear in a; pivot it on the smallest index j with
|
||||
# a nonzero coefficient there. Since a is built up from a_{m-1} down to a_0, every
|
||||
# other index still appearing in that linear form is then already fixed, turning the
|
||||
# factor's pole at non-positive integers into a concrete bound on a_j.
|
||||
pivots = [[] for _ in range(m)]
|
||||
for i in range(len1):
|
||||
charges = [denominator_charges[k][i] for k in range(m)]
|
||||
j = next(k for k in range(m) if charges[k] != 0)
|
||||
pivots[j].append(charges)
|
||||
|
||||
coefficients = {}
|
||||
a = [0] * m
|
||||
|
||||
def recurse(j, num_partial, den_partial, degree):
|
||||
if j < 0:
|
||||
numerator = prod(factorial(-x) for x in num_partial)
|
||||
denominator = prod(factorial(x) for x in den_partial)
|
||||
coefficients[tuple(a)] = QQ(numerator) / QQ(denominator)
|
||||
return
|
||||
|
||||
lower, upper = 0, order - degree
|
||||
for charges in pivots[j]:
|
||||
c = charges[j]
|
||||
bound = -QQ(sum(charges[k] * a[k] for k in range(j + 1, m))) / QQ(c)
|
||||
if c < 0:
|
||||
upper = min(upper, floor(bound))
|
||||
else:
|
||||
lower = max(lower, ceil(bound))
|
||||
|
||||
for a_j in range(lower, upper + 1):
|
||||
a[j] = a_j
|
||||
recurse(
|
||||
j - 1,
|
||||
[num_partial[jj] + numerator_charges[j][jj] * a_j for jj in range(len0)],
|
||||
[den_partial[ii] + denominator_charges[j][ii] * a_j for ii in range(len1)],
|
||||
degree + a_j,
|
||||
)
|
||||
|
||||
recurse(m - 1, [0] * len0, [0] * len1, 0)
|
||||
|
||||
logger.info("Computed holomorphic period from l-vectors for %d moduli at order %d.", m, order)
|
||||
return cls(no_variables=m, coefficients={tuple([0] * m): coefficients}, order=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()),
|
||||
)
|
||||
@@ -1,48 +1,3 @@
|
||||
"""
|
||||
Topological data of Calabi-Yau hypersurfaces and CICYs in toric ambient spaces.
|
||||
|
||||
Contents
|
||||
--------
|
||||
Polytope
|
||||
__init__ Construct from a list of defining points
|
||||
_validate_points_and_get_dimension Check the points and return their dimension
|
||||
_points_to_polytope Lattice points of the convex hull, minus codim-1 interiors
|
||||
|
||||
ToricPolytope(Polytope)
|
||||
__init__ Construct with triangulation, nef partition and Mori cone
|
||||
set_triangulation Select a fine star triangulation by index
|
||||
set_nef_partition Set and validate the nef partition
|
||||
set_Mori_cone Set the Mori cone (default or from given l-vectors)
|
||||
disc A-discriminant (singular loci in Batyrev coordinates)
|
||||
topdata Compute topological data and write it to JSON
|
||||
_set_fine_star_triangulations Fine star triangulations of the polytope
|
||||
_collect_linear_dependencies Linear relations among points in all faces
|
||||
_append_origin_weight Add the origin weight to each relation
|
||||
_setup_discriminant_symbols Symbolic variables for the discriminant
|
||||
_build_equation_system Equation system for one block of relations
|
||||
_eliminate_and_normalize Eliminate lambda-variables and normalise
|
||||
_append_unique_discriminants Collect discriminant factors without duplicates
|
||||
_compute_Kahler_generators Kähler cone generators dual to the Mori cone
|
||||
_compute_intersection_ring Intersection ring polynomial of the CY
|
||||
_intersection_substitution_rules Kähler monomial -> intersection number rules
|
||||
_replace_intersection_monomials Substitute intersection numbers into an expression
|
||||
_compute_Chern_character Total Chern class via adjunction
|
||||
_compute_Chern_polynomials Chern polynomials with intersection numbers
|
||||
_check_for_fibrations Detect possible elliptic fibrations
|
||||
_get_polytope_and_glsm_table Table of polytope points and GLSM charges
|
||||
_get_lvec_rows GLSM charge l-vectors as lists of ints
|
||||
_get_output Dictionary of all topological data
|
||||
_to_json_safe (staticmethod) Convert values into JSON-safe primitives
|
||||
_write_output_json Write the output dictionary to JSON
|
||||
|
||||
ToricPolytopeCICY(ToricPolytope)
|
||||
__init__ Construct from a CICY configuration matrix
|
||||
_CICY_to_points Configuration matrix -> points and nef partition
|
||||
|
||||
ToricPolytopeProjectiveSpace(ToricPolytope)
|
||||
__init__ Construct from weighted projective space weights
|
||||
_compute_points_from_weights Weights -> polytope points
|
||||
"""
|
||||
import numpy as np
|
||||
import logging
|
||||
import os
|
||||
@@ -146,7 +101,6 @@ class ToricPolytope(Polytope):
|
||||
if no_triangulation < 0 or no_triangulation >= len(self.triangulations):
|
||||
raise IndexError("Invalid triangulation index {} (available: 0..{}).".format(no_triangulation, len(self.triangulations) - 1))
|
||||
|
||||
if len(self.triangulations) > 1:
|
||||
logger.info("Using triangulation index %d of [0..%d].", no_triangulation, len(self.triangulations) - 1)
|
||||
self.triangulation = self.triangulations[no_triangulation]
|
||||
|
||||
|
||||
@@ -1,15 +0,0 @@
|
||||
"""
|
||||
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))
|
||||
@@ -2,10 +2,10 @@
|
||||
# expected_intersection_numbers_CY), where expected_disc is a list of
|
||||
# (str(discriminant_factor), codimension) pairs as returned by ToricPolytope.disc().
|
||||
|
||||
TEST_MODELS = [
|
||||
PLAYGROUND_MODEL_FACTORIES = [
|
||||
(
|
||||
"K3_two_parameter_family",
|
||||
lambda: ToricPolytopeProjectiveSpace([4, 2, 1, 1], model_name="K3_two_parameter_family"),
|
||||
lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 4], model_name="K3_two_parameter_family"),
|
||||
2, # cy_dimension
|
||||
2, # no_divs
|
||||
[("z2 - 1/4", 2), ("4096*z1^2*(4*z2 - 1) + 128*z1 - 1", 0)],
|
||||
@@ -13,7 +13,7 @@ TEST_MODELS = [
|
||||
),
|
||||
(
|
||||
"CY3_two_parameter_family",
|
||||
lambda: ToricPolytopeProjectiveSpace([2, 2, 2, 1, 1], model_name="CY3_two_parameter_family"),
|
||||
lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 2, 2], model_name="CY3_two_parameter_family"),
|
||||
3,
|
||||
2,
|
||||
[("z2 - 1/4", 3), ("65536*z1^2*(4*z2 - 1) + 512*z1 - 1", 0)],
|
||||
@@ -21,7 +21,7 @@ TEST_MODELS = [
|
||||
),
|
||||
(
|
||||
"CY4_two_parameter_family",
|
||||
lambda: ToricPolytopeProjectiveSpace([12, 8, 1, 1, 1, 1], model_name="CY4_two_parameter_family"),
|
||||
lambda: ToricPolytopeProjectiveSpace([1, 1, 1, 1, 8, 12], model_name="CY4_two_parameter_family"),
|
||||
4,
|
||||
2,
|
||||
[
|
||||
@@ -1,447 +0,0 @@
|
||||
import logging
|
||||
import os
|
||||
|
||||
import pytest
|
||||
from sage.all import * # noqa: F401
|
||||
|
||||
load(os.path.join(os.path.dirname(__file__), "..", "sage", "periods.sage"))
|
||||
|
||||
|
||||
def _quintic_ideal():
|
||||
# The quintic's Picard-Fuchs operator, of rank 4 with log structure [1, 1, 1, 1] at its MUM point.
|
||||
return PFIdeal(
|
||||
[PFOperator("theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)", no_variables=1)],
|
||||
hodge_numbers=[1, 1, 1, 1],
|
||||
)
|
||||
|
||||
|
||||
def _p22211_ideal():
|
||||
# The Picard-Fuchs ideal of P_{22211}[8] at its MUM point (see
|
||||
# test_ideal_finds_power_series_solution_and_recovers_first_operator), with h^{2,1} = 2.
|
||||
D1 = PFOperator(
|
||||
"theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
|
||||
no_variables=2,
|
||||
)
|
||||
D2 = PFOperator(
|
||||
"theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
|
||||
no_variables=2,
|
||||
)
|
||||
return PFIdeal([D1, D2], hodge_numbers=[1, 2, 2, 1])
|
||||
|
||||
|
||||
def _p333111_ideal():
|
||||
# The Picard-Fuchs ideal of the fourfold P_{333111}[12] at its MUM point: the GKZ operators of
|
||||
# l^(1) = (-4; 1,1,1,0,0,0,1) and l^(2) = (0; 0,0,0,1,1,1,-3), the first without its left factor theta0.
|
||||
# With the intersection numbers kappa_0000 = 36, kappa_0001 = 12, kappa_0011 = 4 (and 0 else) that
|
||||
# toric_topdata computes, the vertical h^{2,2} is 3.
|
||||
D1 = PFOperator(
|
||||
"theta0^2*(theta0 - 3*theta1) - 4*z0*(4*theta0 + 1)*(4*theta0 + 2)*(4*theta0 + 3)",
|
||||
no_variables=2,
|
||||
)
|
||||
D2 = PFOperator(
|
||||
"theta1^3 - z1*(theta0 - 3*theta1)*(theta0 - 3*theta1 - 1)*(theta0 - 3*theta1 - 2)",
|
||||
no_variables=2,
|
||||
)
|
||||
return PFIdeal([D1, D2], hodge_numbers=[1, 2, 3, 2, 1])
|
||||
|
||||
|
||||
def _starting_term(period):
|
||||
# A solution's lowest z-monomial (by total degree, then lexicographically) together with its coefficient
|
||||
# per log-monomial there, e.g. ((1,), {(1,): 1}) for z0*log(z0) + ...
|
||||
z_index = min(
|
||||
(z_index for z_dict in period.coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0),
|
||||
key=lambda z_index: (sum(z_index), z_index),
|
||||
)
|
||||
return z_index, {
|
||||
log_index: z_dict[z_index] for log_index, z_dict in period.coefficients.items() if z_dict[z_index] != 0
|
||||
}
|
||||
|
||||
|
||||
def _nonzero(coefficients):
|
||||
# A coefficients dictionary without its zero entries, to compare bases that store zeros differently.
|
||||
return {
|
||||
log_index: {z_index: coeff for z_index, coeff in z_dict.items() if coeff != 0}
|
||||
for log_index, z_dict in coefficients.items()
|
||||
if any(coeff != 0 for coeff in z_dict.values())
|
||||
}
|
||||
|
||||
|
||||
def _truncated(f, degree):
|
||||
# A polynomial without its terms of total degree above degree.
|
||||
return sum((coeff * monomial for coeff, monomial in f if monomial.degree() <= degree), f.parent()(0))
|
||||
|
||||
|
||||
def _inverse_of_linear(constant, linear, degree):
|
||||
# 1/(constant + linear) for a linear form in eps, as a polynomial truncated at the given degree.
|
||||
return sum(((-linear / constant) ** j for j in range(degree + 1)), linear.parent()(0)) / constant
|
||||
|
||||
|
||||
def _classical_frobenius_periods(coefficient_ratio, starting_terms, order):
|
||||
# The classical Frobenius method. With coefficient_ratio(k) = c(k + eps)/c(eps) for the Gamma-function
|
||||
# coefficients c of the fundamental period (a polynomial in eps, truncated at the starting terms' degree),
|
||||
# the period starting with Q(log z) is Q(d/deps) applied to sum_k c(k + eps)/c(eps) * z^(k + eps) at eps = 0.
|
||||
n = starting_terms[0].parent().ngens()
|
||||
log_ring = PolynomialRing(QQ, n, ["L%d" % i for i in range(n)])
|
||||
logs = log_ring.gens()
|
||||
periods = [{} for _ in starting_terms]
|
||||
for k in _multi_indices(n, order):
|
||||
ratio = coefficient_ratio(k)
|
||||
for period, starting_term in zip(periods, starting_terms):
|
||||
# d^alpha/deps^alpha of ratio * exp(eps . L) at eps = 0, summed over the monomials eps^alpha of Q
|
||||
value = log_ring(0)
|
||||
for q, q_monomial in starting_term:
|
||||
alpha = q_monomial.exponents()[0]
|
||||
for r, r_monomial in ratio:
|
||||
beta = r_monomial.exponents()[0]
|
||||
if all(beta[i] <= alpha[i] for i in range(n)):
|
||||
value += q * r * prod(
|
||||
factorial(alpha[i]) // factorial(alpha[i] - beta[i]) * logs[i] ** (alpha[i] - beta[i])
|
||||
for i in range(n)
|
||||
)
|
||||
for coeff, log_monomial in value:
|
||||
period.setdefault(tuple(log_monomial.exponents()[0]), {})[k] = coeff
|
||||
return periods
|
||||
|
||||
|
||||
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] 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,
|
||||
)
|
||||
D2 = PFOperator(
|
||||
"theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
|
||||
no_variables=2,
|
||||
)
|
||||
moved = PFIdeal([D1, D2]).change_coordinates("z0", "z1 - 1/4")
|
||||
op1, op2 = moved.operators
|
||||
|
||||
solutions = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=14)
|
||||
assert len(solutions) == 1
|
||||
|
||||
period = solutions[0]
|
||||
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, 0)] == 0
|
||||
|
||||
recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
|
||||
assert len(recovered) == 1
|
||||
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():
|
||||
# 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
|
||||
|
||||
|
||||
|
||||
def test_indicials_of_quintic_at_mum_conifold_and_orbifold_points():
|
||||
# At the MUM point all four indicials are 0; at the conifold point they are 0, 1, 1, 2, the doubled 1
|
||||
# carrying the single logarithmic solution; at the orbifold point z0 -> 1/z0 they are 1/5, ..., 4/5.
|
||||
ideal = _quintic_ideal()
|
||||
|
||||
assert ideal.indicials() == [((0,), 4)]
|
||||
assert ideal.change_coordinates("1 - 3125*z0").indicials() == [((0,), 1), ((1,), 2), ((2,), 1)]
|
||||
assert ideal.change_coordinates("1/z0").indicials() == [((QQ(k) / 5,), 1) for k in range(1, 5)]
|
||||
|
||||
|
||||
def test_indicials_reject_hodge_numbers_inconsistent_with_the_operators():
|
||||
# The quintic's operator has rank 4, so its indicials can't add up to the rank 5 of these Hodge numbers.
|
||||
ideal = PFIdeal(_quintic_ideal().operators, hodge_numbers=[1, 1, 1, 1, 1])
|
||||
|
||||
with pytest.raises(ValueError):
|
||||
ideal.indicials()
|
||||
|
||||
|
||||
def test_frobenius_basis_of_quintic_at_mum_point():
|
||||
# The solutions start with 1, log(z0), log(z0)^2 and log(z0)^3, the first being the classic holomorphic
|
||||
# period. Each solution's top-log coefficient is a solution without logs, hence that period again.
|
||||
ideal = _quintic_ideal()
|
||||
|
||||
basis = ideal.frobenius_basis(order=4)
|
||||
|
||||
assert [_starting_term(period) for period in basis] == [((0,), {(j,): 1}) for j in range(4)]
|
||||
holomorphic = {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000, (4,): 305540235000}
|
||||
assert basis[0].coefficients == {(0,): holomorphic}
|
||||
assert all(period.coefficients[(j,)] == holomorphic for j, period in enumerate(basis))
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
|
||||
|
||||
|
||||
def test_frobenius_basis_of_quintic_at_conifold_point():
|
||||
# The indicials 0, 1, 1, 2 differ by integers, so all four solutions come out of one solve from exponent 0.
|
||||
# Picard-Lefschetz: exactly one of them is logarithmic, z0*log(z0) + ..., starting where the vanishing
|
||||
# period z0 + ... does.
|
||||
conifold = _quintic_ideal().change_coordinates("1 - 3125*z0")
|
||||
|
||||
basis = conifold.frobenius_basis(order=4)
|
||||
|
||||
assert [_starting_term(period) for period in basis] == [
|
||||
((0,), {(0,): 1}),
|
||||
((1,), {(0,): 1}),
|
||||
((1,), {(1,): 1}),
|
||||
((2,), {(0,): 1}),
|
||||
]
|
||||
assert [max(log_index[0] for log_index in period.coefficients) for period in basis] == [0, 0, 1, 0]
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in conifold.operators)
|
||||
|
||||
|
||||
def test_frobenius_basis_of_quintic_at_orbifold_point():
|
||||
# Four exponent classes 1/5, ..., 4/5 mod 1, each with a single solution without logs. Order 6 reaches
|
||||
# the step at which Sage's echelon form of sparse rational matrices went wrong (see operators.sage).
|
||||
orbifold = _quintic_ideal().change_coordinates("1/z0")
|
||||
|
||||
basis = orbifold.frobenius_basis(order=6)
|
||||
|
||||
assert [period.indicials for period in basis] == [[QQ(k) / 5] for k in range(1, 5)]
|
||||
assert all(list(period.coefficients) == [(0,)] for period in basis)
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in orbifold.operators)
|
||||
|
||||
|
||||
def test_frobenius_basis_of_two_parameter_model_at_mum_point():
|
||||
# All six indicials are (0, 0), and the solutions start with 1; log(z0), log(z1); log(z0)^2,
|
||||
# log(z0)*log(z1); and a triple-log term proportional to kappa_ijk*log(z_i)*log(z_j)*log(z_k) with the
|
||||
# intersection numbers kappa_000 = 8, kappa_001 = 4: the log structure [1, 2, 2, 1] of the Hodge numbers.
|
||||
ideal = _p22211_ideal()
|
||||
|
||||
assert ideal.indicials() == [((0, 0), 6)]
|
||||
basis = ideal.frobenius_basis(order=3)
|
||||
|
||||
assert [_starting_term(period) for period in basis] == [
|
||||
((0, 0), {(0, 0): 1}),
|
||||
((0, 0), {(1, 0): 1}),
|
||||
((0, 0), {(0, 1): 1}),
|
||||
((0, 0), {(2, 0): 1}),
|
||||
((0, 0), {(1, 1): 1}),
|
||||
((0, 0), {(3, 0): 1, (2, 1): QQ(3) / 2}),
|
||||
]
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
|
||||
|
||||
|
||||
def test_frobenius_basis_of_two_parameter_model_at_strong_coupling():
|
||||
# At (z0, z1) = (0, 1/4) the generators' z-free parts theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1)
|
||||
# vanish on the whole line theta1 = 0. Completing the indicial ideal pins the indicials down to (0, 0)
|
||||
# with multiplicity 4 and (0, 1/2) with multiplicity 2. The latter class holds sqrt(z1)*(1 + ...), the
|
||||
# power series solution of test_ideal_finds_power_series_solution_and_recovers_first_operator, and a
|
||||
# logarithmic partner sqrt(z1)*(log(z0) + log(z1) + ...).
|
||||
moved = _p22211_ideal().change_coordinates("z0", "z1 - 1/4")
|
||||
|
||||
assert moved.indicials() == [((0, 0), 4), ((0, QQ(1) / 2), 2)]
|
||||
basis = moved.frobenius_basis(order=4)
|
||||
|
||||
assert [period.indicials for period in basis] == [[0, 0]] * 4 + [[0, QQ(1) / 2]] * 2
|
||||
assert [_starting_term(period) for period in basis] == [
|
||||
((0, 0), {(0, 0): 1}),
|
||||
((0, 0), {(1, 0): 1}),
|
||||
((0, 0), {(2, 0): 1}),
|
||||
((0, 0), {(3, 0): 1}),
|
||||
((0, 0), {(0, 0): 1}),
|
||||
((0, 0), {(1, 0): 1, (0, 1): 1}),
|
||||
]
|
||||
power_series = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=4)
|
||||
assert basis[4].coefficients == power_series[0].coefficients
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
|
||||
|
||||
|
||||
def test_frobenius_basis_warns_when_log_degree_is_capped_too_low(caplog):
|
||||
# The triple-log solution at the MUM point needs log-degree 3, so with max_log_degree=2 only five of the
|
||||
# six solutions survive - which is reported rather than silently returned.
|
||||
ideal = _p22211_ideal()
|
||||
|
||||
with caplog.at_level(logging.WARNING):
|
||||
basis = ideal.frobenius_basis(order=2, max_log_degree=2)
|
||||
|
||||
assert len(basis) == 5
|
||||
assert any("rank sum(hodge_numbers) is 6" in record.getMessage() for record in caplog.records)
|
||||
|
||||
|
||||
def test_frobenius_basis_of_two_parameter_fourfold(caplog):
|
||||
# At the MUM point the nine solutions have the log structure [1, 2, 3, 2, 1], and the top ones are given by
|
||||
# the intersection numbers: the quartic one by kappa_ijkl*log(z_i)*...*log(z_l), i.e. 36*log(z0)^4 +
|
||||
# 48*log(z0)^3*log(z1) + 24*log(z0)^2*log(z1)^2 up to scale, a cubic one by its derivative in log(z1).
|
||||
ideal = _p333111_ideal()
|
||||
|
||||
assert ideal.indicials() == [((0, 0), 9)]
|
||||
basis = ideal.frobenius_basis(order=3)
|
||||
|
||||
assert [_starting_term(period) for period in basis] == [
|
||||
((0, 0), {(0, 0): 1}),
|
||||
((0, 0), {(1, 0): 1}),
|
||||
((0, 0), {(0, 1): 1}),
|
||||
((0, 0), {(2, 0): 1}),
|
||||
((0, 0), {(1, 1): 1}),
|
||||
((0, 0), {(0, 2): 1}),
|
||||
((0, 0), {(3, 0): 1, (2, 1): 1}),
|
||||
((0, 0), {(1, 2): 1}),
|
||||
((0, 0), {(4, 0): 1, (3, 1): QQ(4) / 3, (2, 2): QQ(2) / 3}),
|
||||
]
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
|
||||
|
||||
# At (0, -1/27), where the strong coupling divisor 1 + 27*z1 = 0 crosses the MUM divisor, all indicials lie
|
||||
# in one exponent class and there is a log(z0)^4 solution. But some indicials are (0, 1): this is no MUM
|
||||
# point, and its log structure [2, 3, 2, 1, 1] must not be compared with the Hodge numbers.
|
||||
moved = ideal.change_coordinates("z0", "z1 + 1/27")
|
||||
with caplog.at_level(logging.WARNING):
|
||||
basis = moved.frobenius_basis(order=3)
|
||||
|
||||
assert moved.indicials() == [((0, 0), 5), ((0, 1), 4)]
|
||||
assert len(basis) == 9
|
||||
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
|
||||
assert not caplog.records
|
||||
|
||||
|
||||
def test_frobenius_bases_at_mum_points_agree_with_the_classical_frobenius_method():
|
||||
# An independent construction of every period at the MUM points, logarithmic ones included: the classical
|
||||
# Frobenius method, from nothing but the Gamma-function coefficients c(k) of the fundamental period. The
|
||||
# starting terms are those of the canonical bases; for P_{22211}[8] the cubic one is
|
||||
# kappa_ijk*eps_i*eps_j*eps_k/8 with the intersection numbers kappa_000 = 8, kappa_001 = 4.
|
||||
eps_ring = PolynomialRing(QQ, 1, ["e0"])
|
||||
e = eps_ring.gen(0)
|
||||
|
||||
def quintic_ratio(k):
|
||||
# c(n + eps)/c(eps) for c(n) = (5n)!/(n!)^5
|
||||
ratio = eps_ring(1)
|
||||
for i in range(1, 5 * k[0] + 1):
|
||||
ratio = _truncated(ratio * (5 * e + i), 3)
|
||||
for i in range(1, k[0] + 1):
|
||||
ratio = _truncated(ratio * _inverse_of_linear(i, e, 3) ** 5, 3)
|
||||
return ratio
|
||||
|
||||
classical = _classical_frobenius_periods(quintic_ratio, [e**j for j in range(4)], order=8)
|
||||
basis = _quintic_ideal().frobenius_basis(order=8)
|
||||
assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
|
||||
|
||||
eps_ring = PolynomialRing(QQ, 2, ["e0", "e1"])
|
||||
e0, e1 = eps_ring.gens()
|
||||
|
||||
def p22211_ratio(k):
|
||||
# c(n + eps0, m + eps1)/c(eps0, eps1) for c(n, m) = (4n)!/((n!)^3 (m!)^2 (n - 2m)!). For n < 2m the
|
||||
# last factor is Gamma(1 + delta)/Gamma(1 + delta + n - 2m) with delta = eps0 - 2*eps1, a polynomial
|
||||
# vanishing at eps = 0, which only the logarithmic periods see.
|
||||
n, m = k
|
||||
delta = e0 - 2 * e1
|
||||
ratio = eps_ring(1)
|
||||
for i in range(1, 4 * n + 1):
|
||||
ratio = _truncated(ratio * (4 * e0 + i), 3)
|
||||
for i in range(1, n + 1):
|
||||
ratio = _truncated(ratio * _inverse_of_linear(i, e0, 3) ** 3, 3)
|
||||
for i in range(1, m + 1):
|
||||
ratio = _truncated(ratio * _inverse_of_linear(i, e1, 3) ** 2, 3)
|
||||
for i in range(1, n - 2 * m + 1):
|
||||
ratio = _truncated(ratio * _inverse_of_linear(i, delta, 3), 3)
|
||||
for i in range(2 * m - n):
|
||||
ratio = _truncated(ratio * (delta - i), 3)
|
||||
return ratio
|
||||
|
||||
starting_terms = [eps_ring(1), e0, e1, e0**2, e0 * e1, e0**3 + QQ(3) / 2 * e0**2 * e1]
|
||||
classical = _classical_frobenius_periods(p22211_ratio, starting_terms, order=6)
|
||||
basis = _p22211_ideal().frobenius_basis(order=6)
|
||||
assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
|
||||
|
||||
|
||||
def test_quintic_bases_at_mum_and_conifold_point_continue_into_each_other():
|
||||
# Both expansions converge on the lens between z0 = 0 and z0 = 1/3125, where the MUM periods must be
|
||||
# constant linear combinations T of the conifold periods. Through T the conifold monodromy
|
||||
# (log w -> log w + 2*pi*i) becomes a Picard-Lefschetz transvection in the MUM basis, and its product with
|
||||
# the MUM monodromy is the monodromy around z0 = oo, of order 5.
|
||||
CF = ComplexField(200)
|
||||
two_pi_i = 2 * CF.pi() * CF.gen()
|
||||
ideal = _quintic_ideal()
|
||||
mum = ideal.frobenius_basis(order=100)
|
||||
conifold = ideal.change_coordinates("1 - 3125*z0").frobenius_basis(order=100)
|
||||
|
||||
def values(basis, points, shifted=False):
|
||||
# One row per point z0, one column per period, continued once around w = 0 if shifted.
|
||||
rows = []
|
||||
for z0 in points:
|
||||
w = z0 if basis is mum else 1 - 3125 * z0
|
||||
log_w = CF(w).log() + (two_pi_i if shifted else 0)
|
||||
rows.append([
|
||||
sum(
|
||||
CF(coeff) * log_w ** log_index[0] * ((z_index[0] + period.indicials[0]) * log_w).exp()
|
||||
for log_index, z_dict in period.coefficients.items()
|
||||
for z_index, coeff in z_dict.items()
|
||||
if coeff != 0
|
||||
)
|
||||
for period in basis
|
||||
])
|
||||
return matrix(CF, rows)
|
||||
|
||||
fit = [CF(x, y) / 3125 for x, y in [(0.4, -0.15), (0.5, 0.1), (0.6, -0.1), (0.45, 0.15)]]
|
||||
test = [CF(x, y) / 3125 for x, y in [(0.55, 0.05), (0.42, -0.05), (0.6, 0.12)]]
|
||||
T = values(conifold, fit).solve_right(values(mum, fit)).transpose()
|
||||
assert (values(mum, test) - values(conifold, test) * T.transpose()).norm() < 1e-8 * values(mum, test).norm()
|
||||
|
||||
M_0 = values(mum, fit).solve_right(values(mum, fit, shifted=True)).transpose()
|
||||
M_c = T * values(conifold, fit).solve_right(values(conifold, fit, shifted=True)).transpose() * T.inverse()
|
||||
identity = identity_matrix(CF, 4)
|
||||
singular_values = (M_c - identity).change_ring(CDF).SVD()[1].diagonal()
|
||||
assert singular_values[1] < 1e-8 * singular_values[0]
|
||||
assert ((M_0 * M_c) ** 5 - identity).norm() < 1e-8
|
||||
@@ -1,22 +1,22 @@
|
||||
import os
|
||||
|
||||
import pytest
|
||||
from sage.all import * # noqa: F401
|
||||
from sage.all import * # noqa: F401,F403
|
||||
from sage.geometry.lattice_polytope import set_palp_dimension
|
||||
|
||||
load(os.path.join(os.path.dirname(__file__), "..", "sage", "toric_topdata.sage"))
|
||||
load(os.path.join(os.path.dirname(__file__), "test_models.sage"))
|
||||
load(os.path.join(os.path.dirname(__file__), "playground_models.sage"))
|
||||
|
||||
set_palp_dimension(11)
|
||||
|
||||
TEST_MODELS_PARAMETRISED = [pytest.param(*row[1:], id=row[0]) for row in TEST_MODELS]
|
||||
PLAYGROUND_MODELS = [pytest.param(*row[1:], id=row[0]) for row in PLAYGROUND_MODEL_FACTORIES] # noqa: F821
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"make_model, expected_cy_dimension, expected_no_divs, expected_disc, expected_intersection_numbers_CY",
|
||||
TEST_MODELS_PARAMETRISED,
|
||||
PLAYGROUND_MODELS,
|
||||
)
|
||||
def test_topdata_and_disc(
|
||||
def test_playground_models_disc_and_topdata(
|
||||
make_model,
|
||||
expected_cy_dimension,
|
||||
expected_no_divs,
|
||||
Reference in New Issue
Block a user