Implementing moving origin of ideal (#19)

---------

Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com>
Reviewed-on: #19

.gitea/ci-image/Dockerfile

@@ -1,26 +1,16 @@
-# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner.
+# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner,
+# since the official sagemath/sagemath image is amd64-only.
 #
-# The official sagemath/sagemath image is amd64-only; conda-forge's `sage`
-# package is the one place SageMath is actually published for linux-aarch64,
-# so this bakes it (plus the lint/test tools) into a single image, built and
-# pushed once rather than reinstalled on every CI run.
-#
-# Build & push (run directly on the Pi, or any arm64 machine with Docker):
+# Build & push commands:
 #   docker build -t gitea.piribauer.ch/julian/sage-ci:latest -f .gitea/ci-image/Dockerfile .
 #   docker login gitea.piribauer.ch -u julian
 #   docker push gitea.piribauer.ch/julian/sage-ci:latest
-#
-# Rebuild and re-push whenever this Dockerfile changes (e.g. bumping the sage
-# version) or the pinned tool versions need updating.
 
 FROM condaforge/miniforge3:latest
 
 RUN mamba install -y -c conda-forge sage ruff pytest \
     && mamba clean -afy
 
-# actions/checkout@v4 (and other JS-based actions) run via `node` inside this
-# container -- act_runner doesn't inject a runtime of its own when a custom
-# `container:` image is set, so one has to be present here.
 RUN mamba install -y -c conda-forge nodejs \
     && mamba clean -afy
 

.gitea/preparse.sh

@@ -1,18 +1,7 @@
 #!/usr/bin/env bash
-# Reproduces the CI "Preparse .sage files" step locally, so `ruff check` sees
-# the same generated *.sage.py files that the pipeline lints. Run this before
-# `ruff check --no-respect-gitignore .` to catch issues that only show up in
-# the generated output.
-#
-# Uses the sage.repl.preparse.preparse_file() Python API directly rather than
-# the `sage --preparse` CLI flag: the conda-forge `sage` binary used in CI is
-# a cut-down entry point that doesn't support `--preparse` (or `--python`) at
-# all, unlike the official sagemath CLI. Whichever `python` on PATH can
-# actually `import sage.repl` is used to run it -- in CI that's the conda
-# env's own `python`; locally (official sage install) it's `sage --python`.
-#
-# The generated files are gitignored; clean them up afterwards with:
-#   git clean -x sage/ playground/
+
+# This script is run in the CI container to preparse all .sage files into .sage.py files
+# and can also be run locally for local linting/testing.
 
 set -e
 

.gitea/workflows/ci.yml

@@ -9,9 +9,7 @@ on:
 jobs:
   lint:
     runs-on: ubuntu-latest
-    # Custom image (see .gitea/ci-image/Dockerfile): the official
-    # sagemath/sagemath image is amd64-only and this runner is arm64, and
-    # ruff/pytest are baked in here so jobs don't reinstall them every run.
+    # Custom image (see .gitea/ci-image/Dockerfile)
     container:
       image: gitea.piribauer.ch/julian/sage-ci:latest
 
@@ -20,18 +18,9 @@ jobs:
         uses: actions/checkout@v4
 
       - name: Preparse .sage files
-        # Turns each *.sage file into real Python (*.sage.py) so ruff can
-        # parse it, and prepends the `from sage.all import *` that `sage`
-        # normally injects at runtime, so Sage's globals (ZZ, var, matrix, ...)
-        # resolve instead of looking like undefined names.
         run: .gitea/preparse.sh
 
       - name: ruff check
-        # --no-respect-gitignore: *.sage.py is gitignored (it's generated, see
-        # the previous step) but that's exactly what we need to lint here.
-        # `python` here is the conda env's own interpreter -- it already has
-        # sage/ruff/pytest importable, `sage --python` isn't a real flag on
-        # the conda-forge sage CLI.
         run: python -m ruff check --no-respect-gitignore .
 
   test:

README.md

@@ -1,2 +1,151 @@
 # 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([1, 2, 3], 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
+
+## period_computation
+
+`period_computation.sage` provides classes for working with Picard–Fuchs operators and their
+period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and
+last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given
+z- and theta-degree.
+
+A `PFOperator` is parsed from a string in the variables `z0, ..., z<n-1>` and `theta0, ..., theta<n-1>`,
+the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard&ndash;Fuchs operator:
+
+```python
+L = PFOperator(
+    "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0",
+    no_variables=1,
+)
+
+L.simplify().operator_string
+```
+
+```term
+'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4'
+```
+
+An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial
+exponents and order.
+
+```python
+ideal = PFIdeal([L])
+period = ideal.find_power_series_solution(indicials=[0], order=3)[0]
+
+period.period_string
+```
+
+```term
+'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1'
+```
+
+The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for
+`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown
+coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter
+model P_{2,2,2,1,1}[8]:
+
+```python
+M1 = PFOperator(
+    "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1",
+    no_variables=2,
+)
+M2 = PFOperator(
+    "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1",
+    no_variables=2,
+)
+ideal = PFIdeal([M1, M2])
+
+period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0]
+recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
+
+period.period_string
+recovered[0].simplify().operator_string
+```
+
+```term
+'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)'
+'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1'
+```
+
+`recovered[0]` is, up to scale, `M1` &mdash; recovered purely from `M1`, `M2`'s shared power series
+solution. 

pyproject.toml

@@ -3,33 +3,18 @@ line-length = 120
 target-version = "py311"
 
 [tool.ruff.lint]
+# E: pycodestyle errors, F: pyflakes (undefined/unused names), I: isort (import order), W: pycodestyle warnings
 select = ["E", "F", "I", "W"]
 ignore = [
     "E501", # symbolic-math expressions routinely exceed a "normal" line length
 ]
 
 [tool.ruff.lint.per-file-ignores]
-# *.sage.py is generated by `sage --preparse` from *.sage sources.
-# - F403/F405/F821/E741: Sage's runtime injects hundreds of globals (ZZ, var,
-#   matrix, LatticePolytope, ...) via `from sage.all import *` before
-#   executing these files, so plain static analysis can't see where names
-#   come from.
-# - E402/E702/I001: Sage's preparser itself emits an import-then-semicolon-
-#   joined-constants preamble (e.g. `_sage_const_1 = Integer(1); ...`) ahead
-#   of the file's own imports; that's Sage's boilerplate, not this project's
-#   code style.
-# - W291/W293: the preparser replaces every integer literal with
-#   `_sage_const_N ` (trailing space included) to preserve token boundaries,
-#   so trailing-whitespace warnings fire mechanically on almost every line
-#   with a number in it -- not something `ruff format` could fix anyway,
-#   since it wouldn't change the .sage source that generated it.
-# - W292: whether the reconstructed file ends in a real trailing newline
-#   depends on how many blank lines the preparser appends, which varies by
-#   Sage version -- not something worth pinning tool versions over.
+# *.sage.py is preparser output: star-import globals, semicolon preamble, and literal-substitution artifacts trip static analysis.
 "*.sage.py" = ["F403", "F405", "F821", "E741", "E402", "E702", "I001", "W291", "W293", "W292"]
 
-# tests/test_smoke.py does `from sage.all import *` and `load(...)` a .sage
-# file to get at its classes -- the same dynamic-namespace situation as
-# *.sage.py above, just in a hand-written file: ruff can't see that
-# Polytope/ToricPolytope/etc. come from the loaded file.
-"tests/test_smoke.py" = ["F403", "F405"]
+# test_smoke.py star-imports sage.all and loads a .sage file, so ruff can't see where its names come from.
+"tests/test_topdata_and_disc.py" = ["F403", "F405"]
+
+# Same star-import + load() pattern as above.
+"tests/test_period_computation.py" = ["F403", "F405"]

sage/period_computation.sage

@@ -0,0 +1,646 @@
+import copy
+import logging
+from sage.all import sage_eval, PolynomialRing, QQ, SR, 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
+
+class PFOperatorAnsatz(PFOperator):
+    """
+    A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees.
+    """
+    def __init__(self, no_variables: int, theta_degree: int, z_degree: int):
+        self.no_variables = no_variables
+        self.theta_degree = theta_degree
+        self.z_degree = z_degree
+
+        z_indices = _multi_indices(self.no_variables, z_degree)
+        theta_indices = _multi_indices(self.no_variables, theta_degree)
+
+        # Every (z_index, theta_index) pair allowed by the degree bounds gets its own
+        # fresh unknown, to be solved for once the Ansatz is applied to a period.
+        operator_terms = [
+            (var("b_" + "_".join(str(x) for x in z_index + theta_index)), (z_index, theta_index))
+            for z_index in z_indices
+            for theta_index in theta_indices
+        ]
+        self.unknowns = [coeff for coeff, _ in operator_terms]
+
+        def _monomial_factors(index, name):
+            return [f"{name}{i}^{exp}" for i, exp in enumerate(index) if exp > 0]
+
+        operator_string = " + ".join(
+            "*".join([str(coeff), *_monomial_factors(z_index, "z"), *_monomial_factors(theta_index, "theta")])
+            for coeff, (z_index, theta_index) in operator_terms
+        )
+
+        extra_locals = {str(coeff): coeff for coeff in self.unknowns}
+        super().__init__(operator_string=operator_string, no_variables=no_variables, extra_locals=extra_locals)
+
+
+class PFIdeal:
+    """
+    A class representing a Picard-Fuchs ideal, which is a collection of PFOperators.
+    """
+    def __init__(self, operators: list):
+        if not all(op.no_variables == operators[0].no_variables for op in operators):
+            raise ValueError("All operators must have the same number of variables.")
+        self.no_variables = operators[0].no_variables
+        self.operators = operators
+        logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators))
+    
+    def add_operator(self, pf_operator: PFOperator):
+        self.operators.append(pf_operator)
+        logger.info("Added PFOperator to PFIdeal: %s", pf_operator.operator_string)
+    
+    def remove_operator(self, pf_operator: PFOperator):
+        self.operators.remove(pf_operator)
+        logger.info("Removed PFOperator from PFIdeal: %s", pf_operator.operator_string)
+
+    def change_coordinates(self, *new_coordinates: str) -> "PFIdeal":
+        """
+        Moves the origin of the whole system of differential equations to a new point, by
+        returning a new PFIdeal whose operators are each rewritten in new_coordinates; see
+        PFOperator.change_coordinates for what new_coordinates means.
+        """
+        new_operators = [pf_operator.change_coordinates(*new_coordinates) for pf_operator in self.operators]
+        logger.info("Changed coordinates of PFIdeal with %d operator(s) to %s.", len(new_operators), new_coordinates)
+        return PFIdeal(new_operators)
+
+    def find_power_series_solution(self, indicials: list, order: int) -> list:
+        """
+        Finds power series solutions (no logs) to this PFIdeal at given indicial exponents/order.
+
+        A term c*z^p*theta^q sends a_k*z^k (true exponent k+indicials) to
+        c*(k+indicials)^q*a_k*z^(k+p): it shifts index k up by p (p>=0), never down or sideways.
+        E.g. z0*theta0 sends a_k*z0^k to (k+rho0)*a_k*z0^(k+1).
+
+        This is the multivariate Frobenius method: canonical series solutions of a regular
+        holonomic D-ideal via its indicial ideal. See M. Saito, B. Sturmfels, N. Takayama,
+        "Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and
+        Computation in Mathematics vol. 6, Springer, 2000, chs. 2-3.
+        """
+        if len(indicials) != self.no_variables:
+            raise ValueError(
+                "Indicials must have length no_variables=%d, got %d." % (self.no_variables, len(indicials))
+            )
+        indicials = [QQ(rho) for rho in indicials]
+
+        operator_terms = [
+            [
+                (monomial.exponents()[0][:self.no_variables], monomial.exponents()[0][self.no_variables:], coeff)
+                for coeff, monomial in pf_operator.operator
+            ]
+            for pf_operator in self.operators
+        ]
+
+        def eigenvalue(k, q):
+            # theta_i^q_i acts on z_i^(k_i + indicials[i]) as multiplication by
+            # (k_i + indicials[i])^q_i; theta^q's combined eigenvalue is the product over i.
+            value = QQ(1)
+            for i in range(self.no_variables):
+                if q[i]:
+                    value *= (k[i] + indicials[i]) ** q[i]
+            return value
+
+        # Process multi-indices in order of increasing total degree: since every operator
+        # monomial has p >= 0 (componentwise), the z^m coefficient of L(y) only ever
+        # depends on a_k for k <= m, so by this point every k < m has already been solved.
+        z_indices = sorted(_multi_indices(self.no_variables, order), key=sum)
+
+        # solved[k] holds a_k written as a vector of coefficients over the `dimension`
+        # independent solutions found so far (a basis of the solution space up to k).
+        dimension = 0
+        solved = {}
+
+        for m in z_indices:
+            # For each operator, "coefficient of z^m in L(y) = 0" splits into a diagonal
+            # part (the p=0, theta-only monomials, whose unknown is a_m itself) plus a
+            # known part contributed by already-solved a_k with k = m - p, p > 0.
+            diagonals = []
+            known_contributions = []
+            for terms in operator_terms:
+                diagonal = QQ(0)
+                contribution = [QQ(0)] * dimension
+                for p, q, coeff in terms:
+                    k = tuple(m[i] - p[i] for i in range(self.no_variables))
+                    if any(ki < 0 for ki in k):
+                        continue  # this monomial would need a_k for a negative multi-index k: no such term
+                    ev = eigenvalue(k, q)
+                    if ev == 0:
+                        continue  # theta^q kills z_i^(k_i + indicials[i]) here, so this monomial contributes nothing
+                    if k == m:
+                        diagonal += coeff * ev  # p = 0: coefficient multiplying the still-unknown a_m
+                    else:
+                        k_vector = solved[k]  # p > 0: a_k is already known, add its contribution
+                        for i in range(dimension):
+                            contribution[i] += coeff * ev * k_vector[i]
+                diagonals.append(diagonal)
+                known_contributions.append(contribution)
+
+            # An operator with a nonzero diagonal lets us solve a_m = -(known part)/diagonal
+            # directly; this is exactly the indicial equation being nonzero at m + indicials.
+            active = next((r for r, d in enumerate(diagonals) if d != 0), None)
+
+            if active is not None:
+                value = [-known_contributions[active][i] / diagonals[active] for i in range(dimension)]
+                # Every operator's equation at m must independently be satisfied by this
+                # same a_m; disagreement means the ideal is inconsistent with these indicials.
+                for r, d in enumerate(diagonals):
+                    if any(d * value[i] + known_contributions[r][i] != 0 for i in range(dimension)):
+                        raise ValueError(
+                            "Inconsistent Picard-Fuchs ideal or indicial exponents %s at multidegree %s."
+                            % (indicials, m)
+                        )
+                solved[m] = value
+            else:
+                # Resonance: every operator's indicial part vanishes at m + indicials, so a_m
+                # cannot be pinned down by this equation. If the already-known lower-degree
+                # data still forces a nonzero constraint here, satisfying it would require a
+                # log(z)-term solution, which this method (deliberately) does not compute.
+                if any(x != 0 for contribution in known_contributions for x in contribution):
+                    raise ValueError(
+                        "Resonance at multidegree %s for indicials %s would require a logarithmic "
+                        "solution, which find_power_series_solution does not compute." % (m, indicials)
+                    )
+                # Otherwise a_m is genuinely free: it starts a new independent solution, so
+                # extend every previously solved coefficient with a 0 in this new direction.
+                dimension += 1
+                for v in solved.values():
+                    v.append(QQ(0))
+                solved[m] = [QQ(0)] * (dimension - 1) + [QQ(1)]
+
+        log_index = tuple([0] * self.no_variables)
+        solutions = [
+            Period(
+                no_variables=self.no_variables,
+                coefficients={log_index: {m: solved[m][i] for m in z_indices}},
+                order=order,
+                indicials=indicials,
+            )
+            for i in range(dimension)
+        ]
+
+        logger.info(
+            "Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order
+        )
+        return solutions
+
+
+class Period:
+    """
+    A class representing a period as a formal power series in z-variables and their logs.
+
+    The coefficients are stored in a dictionary of dictionaries, where the first key is the multi-index
+    of the logarithmic part and the second key is the multi-index of the z-variables. 
+    So, for example, the coefficient of
+
+        log(z0)^2 * log(z1) * z0^3 * z1^2 would be stored as coefficients[(2, 1)][(3, 2)]
+
+    """
+
+    def _initialise_ring(self):
+        z_names = ['z%d' % i for i in range(self.no_variables)]
+        log_names = ['L%d' % i for i in range(self.no_variables)]
+        ring = PolynomialRing(QQ, z_names + log_names)
+        z_gens = ring.gens()[:self.no_variables]
+        log_gens = ring.gens()[self.no_variables:]
+
+        return ring, z_gens, log_gens  # type: (PolynomialRing, list, list)
+
+    def _period_from_string(self, period_string: str) -> dict:
+        ring, z_gens, log_gens = self._initialise_ring()
+        log_of_z = dict(zip(z_gens, log_gens))
+
+        def log(zi):
+            try:
+                return log_of_z[zi]
+            except (KeyError, TypeError):
+                raise ValueError(
+                    "log(...) may only be applied to one of the z-variables z0, ..., z%d"
+                    % (self.no_variables - 1)
+                )
+
+        parse_locals = dict(ring.gens_dict())
+        parse_locals['log'] = log
+
+        try:
+            expression = ring(sage_eval(period_string, locals=parse_locals))
+        except Exception as e:
+            raise ValueError("Invalid period string: %s" % e)
+
+        coefficients = {}
+        for coeff, monomial in expression:
+            exponents = monomial.exponents()[0]
+            z_index = tuple(exponents[:self.no_variables])
+            log_index = tuple(exponents[self.no_variables:])
+            coefficients.setdefault(log_index, {})[z_index] = coeff
+
+        return coefficients
+    
+    def _period_to_string(self) -> str:
+        ring, z_gens, log_gens = self._initialise_ring()
+        expression = ring(0)
+        for log_index, z_dict in self.coefficients.items():
+            for z_index, coeff in z_dict.items():
+                monomial = coeff
+                for i in range(self.no_variables):
+                    monomial *= (log_gens[i] ** log_index[i]) * (z_gens[i] ** z_index[i])
+                expression += monomial
+
+        log_substitutions = {log_gens[i]: log(SR(z_gens[i])) for i in range(self.no_variables)}
+        result = SR(expression).subs(log_substitutions)
+
+        # A nonzero indicial ρ_i means the coefficients above are for z_i^k, but the
+        # actual solution is z_i^(ρ_i + k); make that explicit in the printed form.
+        indicial_prefactor = prod(
+            (SR(z_gens[i]) ** self.indicials[i] for i in range(self.no_variables) if self.indicials[i] != 0),
+            SR(1),
+        )
+        if indicial_prefactor != 1:
+            result *= indicial_prefactor
+
+        return str(result)
+
+    def _max_z_degree(self) -> int:
+        # Highest total z-degree (sum of the z-multi-index) among all coefficients.
+        if not self.coefficients:
+            return 0
+        return max(
+            sum(z_index)
+            for z_dict in self.coefficients.values()
+            for z_index in z_dict
+        )
+
+    def _truncate_coefficients(self, order: int) -> dict:
+        # Drop every (log_index, z_index) entry whose total z-degree exceeds order,
+        # removing it from the dictionary rather than merely zeroing it out.
+        truncated = {}
+        for log_index, z_dict in self.coefficients.items():
+            kept = {z_index: coeff for z_index, coeff in z_dict.items() if sum(z_index) <= order}
+            if kept:
+                truncated[log_index] = kept
+        return truncated
+
+    def _apply_theta(self, coefficients: dict, index: int) -> dict:
+        """
+        Applies the logarithmic derivative theta_i = z_i * d/dz_i once to a coefficients dict of the same
+        shape as self.coefficients. It uses the product rule
+
+            theta_i(z^a log(z)^k) = a_i * z^a log(z)^k + k_i * z^a log(z)^(k - e_i).
+
+        z_index entries are offsets from self.indicials: a stored z_index of a really means
+        z^(a + self.indicials[index]), so theta_i's eigenvalue is a_i + self.indicials[index]
+        rather than the bare a_i (self.indicials is all-zero unless explicitly given, in which
+        case this reduces to the ordinary power-series rule).
+        """
+        result = {}
+        for log_index, z_dict in coefficients.items():
+            for z_index, coeff in z_dict.items():
+                a_i = z_index[index] + self.indicials[index]
+                if a_i != 0:
+                    inner = result.setdefault(log_index, {})
+                    inner[z_index] = inner.get(z_index, 0) + a_i * coeff
+
+                k_i = log_index[index]
+                if k_i != 0:
+                    lowered_log_index = log_index[:index] + (k_i - 1,) + log_index[index + 1:]
+                    inner = result.setdefault(lowered_log_index, {})
+                    inner[z_index] = inner.get(z_index, 0) + k_i * coeff
+        return result
+
+    def apply_operator(self, pf_operator: PFOperator) -> "Period":
+        """
+        Applies a PFOperator to this period and returns the result as a new Period,
+        truncated to self.order. Each monomial of the operator is normalised as
+        coeff * z^p * theta^q (see PFOperator's docstring), so theta^q is applied
+        to the period first and the result is then multiplied by coeff * z^p.
+        """
+        if pf_operator.no_variables != self.no_variables:
+            raise ValueError(
+                "Variable count mismatch: period has %d variable(s), operator has %d."
+                % (self.no_variables, pf_operator.no_variables)
+            )
+
+        result_coefficients = {}
+        for monomial_coeff, monomial in pf_operator.operator:
+            exponents = monomial.exponents()[0]
+            z_exponents = exponents[:self.no_variables]
+            theta_exponents = exponents[self.no_variables:]
+
+            term = self.coefficients
+            for i, power in enumerate(theta_exponents):
+                for _ in range(power):
+                    term = self._apply_theta(term, i)
+
+            for log_index, z_dict in term.items():
+                inner = result_coefficients.setdefault(log_index, {})
+                for z_index, coeff in z_dict.items():
+                    shifted_z_index = tuple(z_index[i] + z_exponents[i] for i in range(self.no_variables))
+                    inner[shifted_z_index] = inner.get(shifted_z_index, 0) + monomial_coeff * coeff
+
+        result_coefficients = {
+            log_index: {z_index: c for z_index, c in z_dict.items() if c != 0}
+            for log_index, z_dict in result_coefficients.items()
+        }
+        result_coefficients = {log_index: z_dict for log_index, z_dict in result_coefficients.items() if z_dict}
+
+        return Period(
+            no_variables=self.no_variables,
+            coefficients=result_coefficients,
+            order=self.order,
+            indicials=self.indicials,
+        )
+
+    def find_annihilating_operators(self, z_degree: int, theta_degree: int) -> list:
+        """
+        Finds operators annihilating this period among PFOperator Ansätze of the given
+        z- and theta-degree.
+        """
+        ansatz = PFOperatorAnsatz(self.no_variables, theta_degree, z_degree)
+        unknowns = ansatz.unknowns
+
+        equations = [
+            coeff
+            for z_dict in self.apply_operator(ansatz).coefficients.values()
+            for coeff in z_dict.values()
+        ]
+
+        coefficient_rows = [[SR(equation).coefficient(b) for b in unknowns] for equation in equations]
+        kernel_basis = matrix(QQ, coefficient_rows, ncols=len(unknowns)).right_kernel().basis()
+
+        operators = []
+        for basis_vector in kernel_basis:
+            solution = {unknowns[j]: basis_vector[j] for j in range(len(unknowns))}
+            solved_operator = ansatz.operator.map_coefficients(lambda c: SR(c).subs(solution))
+            operators.append(PFOperator(str(solved_operator), no_variables=self.no_variables))
+
+        return operators
+
+    def __init__(
+        self,
+        no_variables: int = 1,
+        coefficients: dict = None,
+        period_string: str = None,
+        order: int = None,
+        indicials: list = None,
+    ):
+        self.no_variables = no_variables
+        # indicials[i] is the (rational) Frobenius exponent ρ_i of z_i: a stored
+        # z-index of a really represents z^(a + indicials[i]).
+        self.indicials = list(indicials) if indicials is not None else [0] * no_variables
+
+        # Use coefficients or period_string to initialize the period
+        if period_string is not None:
+            if coefficients is not None:
+                raise ValueError("Provide either coefficients or period_string, not both.")
+            self.coefficients = self._period_from_string(period_string)
+            self.period_string = period_string
+            logger.debug("Using string for initialisation.")
+        elif coefficients is None:
+            self.coefficients = {}
+        else:
+            self.coefficients = coefficients
+            self.period_string = self._period_to_string()
+            logger.debug("Using coefficients for initialisation.")
+
+        if order is None:
+            self.order = self._max_z_degree()
+            logger.debug("Order not provided, using maximum z-degree: %d", self.order)
+        else:
+            self.coefficients = self._truncate_coefficients(order)
+            self.order = order
+            self.period_string = self._period_to_string()
+            logger.debug("Truncated coefficients to order %d.", self.order)
+
+        logger.info("Initialised Period in %d variables at order %d.", self.no_variables, self.order)
+
+
+class PeriodAnsatz(Period):
+    """
+    A class for period Ansätze, characterised by number of variables and their z- and log-multi-degrees.
+    """
+
+    def __init__(self, no_variables: int, z_degree: int, log_degree: int, indicials: list = None):
+        self.no_variables = no_variables
+        self.z_degree = z_degree
+        self.log_degree = log_degree
+
+        log_indices = _multi_indices(self.no_variables, log_degree)
+        z_indices = _multi_indices(self.no_variables, z_degree)
+
+        # Every (log_index, z_index) pair allowed by the degree bounds gets its own
+        # fresh unknown, to be solved for once a PFOperator is applied to the Ansatz.
+        coefficients = {
+            log_index: {
+                z_index: var("a_" + "_".join(str(x) for x in log_index + z_index))
+                for z_index in z_indices
+            }
+            for log_index in log_indices
+        }
+
+        super().__init__(no_variables=no_variables, coefficients=coefficients, order=z_degree, indicials=indicials)
+        self.expansion_coefficients = self.coefficients
+        logger.info(
+            "Initialised PeriodAnsatz with %d unknown coefficient(s).",
+            sum(len(z_dict) for z_dict in self.expansion_coefficients.values()),
+        )

sage/toric_topdata.sage

@@ -101,7 +101,8 @@ 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))
         
-        logger.info("Using triangulation index %d of [0..%d].", 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]
     
     def set_nef_partition(self, nef_partition):
@@ -739,3 +740,34 @@ class ToricPolytopeCICY(ToricPolytope):
             lvec=lvec
         )
     
+class ToricPolytopeProjectiveSpace(ToricPolytope):
+    """
+    Computes topological data for hypersurfaces in toric ambient spaces.
+    """
+    def _compute_points_from_weights(self, weights):
+        # Convert the weights of the projective space to the corresponding points
+        dimension = len(weights) - 1
+        points = identity_matrix(dimension)
+        try:
+            weights.remove(1)
+        except ValueError:
+            raise ValueError("The weights must include a 1 for the projective space.")
+
+        points = points.insert_row(0, [-w for w in weights])
+        return list(points)
+    
+    def __init__(
+            self,
+            weights,
+            no_triangulation = 0,
+            model_name = None,
+            nef_partition = None,
+            lvec = None, # Note that the inner point is the last entry
+        ):
+        super().__init__(
+            self._compute_points_from_weights(weights),
+            no_triangulation=no_triangulation,
+            model_name=model_name,
+            nef_partition=nef_partition,
+            lvec=lvec
+        )

sage/util.py

@@ -0,0 +1,15 @@
+"""
+Utility functions for the Calabi-Yau period geometry project.
+"""
+
+def _multi_indices(no_variables, total_degree: int) -> list:
+    # All no_variables-tuples of non-negative integers whose entries sum to at most total_degree.
+    def helper(remaining_vars, remaining_degree):
+        if remaining_vars == 0:
+            yield ()
+            return
+        for first in range(remaining_degree + 1):
+            for rest in helper(remaining_vars - 1, remaining_degree - first):
+                yield (first,) + rest
+
+    return list(helper(no_variables, total_degree))

tests/test_models.sage

@@ -0,0 +1,104 @@
+# Entries are (id, factory, expected_cy_dimension, expected_no_divs, expected_disc,
+# expected_intersection_numbers_CY), where expected_disc is a list of
+# (str(discriminant_factor), codimension) pairs as returned by ToricPolytope.disc().
+
+TEST_MODELS = [
+    (
+        "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)],
+        {(0, 0): 4, (0, 1): 2, (1, 1): 0},
+    ),
+    (
+        "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)],
+        {(0, 0, 0): 8, (0, 0, 1): 4, (0, 1, 1): 0, (1, 1, 1): 0},
+    ),
+    (
+        "CY4_two_parameter_family",
+        lambda: ToricPolytopeProjectiveSpace([1, 1, 1, 1, 8, 12], model_name="CY4_two_parameter_family"),
+        4,
+        2,
+        [
+            ("z2 - 1/256", 2),
+            ("34828517376*z1^4*(256*z2 - 1) + 322486272*z1^3 - 1119744*z1^2 + 1728*z1 - 1", 0),
+        ],
+        {(0, 0, 0, 0): 64, (0, 0, 0, 1): 16, (0, 0, 1, 1): 4, (0, 1, 1, 1): 1, (1, 1, 1, 1): 0},
+    ),
+    (
+        "CICY3_one_parameter",
+        lambda: ToricPolytopeCICY([[3, 3]], model_name="CICY3_one_parameter"),
+        3,
+        1,
+        [("z1 - 1/46656", 0)],
+        {(0, 0, 0): 9},
+    ),
+    (
+        "CICY3_two_parameter",
+        lambda: ToricPolytopeCICY([[3], [3]], model_name="CICY3_two_parameter"),
+        3,
+        2,
+        [
+            (
+                "-19683*z1^3 - 2187*(27*z1 + 1)*z2^2 - 19683*z2^3 - 2187*z1^2 "
+                "- 81*(729*z1^2 - 189*z1 + 1)*z2 - 81*z1 - 1",
+                0,
+            )
+        ],
+        {(0, 0, 0): 0, (0, 0, 1): 3, (0, 1, 1): 3, (1, 1, 1): 0},
+    ),
+    (
+        "CICY5_two_parameter",
+        lambda: ToricPolytopeCICY([[6, 1], [0, 2]], model_name="CICY5_two_parameter"),
+        5,
+        2,
+        [
+            (
+                "16384*z2^7 - 28672*z2^6 + 21504*z2^5 - 448*(1647086*z1 - 5)*z2^3 - 8960*z2^4 "
+                "- 112*(8235430*z1 + 3)*z2^2 - 678223072849*z1^2 - 28*(4941258*z1 - 1)*z2 - 1647086*z1 - 1",
+                0,
+            )
+        ],
+        {
+            (0, 0, 0, 0, 0): 6,
+            (0, 0, 0, 0, 1): 12,
+            (0, 0, 0, 1, 1): 0,
+            (0, 0, 1, 1, 1): 0,
+            (0, 1, 1, 1, 1): 0,
+            (1, 1, 1, 1, 1): 0,
+        },
+    ),
+    (
+        "CICY3_two_parameter_manual",
+        lambda: 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]],
+            model_name="CICY3_two_parameter_manual",
+        ),
+        4,
+        2,
+        [
+            (
+                "-14348907*z1^5 - 2657205*z1^4 - 196830*z1^3 - 29296875*(270*z1 + 1)*z2^2 "
+                "- 30517578125*z2^3 - 7290*z1^2 + 9375*(98415*z1^3 - 32805*z1^2 + 810*z1 - 1)*z2 "
+                "- 135*z1 - 1",
+                0,
+            )
+        ],
+        {(0, 0, 0, 0): 0, (0, 0, 0, 1): 0, (0, 0, 1, 1): 6, (0, 1, 1, 1): 8, (1, 1, 1, 1): 2},
+    ),
+]

tests/test_period_computation.py

@@ -0,0 +1,203 @@
+import os
+
+from sage.all import *  # noqa: F401
+
+load(os.path.join(os.path.dirname(__file__), "..", "sage", "period_computation.sage"))
+
+
+def test_apply_operator_univariate():
+    # theta^2 - z applied to log(z) should give -z*log(z), since theta(log z) = 1
+    # and theta^2(log z) = theta(1) = 0.
+    period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5)
+    op = PFOperator("theta0^2 - z0", no_variables=1)
+
+    result = period.apply_operator(op)
+
+    assert result.coefficients == {(1,): {(1,): -1}}
+
+
+def test_apply_operator_mixes_variables():
+    # z0*theta1 applied to z0*log(z1): theta1 strips log(z1) down to a bare 1
+    # (leaving z0 untouched), then multiplying by z0 gives z0^2.
+    period = Period(no_variables=2, coefficients={(0, 1): {(1, 0): 1}}, order=5)
+    op = PFOperator("z0*theta1", no_variables=2)
+
+    result = period.apply_operator(op)
+
+    assert result.coefficients == {(0, 0): {(2, 0): 1}}
+
+
+def test_apply_operator_truncates_to_order():
+    # Multiplying by z0^2 pushes some terms above the period's order, so they
+    # should be dropped rather than kept with a nonzero coefficient.
+    period = Period(no_variables=1, coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1}}, order=2)
+    op = PFOperator("z0^2", no_variables=1)
+
+    result = period.apply_operator(op)
+
+    assert result.coefficients == {(0,): {(2,): 1}}
+    assert result.order == 2
+
+
+def test_find_annihilating_operators_recovers_theta_squared():
+    # theta0^2 annihilates log(z0): theta0(log z0) = 1, theta0^2(log z0) = theta0(1) = 0.
+    # Among degree-(0, 2) Ansatze this should be the only solution, up to scaling.
+    period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5)
+
+    ops = period.find_annihilating_operators(z_degree=0, theta_degree=2)
+
+    # A one-dimensional solution space means exactly one basis operator.
+    assert len(ops) == 1
+    assert period.apply_operator(ops[0]).coefficients == {}
+    assert str(ops[0].operator) == "theta0^2"
+
+
+def test_find_annihilating_operators_recovers_geometric_series_operator():
+    # sum_{k=0}^{4} z0^k is annihilated (up to truncation order) by
+    # (1 - z0)*theta0 - z0, i.e. theta0 - z0*theta0 - z0.
+    period = Period(
+        no_variables=1,
+        coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1, (3,): 1, (4,): 1}},
+        order=4,
+    )
+
+    ops = period.find_annihilating_operators(z_degree=1, theta_degree=1)
+
+    assert len(ops) == 1
+    assert period.apply_operator(ops[0]).coefficients == {}
+
+
+def test_find_annihilating_operators_returns_empty_list_when_no_solution_exists():
+    # No degree-0 (constant) operator other than the zero operator can annihilate a
+    # nonzero constant period, so the linear system's only solution is trivial.
+    period = Period(no_variables=1, coefficients={(0,): {(0,): 1}}, order=0)
+
+    ops = period.find_annihilating_operators(z_degree=0, theta_degree=0)
+
+    assert ops == []
+
+
+def test_simplify_factorises_theta_polynomial_per_z_monomial():
+    # theta0^4 - 5*z0*(5*theta0+1)*(5*theta0+2)*(5*theta0+3)*(5*theta0+4), expanded, is the
+    # quintic's Picard-Fuchs operator. simplify() should recover the factorised form: the
+    # z0^0 part (theta0^4) has no theta-factor to pull out, while the z0^1 part factorises
+    # into the four linear pieces.
+    expanded = "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0"
+    op = PFOperator(expanded, no_variables=1)
+
+    simplified = op.simplify()
+
+    # The underlying (expanded) operator is unchanged - only the display string differs.
+    assert simplified.operator == op.operator
+    assert simplified.operator_string == "-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4"
+
+
+def test_find_power_series_solution_recovers_quintic_period():
+    # The quintic's Picard-Fuchs operator has indicial equation theta0^4 = 0 at z0 = 0
+    # (a quadruple root at 0), so its holomorphic (non-logarithmic) power series solution
+    # is found at indicial 0. Up to normalisation this is the classic quintic period
+    # 1 + 120*z0 + 113400*z0^2 + ...
+    op = PFOperator(
+        "theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)",
+        no_variables=1,
+    )
+    ideal = PFIdeal([op])
+
+    solutions = ideal.find_power_series_solution(indicials=[0], order=3)
+
+    assert len(solutions) == 1
+    assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000}}
+
+
+def test_find_power_series_solution_handles_rational_indicial():
+    # theta0 - 1/2 kills z0^(1/2 + k) only when k = 0, since its eigenvalue is 1/2 + k;
+    # so at indicial 1/2 there is exactly one solution (a constant multiple of sqrt(z0)),
+    # while at indicial 0 no power series solution exists at all.
+    op = PFOperator("theta0 - 1/2", no_variables=1)
+    ideal = PFIdeal([op])
+
+    solutions = ideal.find_power_series_solution(indicials=[1/2], order=3)
+    assert len(solutions) == 1
+    assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 0, (2,): 0, (3,): 0}}
+    assert solutions[0].period_string == "sqrt(z0)"
+
+    assert ideal.find_power_series_solution(indicials=[0], order=3) == []
+
+
+def test_ideal_finds_power_series_solution_and_recovers_first_operator():
+    # Picard--Fuchs ideal for P_{22211}[8] 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_change_coordinates_handles_non_monomial_denominator():
+    op = PFOperator("theta0^2 - z0", no_variables=1)
+
+    result = op.change_coordinates("z0/(1-z0)")
+
+    expected = PFOperator(
+        "(1+z0)^3*theta0^2 + z0*(1+z0)^2*theta0 - z0",
+        no_variables=1,
+    )
+    assert result.operator == expected.operator

tests/test_smoke.py

@@ -1,32 +0,0 @@
-import os
-
-import pytest
-from sage.all import *  # noqa: F401,F403
-
-load(os.path.join(os.path.dirname(__file__), "..", "sage", "toric_topdata.sage"))
-
-
-def test_polytope_rejects_too_few_points():
-    with pytest.raises(ValueError):
-        Polytope([[0, 0], [1, 0]])  # 2 points in dimension 2: not enough
-
-
-def test_mirror_quintic_topdata():
-    # Vertices e_1, ..., e_4, -e_1-...-e_4 -- the worked example from
-    # ToricPolytope's own docstring.
-    points = [
-        [1, 0, 0, 0],
-        [0, 1, 0, 0],
-        [0, 0, 1, 0],
-        [0, 0, 0, 1],
-        [-1, -1, -1, -1],
-    ]
-
-    polytope = ToricPolytope(points, model_name="ci_smoke_test_quintic")
-    assert polytope.dimension == 4
-    assert len(polytope.triangulations) >= 1
-
-    polytope.topdata()
-    assert polytope.cy_dimension == 3
-    assert polytope.no_divs == 1  # P^4 has Picard rank 1
-    assert polytope.intersection_numbers_CY == {(0, 0, 0): 5}  # classical quintic self-intersection

tests/test_topdata_and_disc.py

@@ -0,0 +1,34 @@
+import os
+
+import pytest
+from sage.all import *  # noqa: F401
+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"))
+
+set_palp_dimension(11)
+
+TEST_MODELS_PARAMETRISED = [pytest.param(*row[1:], id=row[0]) for row in TEST_MODELS]
+
+
+@pytest.mark.parametrize(
+    "make_model, expected_cy_dimension, expected_no_divs, expected_disc, expected_intersection_numbers_CY",
+    TEST_MODELS_PARAMETRISED,
+)
+def test_topdata_and_disc(
+    make_model,
+    expected_cy_dimension,
+    expected_no_divs,
+    expected_disc,
+    expected_intersection_numbers_CY,
+):
+    polytope = make_model()
+
+    discriminants = polytope.disc()
+    assert [(str(factor), codim) for factor, codim in discriminants] == expected_disc
+
+    polytope.topdata()
+    assert polytope.cy_dimension == expected_cy_dimension
+    assert polytope.no_divs == expected_no_divs
+    assert polytope.intersection_numbers_CY == expected_intersection_numbers_CY