Implementing moving origin of ideal (#19)

---------

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

.gitea/ci-image/Dockerfile

@@ -0,0 +1,17 @@
+# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner,
+# since the official sagemath/sagemath image is amd64-only.
+#
+# 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
+
+FROM condaforge/miniforge3:latest
+
+RUN mamba install -y -c conda-forge sage ruff pytest \
+    && mamba clean -afy
+
+RUN mamba install -y -c conda-forge nodejs \
+    && mamba clean -afy
+
+WORKDIR /workspace

.gitea/preparse.sh

@@ -0,0 +1,26 @@
+#!/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.
+
+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"
+done

.gitea/workflows/ci.yml

@@ -0,0 +1,37 @@
+name: CI
+
+on:
+  push:
+    branches: ["**"]
+  pull_request:
+    branches: ["**"]
+
+jobs:
+  lint:
+    runs-on: ubuntu-latest
+    # Custom image (see .gitea/ci-image/Dockerfile)
+    container:
+      image: gitea.piribauer.ch/julian/sage-ci:latest
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Preparse .sage files
+        run: .gitea/preparse.sh
+
+      - name: ruff check
+        run: python -m ruff check --no-respect-gitignore .
+
+  test:
+    runs-on: ubuntu-latest
+    needs: lint
+    container:
+      image: gitea.piribauer.ch/julian/sage-ci:latest
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Run tests
+        run: python -m pytest tests/

.gitignore

@@ -0,0 +1,10 @@
+# Sage preparser output (regenerated from .sage sources, not hand-maintained)
+*.sage.py
+
+# Python / tooling caches
+__pycache__/
+.mypy_cache/
+.ruff_cache/
+.pytest_cache/
+
+.venv/

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

@@ -0,0 +1,20 @@
+[tool.ruff]
+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" = ["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"]

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

@@ -2,23 +2,23 @@ import numpy as np
 import logging
 import os
 import json
-from datetime import datetime, timezone
+import re
 
+# Logger
 logger = logging.getLogger(__name__)
-logging.basicConfig(level=logging.DEBUG)  # or logging.DEBUG to see both
+logging.basicConfig(level=logging.DEBUG)
 
 # For grep-commands
-import re 
+find_zs = re.compile(r'z\d+')
-find_zs = re.compile('z\d+')
+find_ls = re.compile(r'l\d+')
-find_ls = re.compile('l\d+')
+find_index_z = re.compile(r'\d+')
-find_index_z = re.compile('\d+')
 
 
 class Polytope:
     """
     Class representing a lattice polytope. Objects contain a list of defining points
     and a list of points in the convex hull of the polytope. Points inside codimension
-    are excluded from the convex hull and stored in a separate list.
+    one are excluded from the convex hull and stored in a separate list.
     """
 
     #############################
@@ -70,6 +70,8 @@ class ToricPolytope(Polytope):
     toric polytope describing the original family's ambient space.
 
     Example: For the mirror quintic, the vertices are e_1, ..., e_4, -e_1-...-e_4.
+
+    For CICYs, a nef-partition must be provided, defaulting to hypersurface if none given.
     """
 
     ##############################
@@ -99,6 +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)
         self.triangulation = self.triangulations[no_triangulation]
     
@@ -111,7 +114,7 @@ class ToricPolytope(Polytope):
                 raise ValueError("Invalid nef partition: {}".format(nef_partition))
             self.nef_partition = nef_partition
         else:
-            self.nef_partition = [range(len(self.relevant_lattice_points) - 1)]  # Default to a trivial partition if none is provided
+            self.nef_partition = [list(range(len(self.relevant_lattice_points) - 1))]  # Default to a trivial partition if none is provided
 
     def set_Mori_cone(self, lvec = None):
         """
@@ -127,14 +130,13 @@ class ToricPolytope(Polytope):
                 raise ValueError("All rows in lvec must have the same length.")
             elif len(lvec[0]) != len(self.relevant_lattice_points):
                 raise ValueError("The length of the l-vectors is invalid.")
-            elif matrix(matrix(default_Mori_cone.rays()).stack(matrix(lvec))).rank() != matrix(default_Mori_cone.rays()).rank():
+            elif matrix(default_Mori_cone.rays()).stack(matrix(lvec)).rank() != matrix(default_Mori_cone.rays()).rank():
                 raise ValueError("The provided l-vectors are not in the linear span of the original Mori cone.")
             else:
                 self.Mori_cone = Cone(lvec)
         
         if len(matrix(self.Mori_cone.rays()).kernel().gens()) > 0:
             logger.warning("Non-simplicial Mori-cone encountered.")
-            logger.info("Using the first %d linearly independent vectors.", len(self.Mori_cone.rays()))
             self.Mori_cone = Cone(
                                 transpose(
                                     transpose(self.Mori_cone.rays())
@@ -144,6 +146,7 @@ class ToricPolytope(Polytope):
                                     ).kernel().gens()
                                 ).transpose()
                              )
