Merge branch '5-pipeline-with-unit-tests' of https://gitea.piribauer.ch/julian/Calabi-Yau-Period-Geometry into 5-pipeline-with-unit-tests

.gitea/ci-image/Dockerfile

@@ -1,16 +1,26 @@
-# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner,
-# since the official sagemath/sagemath image is amd64-only.
+# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner.
 #
-# Build & push commands:
+# 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):
 #   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,26 +1,15 @@
 #!/usr/bin/env bash
-
-# 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.
+# 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
+# `sage --python -m ruff check --no-respect-gitignore .` to catch issues that
+# only show up in the generated output (e.g. missing trailing newlines).
+#
+# The generated files are gitignored; clean them up afterwards with:
+#   git clean -x sage/ playground/
 
 set -e
 
-if command -v python >/dev/null 2>&1 && python -c "import sage.repl" >/dev/null 2>&1; then
-  run_python() { python "$@"; }
-else
-  run_python() { sage --python "$@"; }
-fi
-
 for f in $(find . -name "*.sage" -not -path "./playground/*"); do
-  run_python -c "
-import sys
-from sage.repl.preparse import preparse_file
-
-path = sys.argv[1]
-with open(path) as fh:
-    src = fh.read()
-with open(path + '.py', 'w') as fh:
-    fh.write('from sage.all import *  # noqa: F401,F403\n')
-    fh.write(preparse_file(src))
-" "$f"
+  sage --preparse "$f"
+  printf 'from sage.all import *  # noqa: F401,F403\n%s\n' "$(cat "${f}.py")" > "${f}.py"
 done

.gitea/workflows/ci.yml

@@ -9,7 +9,9 @@ on:
 jobs:
   lint:
     runs-on: ubuntu-latest
-    # Custom image (see .gitea/ci-image/Dockerfile)
+    # 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.
     container:
       image: gitea.piribauer.ch/julian/sage-ci:latest
 
@@ -18,9 +20,18 @@ 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,151 +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([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,18 +3,33 @@ 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 preparser output: star-import globals, semicolon preamble, and literal-substitution artifacts trip static analysis.
+# *.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" = ["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 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"]

sage/period_computation.sage

@@ -1,646 +0,0 @@
-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,8 +101,7 @@ 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)
+        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):
@@ -740,34 +739,3 @@ 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

@@ -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))

tests/test_models.sage

@@ -1,104 +0,0 @@
-# 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

@@ -1,203 +0,0 @@
-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

@@ -0,0 +1,32 @@
+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

@@ -1,34 +0,0 @@
-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