import numpy as np
import logging
import os
import json
import re

# Logger
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.DEBUG)

# For grep-commands
find_zs = re.compile(r'z\d+')
find_ls = re.compile(r'l\d+')
find_index_z = re.compile(r'\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
    one are excluded from the convex hull and stored in a separate list.
    """

    #############################
    # Initialisation
    #############################

    def _validate_points_and_get_dimension(self) -> int:
        if len(self.defining_points) < 1:
            raise ValueError("At least one point is required.")
        elif len(set([len(p) for p in self.defining_points])) != 1:
            raise ValueError("All points must have the same dimension.")
        elif len(self.defining_points) <= len(self.defining_points[0]):
            raise ValueError("The number of points must be greater than the dimension of the points.")
        return len(self.defining_points[0])

    def _points_to_polytope(self) -> tuple:
        """
        Takes the defining points of the polytope and returns the points inside its convex hull
        while removing and returning points inside faces of co-dimension one
        as a second return value.
        """
        lattice_polytope = LatticePolytope(self.defining_points)
        points_in_convex_hull = [list(m) for m in lattice_polytope.points()]
    
        points_in_codim1_faces = [] # is populated below
        for facet in lattice_polytope.facets():
            for i in facet.interior_point_indices():
                points_inside_facet = [list(m) for m in facet.points(i)][0]
                if points_inside_facet in points_in_convex_hull:
                    points_in_convex_hull.remove(points_inside_facet) # remove those inside codim1 faces
                    points_in_codim1_faces.append(points_inside_facet) # save omitted points in pcodim1
    
        # Moving zeros to the end
        points_in_convex_hull.remove(self.origin)
        points_in_convex_hull.append(self.origin)

        return lattice_polytope, points_in_convex_hull, points_in_codim1_faces

    def __init__(self, points):
        self.defining_points = points
        self.dimension = self._validate_points_and_get_dimension()
        self.origin = [0] * self.dimension
        self.lattice_polytope, self.relevant_lattice_points, self.points_in_codim1_faces = self._points_to_polytope()

class ToricPolytope(Polytope):
    """
    Used to represent the anti-canonical Weyl divisor -K of a mirror family.
    By definition, it is also given by the polar dual polytope to the 
    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.
    """

    ##############################
    # Initialisation
    ##############################

    def _set_fine_star_triangulations(self):
        """
        Returns a list of fine star triangulations of the polytope.
        """
        self.point_configuration = PointConfiguration(self.relevant_lattice_points)
        star_triangulations = self.point_configuration.restrict_to_star_triangulations(self.origin)
        fine_star_triangulations = star_triangulations.restrict_to_fine_triangulations()

        if len(fine_star_triangulations.triangulations_list()) == 0:
            raise ValueError("No fine star triangulation found for the given polytope.")
        if len(fine_star_triangulations.triangulations_list()) > 1:
            logger.warning("More than one fine star triangulation! (%d)", len(fine_star_triangulations.triangulations_list()))
        
        self.triangulations = fine_star_triangulations.triangulations_list()
        
    def set_triangulation(self, no_triangulation):
        """
        Sets the triangulation of the polytope based on the provided index,
        ensuring that the index is valid.
        """
        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]
    
    def set_nef_partition(self, nef_partition):
        """
        Sets the nef partition for the polytope, ensuring that it is valid.
        """
        if nef_partition is not None:
            if sorted(flatten(nef_partition)) != list(range(len(self.relevant_lattice_points) - 1)):
                raise ValueError("Invalid nef partition: {}".format(nef_partition))
            self.nef_partition = nef_partition
        else:
            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):
        """
        Sets the Mori cone of the toric variety associated with the polytope.
        """
        default_Mori_cone = self.ToricVariety.Mori_cone()
        if lvec is None:
            self.Mori_cone = default_Mori_cone
        else:
            if not isinstance(lvec, (list, tuple, Matrix)):
                raise TypeError("lvec must be a list, tuple, or Matrix.")
            elif len(set([len(row) for row in lvec])) != 1:
                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(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.")
            self.Mori_cone = Cone(
                                transpose(
                                    transpose(self.Mori_cone.rays())
                                    *
                                    matrix(
                                            transpose(matrix(self.Mori_cone.rays()).kernel().gens())
                                    ).kernel().gens()
                                ).transpose()
                             )
            logger.info("Using the first %d linearly independent vectors.", len(self.Mori_cone.rays()))

        self.lvec = matrix(self.Mori_cone.rays())

    def __init__(
            self,
            points,
            no_triangulation = 0,
            model_name = None,
            nef_partition = None,
            lvec = None, # Note that the inner point is the last entry
        ):
        super().__init__(points)
        self.model_name = model_name
        self._set_fine_star_triangulations()

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

    def _collect_linear_dependencies(self, only_strong_coupling):
        # Compute linear relations in all faces and map them to global coordinates.
        all_linear_dependencies_among_points = []

        polytope_faces = self.lattice_polytope.faces()
        relevant_faces = polytope_faces[:-1] if only_strong_coupling else polytope_faces

        for faces_of_given_dimension in relevant_faces:
            relations_at_given_dimension = []
            for face in faces_of_given_dimension:
                points_in_face = [list(point) for point in matrix(face.points())]
                if self.origin in points_in_face:
                    points_in_face.remove(self.origin)

                # Keep only points that survived codim-1 face filtering.
                points_in_face = [point for point in points_in_face if point in self.relevant_lattice_points]
                if len(points_in_face) == 0:
                    continue

                # In principle, the relations are just the kernel of the matrix of points in the face, but
                # to express the discriminants in Batyrev coordinates, we need to map them to the global coordinates
                # of the polytope.
                point_index_dic = {str(point): self.relevant_lattice_points.index(point) for point in points_in_face}
                linear_dependencies_in_face = matrix(points_in_face).kernel().gens()
                for generator in linear_dependencies_in_face:
                    # Exclude origin here; its weight is appended later as minus the sum.
                    generator_as_global_relation = [0] * (len(self.relevant_lattice_points) - 1)
                    for idx in range(len(generator)):
                        generator_as_global_relation[point_index_dic[str(points_in_face[idx])]] = generator[idx]
                    relations_at_given_dimension.append(generator_as_global_relation)

            if len(relations_at_given_dimension) == 0:
                unique_relations = []
            else:
                unique_relations = np.unique(np.matrix(relations_at_given_dimension), axis=0).tolist()
            all_linear_dependencies_among_points.append(unique_relations)

        return all_linear_dependencies_among_points

    def _append_origin_weight(self, all_linear_dependencies_among_points):
        # Adding the origin weight to each relation.
        for relations_at_dim in all_linear_dependencies_among_points:
            for relation in relations_at_dim:
                if len(relation) > 0:
                    relation.append(-sum(relation))

    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))]
        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)
        mori_matrix = matrix(mori_rays)

        return a_vars, z_vars, lambda_vars, a_row, mori_matrix

    def _build_equation_system(self, relations_at_given_dimension, a_vars, z_vars, lambda_vars, a_row, mori_matrix):
        equation_system = []
        number_of_points = len(relations_at_given_dimension[0])

        # Eq. (8)-style building blocks: linear forms in lambda_a for each point index i.
        linear_forms = [
            sum([
                relations_at_given_dimension[relation_id][point_id] * lambda_vars[relation_id]
                for relation_id in range(len(relations_at_given_dimension))
            ])
            for point_id in range(number_of_points)
        ]

        for relation in relations_at_given_dimension:
            # This determines the Batyrev coordinate the dependence represents:
            solution = solve((a_row * mori_matrix - matrix(relation)).list(), a_row.list())
            z_index = [abs(ai.subs(solution)) for ai in a_vars].index(1) + 1

            equation_system.append(
                z_vars[z_index] - prod([linear_forms[point_id] ** relation[point_id] for point_id in range(number_of_points)])
            )

        return equation_system

    def _eliminate_and_normalize(self, equation_system, lambda_vars):
        lambda_names = sorted(list(set(find_ls.findall(str(equation_system)))))
        lambda_symbols = [lambda_vars[int(name[1:])] for name in lambda_names]

        reverse_solve = False
        try:
            polynomials = maxima.eliminate(equation_system, lambda_symbols[1:]).sage()
            if polynomials[0] == 0:
                polynomials = maxima.eliminate(equation_system, lambda_symbols[:-1]).sage()
                reverse_solve = True
        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:
                raise ValueError("No lambda variables found in elimination system.")
            if not reverse_solve:
                lambda0 = lambda_symbols[0]
                polynomials = [poly / lambda0 ** (poly.degree(lambda0)) for poly in polynomials]
            else:
                last_lambda = lambda_symbols[-1]
                polynomials = [poly / last_lambda ** (poly.degree(last_lambda)) for poly in polynomials]
        except Exception:
            pass

        return polynomials

    def _append_unique_discriminants(self, discriminants, seen_discriminants, polynomials, index, total_relation_blocks, only_strong_coupling):
        for polynomial in polynomials:
            if polynomial == 0:
                continue

            polynomial_key = str(polynomial)
            if polynomial_key in seen_discriminants:
                continue

            seen_discriminants.add(polynomial_key)
            if only_strong_coupling:
                codim = total_relation_blocks - index
            else:
                codim = total_relation_blocks - 1 - index
            discriminants.append([polynomial, codim])

    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. There is an option
        to only compute strong coupling discriminant factors, so those coming from dependencies
        in one-dimensional faces.

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

        all_linear_dependencies_among_points = [relations for relations in all_linear_dependencies_among_points if relations]
        discriminants = []
        seen_discriminants = set()
        total_relation_blocks = len(all_linear_dependencies_among_points)

        for index, relations_at_given_dimension in enumerate(all_linear_dependencies_among_points):
            if len(relations_at_given_dimension) == 0:
                continue

            equation_system = self._build_equation_system(
                relations_at_given_dimension,
                a_vars,
                z_vars,
                lambda_vars,
                a_row,
                mori_matrix,
            )
            polynomials = self._eliminate_and_normalize(equation_system, lambda_vars)
            self._append_unique_discriminants(
                discriminants,
                seen_discriminants,
                polynomials,
                index,
                total_relation_blocks,
                only_strong_coupling,
            )

        if only_strong_coupling:
            # 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

    ##############################
    # Topological Data Computation
    ##############################

    def _compute_Kahler_generators(self):
        """
        Computes the Kähler cone generators J_a dual to the Mori cone generators l^a.

        The toric divisors D_i satisfy ∑_i l^a_i D_i = 0 in cohomology, so only a
        subset of them are linearly independent. We identify the free generators via
        the z-variables appearing in the cohomology ring, extract the corresponding
        columns of the Mori matrix l, invert it, and express J = Bs * D_basis.
        """
        z_names = list(set(find_zs.findall(str(self.D))))
        z_indices = sorted([int(find_index_z.search(name).group()) for name in z_names])
        Bsinv = self.lvec.matrix_from_columns(z_indices)
        if Bsinv.det() == 0:
            raise ValueError("l-vectors are not independent in the divisor basis. Supply them manually via lvec=...")
        Bs = transpose(Bsinv**(-1))
        return Bs * vector([self.D[i] for i in z_indices])

    def _compute_intersection_ring(self):
        """
        Builds the intersection ring polynomial by adding up the intersection numbers
        in self.intersection_numbers_CY.
        """
        # Full ring: ordered tuples, so t_i*t_j and t_j*t_i both contribute.
        intring = sum(
            product([self.t[i] for i in tl]) * self.intersection_numbers_CY[tuple(sorted(tl))]
            for tl in Tuples(range(self.no_divs), self.cy_dimension)
        )
        # Ring without multiplicities: each distinct monomial once.
        intringnomults = sum(
            product([self.t[i] for i in tl]) * val
            for tl, val in self.intersection_numbers_CY.items()
        )
        return intring, intringnomults

    def _intersection_substitution_rules(self):
        """
        Builds monomial substitution rules from intersection numbers, e.g.
        t0^3 -> kappa_(0,0,0), t0*t1^2 -> kappa_(0,1,1).
        """
        rules = {}
        for idx_tuple, value in self.intersection_numbers_CY.items():
            monomial = product([self.t[i] for i in idx_tuple]).expand()
            rules[monomial] = value
        return rules

    def _replace_intersection_monomials(self, expression):
        """
        Replaces top-degree Kähler monomials in an expression by CY intersection numbers.
        """
        expanded_expression = expression.expand()
        replaced_expression = expanded_expression

        # Replace each degree-n monomial by its corresponding intersection number,
        # preserving the polynomial coefficient in front of that monomial.
        for idx_tuple, value in self.intersection_numbers_CY.items():
            monomial = product([self.t[i] for i in idx_tuple]).expand()
            coeff = expanded_expression.coefficient(monomial)
            if coeff != 0:
                replaced_expression -= coeff * monomial
                replaced_expression += coeff * value

        return replaced_expression.expand()

    def _compute_Chern_character(self):
        """
        Computes the total Chern class of the CY via the adjunction formula:

            c(Y) = c(TX) / c(N_{Y/X})

        where c(TX) = ∏_i (1 + D_i) and c(N_{Y/X}) = ∏_k (1 + Y_k*) with Y_k*
        the nef partition divisor classes. Everything is expressed in the Kähler
        basis by solving lift(J_i) = t_i for the cohomology ring z-variables.
        """
        z_names = find_zs.findall(str([self.J[i] for i in range(self.no_divs)]))
        var(z_names)
        zs = [eval(z) for z in z_names]
        tsubs = solve([lift(self.J[i]) == self.t[i] for i in range(self.no_divs)], zs)

        numerator = product([1 + eval(str(lift(d))).subs(tsubs[0]) for d in self.D])
        denominator = product([
            1 + sum([eval(str(lift(self.HH(self.D[p])))) for p in part])
            for part in self.nef_partition
        ]).subs(tsubs[0])
        return numerator / denominator

    def _compute_Chern_polynomials(self, Chern_character):
        """
        Computes the Chern polynomials of the CY by expanding the total Chern character
        and substituting intersection numbers.
        Returns both a Chern polynomial for each degree restricted to suitable intersections
        of Kähler divisors and their integrated versions.
        """
        c = var('c')

        Chern_character_taylored = (Chern_character).subs({t: c * t for t in self.t}).taylor(c, 0, self.cy_dimension)

        Chern_polynomials = [
            [
                (prod(self.t[i] for i in tple) * Chern_character_taylored.coefficient(c, k)).expand()
                for tple in UnorderedTuples(range(self.no_divs), self.cy_dimension - k)
            ]
            for k in range(1, self.cy_dimension + 1) # for each degree
        ] 

        integrated_Chern_classes = [
            [
                self._replace_intersection_monomials(poly) for poly in polys
            ]
            for polys in Chern_polynomials
        ]

        return Chern_polynomials, integrated_Chern_classes

    def _check_for_fibrations(self):
        """
        Checks for possible elliptic fibrations via the condition
            elliptic in direction i <=> J_i^n == 0 but J_i^(n-1)*J_j != 0 for some j != i.
        """
        messages = []

        for i in range(self.no_divs):
            self_intersection = self.intersection_numbers_CY[tuple([i] * self.cy_dimension)]
            if self_intersection != 0:
                continue
            has_nonzero_neighbor = any(
                self.intersection_numbers_CY[tuple(sorted([i] * (self.cy_dimension - 1) + [j]))] != 0
                for j in range(self.no_divs) if j != i
            )
            if has_nonzero_neighbor:
                messages.append("Possible elliptic fibration in the divisor class dual to t{}.".format(i))

        if messages:
            logger.info("\n\n--- Elliptic fibrations --------------------------\n%s", "\n".join(messages))

        return messages

    def _get_polytope_and_glsm_table(self):
        """
        Builds a list of lines showing the polytope points as columns of a
        (dimension x n) matrix (one column per point), with the GLSM charge
        l-vectors listed directly underneath: one row per Mori cone
        generator, with each l-vector entry aligned below the point
        coordinate it belongs to. Returned as a list of lines (rather than a
        single newline-joined string) so that each line becomes its own JSON
        array entry and stays on its own line when written to file.
        
        To-do: add inner point and CICY functionality.
        """
        points = list(self.relevant_lattice_points) # create a copy
        points.remove(self.origin)  # Remove the origin point
        n_outer_points = len(points)

        nef_unit_vectors = [tuple([1 if i == j else 0 for j in range(len(self.nef_partition))]) for i in range(len(self.nef_partition))]
        for i, partition in enumerate(self.nef_partition):
            for p in partition:
                points[p] = tuple(list(nef_unit_vectors[i]) + list(points[p]))  # Append nef partition index to each point
            points.append(tuple(list(nef_unit_vectors[i]) + [0] * self.dimension))  # Add a point for the nef partition itself

        point_rows = [[str(point[row]) for point in points] for row in range(self.dimension + len(self.nef_partition))]

        lvec_rows = self._get_lvec_rows()
        lvec_rows_string = [[str(entry) for entry in row] for row in lvec_rows]  # Exclude the last entry of each l-vector (origin point)

        column_widths = [
            max(len(row[col]) for row in point_rows + lvec_rows_string)
            for col in range(n_outer_points + len(self.nef_partition))
        ]

        def format_row(row_entries, separator):
            # Double the separator right before the trailing nef-partition-marker columns.
            row_string = row_entries[0].rjust(column_widths[0])
            for col in range(1, len(row_entries)):
                sep = separator * 2 if col == n_outer_points else separator
                row_string += sep + row_entries[col].rjust(column_widths[col])
            return row_string

        point_lines = [format_row(row, "|") for row in point_rows]
        lvec_lines = [format_row(row, ",") for row in lvec_rows_string]
        lvec_lines = [line.replace(',,', '; ') for line in lvec_lines]

        # Blank line after the leading nef-partition-indicator rows.
        num_nef_rows = len(self.nef_partition)
        point_lines = point_lines[:num_nef_rows] + [""] + point_lines[num_nef_rows:]

        return point_lines + ["-" * len(point_lines[-1])] + lvec_lines

    def _get_lvec_rows(self):
        """
        Returns the GLSM charge l-vectors (one row per Mori cone generator), with the
        nef-partition weight(s) appended, as a plain list of lists of ints - accessible
        as lvec[0], lvec[1], ... once loaded from the JSON output.
        """
        lvec_rows = [list(row[:-1]) for row in self.lvec.rows()]  # strip the total weight entry
        for row in lvec_rows:
            for partition in self.nef_partition:
                row.append(-sum([row[j] for j in partition]))  # Append the weight of the NEF partition
        return lvec_rows

    def _get_output(self):
        """
        Returns a dictionary containing all relevant topological data of the CY.
        """
        output = {
            "polytope_and_glsm_table": self._get_polytope_and_glsm_table(),
            "lvec": self._get_lvec_rows(),
            "kahler_cone_generators": [str(j) for j in self.J],
            "intersection_numbers_CY": dict(self.intersection_numbers_CY),
            "intersection_numbers_ambient_space": dict(self.intersection_numbers_ambient_space),
            "intersection_ring": self.intersection_ring,
            "intersection_ring_no_multiplicities": self.intersection_ring_no_multiplicities,
            "chern_polynomials": self.Chern_polynomials,
            "integrated_Chern_classes": self.integrated_Chern_classes,
        }

        debug_string = "\n\n--- Kähler cone generators (ambient space) ---\n" + str(output["kahler_cone_generators"]) + \
                        "\n\n--- Intersection numbers ambient space -------\n" + str(output["intersection_numbers_ambient_space"])

        info_string = "\n\n--- Polytope and GLSM table ------------------\n" + "\n".join(output["polytope_and_glsm_table"]) + \
                        "\n\n--- L-vectors (GLSM charges) -----------------\n" + str(output["lvec"]) + \
                        "\n\n--- Intersection numbers CY ------------------\n" + str(output["intersection_numbers_CY"]) + \
                        "\n\n--- Intersection ring ------------------------\n" + str(output["intersection_ring"]) + \
                        "\n\n--- Intersection ring no multiplicities ------\n" + str(output["intersection_ring_no_multiplicities"]) + \
                        "\n\n--- Chern polynomials ------------------------\n" + str(output["chern_polynomials"]) + \
                        "\n\n--- Integrated Chern classes -----------------\n" + str(output["integrated_Chern_classes"])

        logger.info(info_string)
        logger.debug(debug_string)

        return output

    @staticmethod
    def _to_json_safe(value):
        """
        Recursively converts a value into JSON-safe primitives: dicts/lists are
        recursed into (tuple dict keys become strings), numeric Sage objects
        (Integer, Rational, ...) become plain Python int/float, and anything else
        (e.g. symbolic expressions) is stringified.
        """
        if isinstance(value, dict):
            return {str(k): ToricPolytope._to_json_safe(v) for k, v in value.items()}
        if isinstance(value, (list, tuple)):
            return [ToricPolytope._to_json_safe(v) for v in value]
        if value is None or isinstance(value, (bool, str, int, float)):
            return value
        try:
            as_int = int(value)
            if as_int == value:
                return as_int
        except (TypeError, ValueError):
            pass
        try:
            return float(value)
        except (TypeError, ValueError):
            pass
        return str(value)

    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."""
        if not self.model_name:
            logger.info("No model_name given, skipping writing topological data to JSON.")
            return None

        repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
        out_dir = os.path.join(repo_root, "data/topdata")
        os.makedirs(out_dir, exist_ok=True)
        out_path = os.path.join(out_dir, "{}.json".format(self.model_name))

        with open(out_path, "w") as f:
            json.dump(self._to_json_safe(output), f, indent=2)

        logger.info("Wrote topological data to %s", out_path)
        return out_path

    def topdata(self):
        """
        Computes topological data for hypersurfaces and CICYs in toric ambient spaces
        and saves the result to /data/topdata/<model_name>.json if a model_name is provided.
        """
        # Each point gives a toric divisor, each "column" gives one linear relation
        self.no_divs = len(self.relevant_lattice_points) - self.dimension - 1

        # Cohomology ring and toric divisors D_0, ..., D_{N-1}.
        self.HH = self.ToricVariety.cohomology_ring()
        self.D = [self.HH(self.ToricVariety.divisor(i)) for i in range(len(self.relevant_lattice_points) - 1)]

        # Kähler cone generators J_a (dual to Mori cone generators l^a).
        self.J = self._compute_Kahler_generators()

        # Poincaré dual of the CY: product over nef partition parts.
        self.Ydualform = product([sum([self.D[p] for p in part]) for part in self.nef_partition])

        # Kähler parameter symbols t_0, ..., t_{nodivs-1}.
        self.t = var('t', n=self.no_divs, latex_name='t')

        # Intersection numbers κ_{i_1...i_n} = ∫_Y J_{i_1}∧...∧J_{i_n}.
        self.intersection_numbers_CY = {
            tuple(tl): self.ToricVariety.integrate(self.Ydualform * product([self.J[tl[j]] for j in range(self.cy_dimension)]))
            for tl in UnorderedTuples(range(self.no_divs), self.cy_dimension) # Unordered means only increasing tuples
        }
        self.intersection_numbers_ambient_space = {
            tuple(tl): self.ToricVariety.integrate(product([self.J[tl[j]] for j in range(self.dimension)]))
            for tl in UnorderedTuples(range(self.no_divs), self.dimension)
        }
        
        # Intersection ring polynomial and version without multiplicities.
        self.intersection_ring, self.intersection_ring_no_multiplicities = self._compute_intersection_ring()

        Chern_character = self._compute_Chern_character()
        
        self.Chern_polynomials, self.integrated_Chern_classes = self._compute_Chern_polynomials(Chern_character)

        output = self._get_output()
        self._write_output_json(output)

        self._check_for_fibrations()

class ToricPolytopeCICY(ToricPolytope):
    """
    Computes topological data for CICYs in toric ambient spaces.
    """

    def _CICY_to_points(self, CICY):
        """
        Converts a CICY configuration matrix to a list of points in the ambient space.
        Each row of the CICY corresponds to a projective space, and each column corresponds to a polynomial.
        The entries give the weights of the polynomials in the respective projective spaces.
        """
        ambients=(np.array(CICY).transpose()).sum(axis=0)-1
        if ambients in ZZ:
            ambients=[ambients]
        dimp=sum(ambients)
        polytope=identity_matrix(int(dimp))
        zeros=[0]*dimp
        polytope=polytope.insert_row(0,zeros)

        offset=0
        for i in range(len(ambients)):
            newrow=[0]*dimp
            entry=[]
            for j in range(ambients[i]):
                newrow[offset]=-1
                offset+=1
                entry.append(offset)
            polytope=polytope.insert_row(sum([ambients[n]+1 for n in [0..i]]),newrow)

        partition=[]
        eqs=matrix(CICY).transpose()
        counters=[0]*len(ambients)
        for eq in eqs:
            part=[]
            for j in range(len(eq)):
                part.append([sum([ambients[m]+1 for m in [0..(j-1)]])+counters[j]+a for a in [0..(eq[j]-1)]])
                counters[j]+=len(part[j])
            partition.append(flatten(part))
        
        return list(polytope), partition


    def __init__(self, CICY, no_triangulation=0, model_name=None, lvec=None):
        points, nef_partition = self._CICY_to_points(CICY)
        super().__init__(
            points,
            no_triangulation=no_triangulation,
            model_name=model_name,
            nef_partition=nef_partition,
            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
        )