+            logger.info("Using the first %d linearly independent vectors.", len(self.Mori_cone.rays()))
 
         self.lvec = matrix(self.Mori_cone.rays())
 
@@ -161,11 +164,22 @@ class ToricPolytope(Polytope):
 
         self.set_triangulation(no_triangulation)
         self.set_nef_partition(nef_partition)
+        self.cy_dimension = self.dimension - len(self.nef_partition)
 
         self.fan = self.triangulation.fan(self.origin)
         self.ToricVariety = ToricVariety(self.fan)
         self.set_Mori_cone(lvec = lvec)
 
+        self.model_name = model_name if model_name else "Unnamed Model"
+        
+        logger.info(
+            "Initialised %s: %s%d with h21 = %d.",
+            self.model_name,
+            "CY" if len(self.nef_partition) == 1 else "CICY",
+            self.cy_dimension,
+            len(list(self.lvec)),
+        )
+
     ###############################
     # Discriminant Computation
     ###############################
@@ -219,7 +233,7 @@ class ToricPolytope(Polytope):
     def _setup_discriminant_symbols(self):
         mori_rays = self.Mori_cone.rays()
 
-        a_vars = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (1..len(mori_rays))]
+        a_vars = [var("a_{}".format(u), latex_name="a_{{{}}}".format(u)) for u in (1..len(mori_rays))]
         z_vars = var('z', n=len(mori_rays)+1, latex_name='z') # z[0] is superfluous
         lambda_vars = var('l', n=len(mori_rays)+1, latex_name='l') # l[0] is superfluous
         a_row = matrix(a_vars)
@@ -261,9 +275,10 @@ class ToricPolytope(Polytope):
             if polynomials[0] == 0:
                 polynomials = maxima.eliminate(equation_system, lambda_symbols[:-1]).sage()
                 reverse_solve = True
-        except:
+        except Exception:
             # In one-parameter cases there may be nothing to eliminate.
             polynomials = equation_system
+            logger.debug("No elimination needed for the equation system: %s", equation_system)
 
         try:
             if len(lambda_symbols) == 0:
@@ -274,7 +289,7 @@ class ToricPolytope(Polytope):
             else:
                 last_lambda = lambda_symbols[-1]
                 polynomials = [poly / last_lambda ** (poly.degree(last_lambda)) for poly in polynomials]
-        except:
+        except Exception:
             pass
 
         return polynomials
@@ -295,20 +310,22 @@ class ToricPolytope(Polytope):
                 codim = total_relation_blocks - 1 - index
             discriminants.append([polynomial, codim])
 
-    def disc(self, only_strong_coupling=False, no_triangulation=0):
+    def disc(self, only_strong_coupling=False):
         """
         Computes the A-discriminant for the toric polytope (following Aspinwall, Plesser, Wang),
-        which gives the singular loci of the moduli space in Batyrev coordinates.
+        which gives the singular loci of the moduli space in Batyrev coordinates. There is an option
-
+        to only compute strong coupling discriminant factors, so those coming from dependencies
-        Parameters:
+        in one-dimensional faces.
-        - only_strong_coupling: If True, only gives loci arising from edges (one-dimensional faces).
-        - no_triangulation: Index of the triangulation to use.
 
         Returns:
         A list of tuples [disc_i, codim_i] of discriminant factors disc_i
         coming from a relation inside a face of codimension codim_i.
         """
 
+        logger.info("Computing discriminant for %s", self.model_name)
+        if only_strong_coupling:
+            logger.info("Restricting to strong-coupling loci (codimension-1 faces).")
+
         all_linear_dependencies_among_points = self._collect_linear_dependencies(only_strong_coupling)
         self._append_origin_weight(all_linear_dependencies_among_points)
         a_vars, z_vars, lambda_vars, a_row, mori_matrix = self._setup_discriminant_symbols()
@@ -344,6 +361,9 @@ class ToricPolytope(Polytope):
             # Insert an empty codimension-0 entry for compatibility with prior behavior.
             discriminants = [[]] + discriminants
 
+        logger.info("Found %d discriminant factors.", len(discriminants))
+        logger.info("Discriminant factors: %s", discriminants)
+
         return discriminants
 
     ##############################
@@ -482,7 +502,7 @@ class ToricPolytope(Polytope):
                 messages.append("Possible elliptic fibration in the divisor class dual to t{}.".format(i))
 
         if messages:
-            logger.info("\n\n--- Elliptic fibrations --------------------------\n" + "\n".join(messages))
+            logger.info("\n\n--- Elliptic fibrations --------------------------\n%s", "\n".join(messages))
 
         return messages
 
@@ -609,11 +629,7 @@ class ToricPolytope(Polytope):
     def _write_output_json(self, output):
         """
         Writes the given output dictionary (from _get_output, i.e. exactly what is
-        logged) to data/<model_name>.json (relative to the repository root), if a
+        logged) to data/<model_name>.json."""
-        model_name is given. Plain JSON, so it is directly loadable in Python
-        (json.load), C++ (e.g. nlohmann::json), or any other language without any
-        custom parsing code; lvec loads back as a plain list of lists, e.g. lvec[0].
-        """
         if not self.model_name:
             logger.info("No model_name given, skipping writing topological data to JSON.")
             return None
@@ -631,24 +647,9 @@ class ToricPolytope(Polytope):
 
     def topdata(self):
         """
-        Computes topological data for hypersurfaces and CICYs in toric ambient spaces.
+        Computes topological data for hypersurfaces and CICYs in toric ambient spaces
-        For hypersurfaces, nef_partition should remain untouched (=0).
+        and saves the result to /data/topdata/<model_name>.json if a model_name is provided.
-        For CICYs with d polynomials, a nef_partition has to be supplied:
-        its format should be a list of d lists as in the examples below
-        which is a decomposition of the N points of $(polytope); the number
-        should corresponds to an enumeration of the points of $(polytope) with
-        the inner point omitted.
-
-        Conditions for possible fibrations are:
-            - elliptic:           if t_i^n = 0 and t_i^{n-1} != 0
-            - K3 (only 3-folds):  if c2.t_i = 24
-
-        Input:  - (for CICYs:) nef-partition
-                - (optional:) set of l-vectors to be used
-                - (optional:) index of triangulation
         """
-        no_polys = len(self.nef_partition)
-        self.cy_dimension = self.dimension - no_polys
         # Each point gives a toric divisor, each "column" gives one linear relation
         self.no_divs = len(self.relevant_lattice_points) - self.dimension - 1
 
@@ -731,7 +732,6 @@ class ToricPolytopeCICY(ToricPolytope):
 
     def __init__(self, CICY, no_triangulation=0, model_name=None, lvec=None):
         points, nef_partition = self._CICY_to_points(CICY)
-        print(points, nef_partition)
         super().__init__(
             points,
             no_triangulation=no_triangulation,
@@ -740,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/toric_topdata.sage.py

@@ -1,348 +0,0 @@
-
-
-# This file was *autogenerated* from the file sage/toric_topdata.sage
-from sage.all_cmdline import *   # import sage library
-
-_sage_const_0 = Integer(0); _sage_const_1 = Integer(1); _sage_const_2 = Integer(2); _sage_const_8 = Integer(8); _sage_const_10 = Integer(10); _sage_const_4 = Integer(4); _sage_const_24 = Integer(24)### FUNCTIONS: ####
-# - pointstopoly: converts list of points to CY-polytope
-# - disc: computes discriminant for polytope
-# - topdata: computes all sorts of topological data for polytope
-# - CICYtopdata: computes topdata in the simple CICY format, e.g. [[3,3]] for two cubics in P5 or [[3,0,1],[0,3,1]] for the Tian--Yau manifold
-
-
-import numpy as np
-# For grep-commands
-import re
-findzs = re.compile('z\d+')
-findls = re.compile('l\d+')
-findindexz = re.compile('\d+')
-
-
-def pointstopoly(points):
-    # Takes a list of points and returns the points inside its convex hull while omitting points inside faces of co-dimension one
-    
-    global pcodim1 #points in co-dimension one
-    
-    zeros = [_sage_const_0  for i in range(len(points[_sage_const_0 ]))] #origin
-    pc = LatticePolytope(points) # all points
-    polytope = [list(m) for m in pc.points()]
-    pcodim1 = []
-    for f in pc.facets():
-        for i in f.interior_point_indices():
-            try:
-                polytope.remove([list(m) for m in f.points(i)][_sage_const_0 ]) # remove those inside codim1 faces
-                pcodim1.append([list(m) for m in f.points(i)][_sage_const_0 ]) # save omitted points in pcodim1
-            except:
-                pass
-    # Moving zeros to the end
-    polytope.remove(zeros)
-    polytope.append(zeros)
-    return polytope
-
-def disc(polytope, only_sc=False, no_triangulation=_sage_const_0 ):
-    # Takes a poltyope and returns a list of tuples [disc_i,codim_i] of discriminant factors disc_i coming from a relation inside a face of codimension codim_i
-    
-    global globalks, pc
-    pc = LatticePolytope(polytope)
-    dimp = len(polytope[_sage_const_0 ])
-    zeros = [_sage_const_0  for i in range(dimp)]
-    # Compute all relations in all faces
-    globalks = []
-    if only_sc:
-        pcfaces = pc.faces()[:-_sage_const_1 ]
-    else:
-        pcfaces = pc.faces()
-     
-    for facesofdim in pcfaces: # for each face-dimension
-        atdim = [] 
-        for face in facesofdim: # for each face of fixed dimension
-            atface = []
-            facepoints = [list(i) for i in matrix(face.points())]
-            try:
-                facepoints.remove(zeros)
-            except:
-                pass
-
-            for p in facepoints.copy():
-                try:
-                    if p not in polytope:
-                        facepoints.remove(p)       # remove points that were omitted for polytope
-                except:
-                    pass
-
-            rels = {str(i):polytope.index(i) for i in facepoints}       # dictionary for getting positions in relations correct
-            kernel = matrix(facepoints).kernel().gens()
-            for k in kernel:                                     # just formatting
-                globalk = [_sage_const_0  for i in range(len(polytope)-_sage_const_1 )]    # MAY HAVE TO BE ADAPTED FOR CASES WITH MORE VERTICES
-                for i in range(len(k)):
-                    globalk[rels[str(facepoints[i])]] = k[i]
-                atdim.append(globalk)
-        globalks.append([[i[j] for j in range(len(i))] for i in np.unique(np.matrix(atdim),axis=_sage_const_0 )])
-    for atdim in globalks: # inserting weight for innter point
-        for k in atdim:
-            if len(k)>_sage_const_0 :
-                k.append(-sum(k))
-    # Computing discriminants 
-    pc = PointConfiguration(polytope)
-    pc_star =  pc.restrict_to_star_triangulations(zeros)
-    pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
-    if(len(pc_star_and_fine.triangulations_list())==_sage_const_0 ):
-        print("No fine star triangulation!")
-        return _sage_const_0 
-    if(len(pc_star_and_fine.triangulations_list())>_sage_const_1 ):
-        print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
-    triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
-    fan=triangulation.fan(zeros)
-    X=ToricVariety(fan)
-    ls = [l for l in X.Mori_cone().rays()]
-    if (len(matrix(ls).kernel().gens())>_sage_const_0 ):
-        print("Non-simplicical Mori-cone")
-
-    As = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (ellipsis_iter(_sage_const_1 ,Ellipsis,len(ls)))]
-    z = var('z', n=len(ls)+_sage_const_1 , latex_name='z') # z[0] is superfluous
-    l = var('l', n=len(ls)+_sage_const_1 , latex_name='l') # l[0] is superfluous
-    #globalks.reverse()
-    globalks = [k for k in globalks if k not in ([[]],)]
-    discs = []
-    for i,atdim in enumerate(globalks):
-        solsys = []
-        for kindex in range(len(atdim)):
-            sol=solve((matrix(As)*matrix(ls)-matrix(atdim[kindex])).list(),matrix(As).list())
-            zindex = [abs(i.subs(sol)) for i in As].index(_sage_const_1 ) +_sage_const_1 
-            solsys.append(z[zindex]-prod([sum([atdim[j][i]*l[j] for j in range(len(atdim))])**(atdim[kindex][i]) for i in range(len(atdim[_sage_const_0 ]))]))
-            lambdas = list(set(findls.findall(str(solsys))))
-            lambdas.sort()
-        try:
-            pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[_sage_const_1 :]]).sage()
-            if pol[_sage_const_0 ]==_sage_const_0 :
-                pol = maxima.eliminate(solsys,[eval(i) for i in lambdas[:-_sage_const_1 ]]).sage()
-                reversesolve = True
-            else:
-                reversesolve = False
-        except:  # this is just for 1-parameter cases where there is nothing to solve
-            pol = solsys
-            pass
-        try:
-            if not reversesolve:
-                pol = [i/l0**(i.degree(l0)) for i in pol]
-            else:
-                pol = [i/eval(lambdas[-_sage_const_1 ])**(i.degree(eval(lambdas[-_sage_const_1 ]))) for i in pol]
-        except:
-            pass
-        for poli in pol:
-            if (poli not in [a[_sage_const_0 ] for a in discs] and poli!=_sage_const_0 ):
-                if only_sc:
-                    discs.append([poli,len(globalks)-i]) # account for offset
-                else:
-                    discs.append([poli,len(globalks)-_sage_const_1 -i])
-                
-    if only_sc: #insert an empty list for the codimension 0 face discriminant
-        discstmp = discs
-        discstmp.reverse()
-        discstmp.append([])
-        discstmp.reverse()
-        discs = discstmp
-    return discs
-
-def topdata(polytope, nef_partition=_sage_const_0 , lvec=_sage_const_0 , no_triangulation=_sage_const_0 , returnvars=False):
-    # Computes topological data for hypersurfaces and CICYs in toric ambient spaces.
-    # For hypersurfaces, nef_partition should remain untouched (=0).
-    # For CICYs with d polynomials, a nef_partition has to be supplied:
-    # its format should be a list of d lists as in the examples below
-    # which is a decomposition of the N points of $(polytope); the number
-    # should corresponds to an enumeration of the points of $(polytope) with
-    # the inner point omitted.
-
-    # Conditions for possible fibrations are:
-    #         - elliptic:           if t_i^n = 0 and t_i^{n-1} != 0
-    #         - K3 (only 3-folds):  if c2.t_i = 24
-
-    # input:  - array of points in polytope
-    #         - (for CICYs:) nef-partition
-    #         - (optional:) set of l-vectors to be used
-    #         - (optional:) index of triangulation
-
-    # output: - Mori-cone generators
-    #         - Intersection rings of CY and ambient space
-    #         - Topological data
-    #         - Possible fibrations
-    #         - GLOBALS: triangulation and Mori-cone generators
-
-    global triangulation, lall, MoriMatrix, fan, X, Ydualform, J, D
-
-    pc = PointConfiguration(polytope)
-
-    dimp=len(polytope[_sage_const_0 ]) # dimension of polytope
-    if( nef_partition == _sage_const_0 ):
-        no_polys = _sage_const_1 
-        nef_partition = [[i for i in range(len(polytope)-_sage_const_1 )]] # trivial partition
-    else:
-        if( sorted(flatten(nef_partition)) != [ i for i in range(len(polytope)-_sage_const_1 )] ):
-            print("Nef-partition not valid!")
-        no_polys = len(nef_partition)
-    zeros=zero_vector(dimp)
-
-    pc_star =  pc.restrict_to_star_triangulations(zeros)
-    pc_star_and_fine=pc_star.restrict_to_fine_triangulations()
-    if(len(pc_star_and_fine.triangulations_list())==_sage_const_0 ):
-        print("No fine star triangulation!")
-        return _sage_const_0 
-    if(len(pc_star_and_fine.triangulations_list())>_sage_const_1 ):
-        print("More than one fine star triangulation! (",len(pc_star_and_fine.triangulations_list()),")")
-    triangulation=pc_star_and_fine.triangulations_list()[no_triangulation]
-    fan=triangulation.fan(zeros)
-    X=ToricVariety(fan)
-
-    # Change l-vectors if given:
-    if lvec==_sage_const_0 :
-        lall=matrix(X.Mori_cone().rays())
-    else:
-        lall=lvec
-    nodivs = len(polytope)-dimp-_sage_const_1 
-    MoriMatrix=matrix([[-sum([lall[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), lall[:,:-_sage_const_1 ]
-    # Check whether Mori-cone is simplicial and continue with first $(nodivs) vectors
-    if ( lall.dimensions()[_sage_const_0 ] > nodivs ):
-        print("Non-simplicial Kähler-cone! (",lall.dimensions()[_sage_const_0 ]," > ",nodivs,")")
-        print("Picking ",nodivs," linearly independent vectors.")
-        l = transpose(transpose(lall)*(matrix(transpose(matrix(lall.kernel().gens())).kernel().gens()).transpose()))
-        MoriMatrix=matrix([[-sum([l[j][i] for i in part]) for part in nef_partition] for j in range(nodivs)]), l[:,:-_sage_const_1 ]
-    else:
-        l = lall
-
-    HH=X.cohomology_ring()
-    D = [HH(X.divisor(i)) for i in (ellipsis_range(_sage_const_0 ,Ellipsis,len(polytope)-_sage_const_2 ))]
-    zs = list(set(findzs.findall(str(D)))) # gives the $(nodivs) z-variables present in divisor classes
-    # Find Kähler-cone generators as duals to l-vectors:
-    Bsinv=l.matrix_from_columns([eval(m) for m in sorted(list(findindexz.findall(str(zs))))])
-    if (Bsinv.det()==_sage_const_0 ):
-        print("l-vectors not independent in divisor basis. Pick them manually with ``lvec=...''")
-        return
-    Bs = transpose(Bsinv**(-_sage_const_1 ))
-    J = Bs*vector([D[i] for i in list(set([eval(i) for i in findindexz.findall(str(findzs.findall(str(D))))]))])
-    # Hyperplane divisor class:
-    Ydualform = product([sum([HH(D[p]) for p in part]) for part in nef_partition])
-
-    # Tuple lists for computations below
-    TupleListCY = UnorderedTuples(range(nodivs),dimp-no_polys)
-    TupleListCYordered = Tuples(range(nodivs),dimp-no_polys) # for intersection ring multiplicities must be included
-    TupleListAm = UnorderedTuples(range(nodivs),dimp)
-
-    # Find intersection numbers on CY and on ambient space:
-    intCY = [list((tl,X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])))) for tl in TupleListCY]
-    intAm = [list((tl,X.integrate(product([J[tl[j]] for j in range(dimp)])))) for tl in TupleListAm]
-
-    # Form intersection ring on CY:
-    t = var('t', n=nodivs, latex_name='t')
-    intring = sum([product([t[index] for index in tl])*X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])) for tl in TupleListCYordered])
-    intringnomults = sum([product([t[index] for index in tl])*X.integrate(Ydualform*product([J[tl[j]] for j in range(dimp-no_polys)])) for tl in set([tuple(sorted(i)) for i in TupleListCYordered])])
-    # Compute Chern-character for topological data:
-    var(findzs.findall(str([J[i] for i in range(nodivs)])))                    # Introduces all z in Kähler-forms as variables...
-    zs = [eval(z) for z in findzs.findall(str([J[i] for i in range(nodivs)]))] # ... and puts them into a vector.
-    tsubs=solve([lift(J[i])==t[i] for i in range(nodivs)],zs)       # Finds substitution rule for zs in terms Kähler-cone generators
-    cc = (product([_sage_const_1 +eval(str((lift(d)))).subs(tsubs[_sage_const_0 ]) for d in D]))/(product([_sage_const_1 +sum([eval(str(lift(HH(D[p])))) for p in part]) for part in nef_partition]).subs(tsubs[_sage_const_0 ]))    # Adjunction-formula
-    
-    print('--- Toric divisors (ambient space) -----------')
-    print(D)
-    print('\n--- Kähler cone generators (ambient space) --- ')
-    print(J)
-    print('\n--- Mori-cone-generators (ambient space) ------')
-    for j in range(nodivs):
-        print(MoriMatrix[_sage_const_0 ][j],MoriMatrix[_sage_const_1 ][j])
-
-    print('\n--- Intersection on CY ------------------------')
-    Ccy=[product([t[iden[_sage_const_0 ][i]] for i in range(dimp-no_polys)])==iden[_sage_const_1 ] for iden in intCY]
-    print(Ccy)
-    print('\nR = ',intring)
-    print('\nR (no multiplicities) = ',intringnomults)
-
-    print('\n--- Intersection in ambient space -------------')
-    Cam=[product([t[iden[_sage_const_0 ][i]] for i in range(dimp)])==iden[_sage_const_1 ] for iden in intAm]
-    print(Cam)
-
-    print('\n--- Topological data --------------------------')
-    #print((product([1+lift(d) for d in D]))/(product([1+sum([lift(HH(D[p])) for p in part]) for part in nef_partition])))
-    chi = X.integrate(Ydualform*(product([_sage_const_1 +lift(d) for d in D]))/(product([_sage_const_1 +sum([lift(HH(D[p])) for p in part]) for part in nef_partition])))
-    print('chi = ',chi)
-
-    c = var('c')
-    if( round((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,_sage_const_1 ).coefficient(c,_sage_const_1 ).subs({t[i]:_sage_const_1  for i in range(nodivs)}),_sage_const_8 ) != _sage_const_0 ):
-        return "First Chern-class not zero!"
-
-    # Print Chern classes:
-    cJ = [[_sage_const_0 ]]*(dimp-no_polys)
-    for i in range(_sage_const_2 ,dimp-no_polys):
-        print('\nc',i,' = ',(cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,i).coefficient(c,i))
-        TupleListi = UnorderedTuples(range(nodivs),i)
-        TupleListz = UnorderedTuples(range(nodivs),dimp-no_polys-i)
-        cJ[i] = [sum([((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,i).coefficient(c,i)).coefficient(product([t[i] for i in tupl]))*((product([t[j] for j in tuplz])*product([t[i] for i in tupl])).subs(Ccy)) for tupl in TupleListi]) for tuplz in TupleListz]
-        print('integrated: '+str([cJ[i][j]*product([t[k] for k in TupleListz[j]]) for j in range(len(TupleListz))]))
-    print('\nc',dimp-no_polys,' = ',(cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,dimp-no_polys).coefficient(c,dimp-no_polys))
-    
-    msg = ""
-    # elliptic:
-    ## Indicates whether J_i^n==0 with J_i^(n-1)!=0 for n-folds and some J_i
-    for i in range(nodivs):
-        for intn in intCY:
-            if( intn[_sage_const_0 ] == [i]*(dimp-no_polys) and intn[_sage_const_1 ] == _sage_const_0 ):
-                if not(all([round((t[i]**(dimp-no_polys-_sage_const_1 )*t[j]).subs(Ccy),_sage_const_10 ) == _sage_const_0  for j in range(nodivs) if j != i]) ):
-                    msg += "Possible elliptic fibration in cycle dual to t"+str(i)+".\n"
-
-    # K3 for 3-folds
-    ## Indicates whether c2.J_i =24 for some J_i
-    if( dimp+_sage_const_1 -no_polys == _sage_const_4 ):
-        TupleList2 = UnorderedTuples(range(nodivs),_sage_const_2 )
-        for j in range(nodivs):
-            if( sum([((cc).subs({t:c*t for t in t}).taylor(c,_sage_const_0 ,_sage_const_2 ).coefficient(c,_sage_const_2 )).coefficient(product([t[i] for i in tupl]))*((product([t[i] for i in tupl])*t[j]).subs(Ccy)) for tupl in TupleList2]) == _sage_const_24  ):
-                msg += "Possible K3-fibration in divisor t"+str(j)+".\n"
-
-    if( msg != "" ):
-        print('\n--- Fibrations --------------------------------')
-        print(msg)
-    if returnvars:
-        return str([(MoriMatrix[_sage_const_0 ][j],MoriMatrix[_sage_const_1 ][j]) for j in range(nodivs)]).replace("(","{").replace(")","}").replace("[","{").replace("]","}") # might return more in the future if necessary
-    
-    
-def CICYtopdata(CICY,justdata=False,justmori=False, lvec=_sage_const_0 ):
-    # Input: a list l with entries l[i,j] that give the weight in the 
-    #        ambient projective space i of polynomial j.
-    #        E.g. (P^3| 3 1)
-    #             (P^2| 2 0)
-    #        corresponds to the list ((3,1),(2,0))
-    
-    ambients=(np.array(CICY).transpose()).sum(axis=_sage_const_0 )-_sage_const_1 
-    if ambients in ZZ:
-        ambients=[ambients]
-    dimp=sum(ambients)
-    polytope=identity_matrix(int(dimp))
-    zeros=[_sage_const_0 ]*dimp
-    polytope=polytope.insert_row(_sage_const_0 ,zeros)
-    
-    offset=_sage_const_0 
-    for i in range(len(ambients)):
-        newrow=[_sage_const_0 ]*dimp
-        entry=[]
-        for j in range(ambients[i]):
-            newrow[offset]=-_sage_const_1 
-            offset+=_sage_const_1 
-            entry.append(offset)
-        polytope=polytope.insert_row(sum([ambients[n]+_sage_const_1  for n in (ellipsis_range(_sage_const_0 ,Ellipsis,i))]),newrow)
-    
-    partition=[]
-    eqs=matrix(CICY).transpose()
-    counters=[_sage_const_0 ]*len(ambients)
-    for eq in eqs:
-        part=[]
-        for j in range(len(eq)):
-            part.append([sum([ambients[m]+_sage_const_1  for m in (ellipsis_range(_sage_const_0 ,Ellipsis,(j-_sage_const_1 )))])+counters[j]+a for a in (ellipsis_range(_sage_const_0 ,Ellipsis,(eq[j]-_sage_const_1 )))])
-            counters[j]+=len(part[j])
-        partition.append(flatten(part))
-    
-    if justdata:
-        return [list(polytope),partition]
-    elif justmori:
-        moricone((list(polytope)))
-    else:
-        topdata(list(polytope),nef_partition=partition,lvec=lvec)
-
-topdata(pointstopoly([(_sage_const_1 ,_sage_const_0 ),(_sage_const_0 ,_sage_const_1 ),(-_sage_const_1 ,-_sage_const_1 )]))
-

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_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