# Source code for libpysal.cg.alpha_shapes

    """
    Computation of alpha shape algorithm in 2-D based on original implementation
    by Tim Kittel (@timkittel) available at:

    Author(s):
        Dani Arribas-Bel daniel.arribas.bel@gmail.com
        Levi John Wolf levi.john.wolf@gmail.com
    """

    import numpy as np
    import scipy.spatial as spat
    import shapely
    from scipy import sparse

    from ..common import HAS_JIT, jit, requires

    if not HAS_JIT:
        from warnings import warn

        NUMBA_WARN = (
            "Numba not imported, so alpha shape construction may be slower than expected."
        )

    EPS = np.finfo(float).eps

    __all__ = ["alpha_shape", "alpha_shape_auto"]

    @jit(nopython=True)
    def nb_dist(x, y):
        """numba implementation of distance between points `x` and `y`

        Parameters
        ----------
        x : ndarray
            Coordinates of point `x`
        y : ndarray
            Coordinates of point `y`

        Returns
        -------
        dist : float
            Distance between `x` and `y`

        Examples
        --------
        >>> x = np.array([0, 0])
        >>> y = np.array([1, 1])
        >>> dist = nb_dist(x, y)
        >>> dist
        1.4142135623730951
        """
        sum_ = 0
        for x_i, y_i in zip(x, y):  # noqa: B905
            sum_ += (x_i - y_i) ** 2
        dist = np.sqrt(sum_)
        return dist

    @jit(nopython=True)
    def r_circumcircle_triangle_single(a, b, c):
        """Computation of the circumcircle of a single triangle

        Parameters
        ----------
        a : ndarray
            (2,) Array with coordinates of vertex `a` of the triangle
        b : ndarray
            (2,) Array with coordinates of vertex `b` of the triangle
        c : ndarray
            (2,) Array with coordinates of vertex `c` of the triangle

        Returns
        -------
        r : float
            Circumcircle of the triangle

        Notes
        -----
        Source for equations:

        [Last accessed July 11th. 2018]

        Examples
        --------
        >>> a = np.array([0, 0])
        >>> b = np.array([0.5, 0])
        >>> c = np.array([0.25, 0.25])
        >>> r = r_circumcircle_triangle_single(a, b, c)
        >>> r
        0.2500000000000001
        """
        ab = nb_dist(a, b)
        bc = nb_dist(b, c)
        ca = nb_dist(c, a)

        num = ab * bc * ca
        den = np.sqrt((ab + bc + ca) * (bc + ca - ab) * (ca + ab - bc) * (ab + bc - ca))
        if den == 0:
            return np.array([ab, bc, ca]).max() / 2.0
        else:
            return num / den

    @jit(nopython=True)
    def r_circumcircle_triangle(a_s, b_s, c_s):
        """Computation of circumcircles for a series of triangles

        Parameters
        ----------
        a_s : ndarray
            (N, 2) array with coordinates of vertices `a` of the triangles
        b_s : ndarray
            (N, 2) array with coordinates of vertices `b` of the triangles
        c_s : ndarray
            (N, 2) array with coordinates of vertices `c` of the triangles

        Returns
        -------
        radii : ndarray
            (N,) array with circumcircles for every triangle

        Examples
        --------
        >>> a_s = np.array([[0, 0], [2, 1], [3, 2]])
        >>> b_s = np.array([[1, 0], [5, 1], [2, 4]])
        >>> c_s = np.array([[0, 7], [1, 3], [4, 2]])
        >>> rs = r_circumcircle_triangle(a_s, b_s, c_s)
        >>> rs
        array([3.53553391, 2.5       , 1.58113883])
        """
        len_a = len(a_s)
        r2 = np.zeros((len_a,))
        for i in range(len_a):
            r2[i] = r_circumcircle_triangle_single(a_s[i], b_s[i], c_s[i])
        return r2

    @jit(nopython=True)
    def get_faces(triangle):
        """Extract faces from a single triangle

        Parameters
        ----------
        triangles : ndarray
            (3,) array with the vertex indices for a triangle

        Returns
        -------
        faces : ndarray
            (3, 2) array with a row for each face containing the indices of the two
            points that make up the face

        Examples
        --------
        >>> triangle = np.array([3, 1, 4], dtype=np.int32)
        >>> faces = get_faces(triangle)
        >>> faces
        array([[3., 1.],
               [1., 4.],
               [4., 3.]])
        """
        faces = np.zeros((3, 2))
        for i, (i0, i1) in enumerate([(0, 1), (1, 2), (2, 0)]):
            faces[i] = triangle[i0], triangle[i1]
        return faces

    @jit(nopython=True)
    def build_faces(faces, triangles_is, num_triangles, num_faces_single):
        """Build facing triangles

        Parameters
        ----------
        faces : ndarray
            (num_triangles * num_faces_single, 2) array of zeroes in int form
        triangles_is : ndarray
            (D, 3) array, where D is the number of Delaunay triangles, with the
            vertex indices for each triangle
        num_triangles : int
            Number of triangles
        num_faces_single : int
            Number of faces a triangle has (i.e. 3)

        Returns
        -------
        faces : ndarray
            Two dimensional array with a row for every facing segment containing
            the indices of the coordinate points

        Examples
        --------
        >>> import scipy.spatial as spat
        >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
        >>> triangulation = spat.Delaunay(pts)
        >>> triangulation.simplices
        array([[3, 1, 4],
               [1, 2, 4],
               [2, 1, 0]], dtype=int32)
        >>> num_faces_single = 3
        >>> num_triangles = triangulation.simplices.shape[0]
        >>> num_faces = num_triangles * num_faces_single
        >>> faces = np.zeros((num_faces, 2), dtype=np.int_)
        >>> mask = np.ones((num_faces,), dtype=np.bool_)
        >>> faces = build_faces(
        ...     faces, triangulation.simplices, num_triangles, num_faces_single
        ... )
        >>> faces
        array([[3, 1],
               [1, 4],
               [4, 3],
               [1, 2],
               [2, 4],
               [4, 1],
               [2, 1],
               [1, 0],
               [0, 2]])
        """
        for i in range(num_triangles):
            from_i = num_faces_single * i
            to_i = num_faces_single * (i + 1)
            faces[from_i:to_i] = get_faces(triangles_is[i])
        return faces

    @jit(nopython=True)
    def nb_mask_faces(mask, faces):
        """Run over each row in `faces`, if the face in the following row is the
        same, then mark both as False on `mask`

        Parameters
        ----------
        mask : ndarray
            One-dimensional boolean array set to True with as many observations as
            rows in `faces`
        faces : ndarray
            Sorted sequence of faces for all triangles (ie. triangles split by each
            segment)

        Returns
        -------
        masked : ndarray
             Sequence of outward-facing faces

        Examples
        --------
        >>> import numpy as np
        >>> faces = np.array(
        ...     [
        ...         [0, 1], [0, 2], [1, 2], [1, 2], [1, 3], [1, 4], [1, 4], [2, 4], [3, 4]
        ...     ]
        ... )
        >>> mask = np.ones((faces.shape[0], ), dtype=np.bool_)
        >>> masked = nb_mask_faces(mask, faces)
        >>> masked
        array([[0, 1],
               [0, 2],
               [1, 3],
               [2, 4],
               [3, 4]])
        """
        for k in range(faces.shape[0] - 1):
            if mask[k] and np.all(faces[k] == faces[k + 1]):
                mask[k] = False
                mask[k + 1] = False
        return faces[mask]

    def get_single_faces(triangles_is):
        """Extract outward facing edges from collection of triangles

        Parameters
        ----------
        triangles_is : ndarray
            (D, 3) array, where D is the number of Delaunay triangles, with the
            vertex indices for each triangle

        Returns
        -------
        single_faces : ndarray

        Examples
        --------
        >>> import scipy.spatial as spat
        >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
        >>> alpha = 0.33
        >>> triangulation = spat.Delaunay(pts)
        >>> triangulation.simplices
        array([[3, 1, 4],
               [1, 2, 4],
               [2, 1, 0]], dtype=int32)
        >>> get_single_faces(triangulation.simplices)
        array([[0, 1],
               [0, 2],
               [1, 3],
               [2, 4],
               [3, 4]])
        """
        num_faces_single = 3
        num_triangles = triangles_is.shape[0]
        num_faces = num_triangles * num_faces_single
        faces = np.zeros((num_faces, 2), dtype=np.int_)
        mask = np.ones((num_faces,), dtype=np.bool_)

        faces = build_faces(faces, triangles_is, num_triangles, num_faces_single)

        orderlist = [f"x{i}" for i in range(faces.shape[1])]
        dtype_list = [(el, faces.dtype.str) for el in orderlist]
        # Arranging each face so smallest vertex is first
        faces.sort(axis=1)
        # Arranging faces in ascending way
        faces.view(dtype_list).sort(axis=0)
        # Masking
        single_faces = nb_mask_faces(mask, faces)
        return single_faces

    @requires("geopandas", "shapely")
    def _alpha_geoms(alpha, triangles, radii, xys):
        """Generate alpha-shape polygon(s) from `alpha` value, vertices of
        `triangles`, the `radii` for all points, and the points themselves

        Parameters
        ----------
        alpha : float
            Alpha value to delineate the alpha-shape
        triangles : ndarray
             (D, 3) array, where D is the number of Delaunay triangles, with the
             vertex indices for each triangle
        radii : ndarray
            (N,) array with circumcircles for every triangle
        xys : ndarray
            (N, 2) array with one point per row and coordinates structured as X and Y

        Returns
        -------
        geoms : GeoSeries
            Polygon(s) resulting from the alpha shape algorithm, in a GeoSeries.
            The output is a GeoSeries even if only a single polygon is returned. There is
            no CRS included in the returned GeoSeries.

        Examples
        --------
        >>> import scipy.spatial as spat
        >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
        >>> alpha = 0.33
        >>> triangulation = spat.Delaunay(pts)
        >>> triangles = pts[triangulation.simplices]
        >>> triangles
        array([[[6, 7],
                [3, 5],
                [9, 3]],
        <BLANKLINE>
               [[3, 5],
                [4, 1],
                [9, 3]],
        <BLANKLINE>
               [[4, 1],
                [3, 5],
                [0, 1]]])
        >>> a_pts = triangles[:, 0, :]
        >>> b_pts = triangles[:, 1, :]
        >>> c_pts = triangles[:, 2, :]
        >>> radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
        >>> geoms = alpha_geoms(alpha, triangulation.simplices, radii, pts)
        >>> geoms
        0    POLYGON ((0.00000 1.00000, 3.00000 5.00000, 4....
        dtype: geometry
        """
        from geopandas import GeoSeries
        from shapely.geometry import LineString
        from shapely.ops import polygonize

        triangles_reduced = triangles[radii < 1 / alpha]
        outer_triangulation = get_single_faces(triangles_reduced)
        face_pts = xys[outer_triangulation]
        geoms = GeoSeries(list(polygonize(list(map(LineString, face_pts)))))
        return geoms

    @requires("geopandas", "shapely")
    def alpha_shape(xys, alpha):
        """Alpha-shape delineation (Edelsbrunner, Kirkpatrick & Seidel, 1983)
        from a collection of points

        Parameters
        ----------
        xys : ndarray
            (N, 2) array with one point per row and coordinates structured as X and Y
        alpha : float
            Alpha value to delineate the alpha-shape

        Returns
        -------
        shapes : GeoSeries
             Polygon(s) resulting from the alpha shape algorithm. The GeoSeries
             object remains so even if only a single polygon is returned. There is
             no CRS included in the object. Note that the returned shape(s) may
             have holes, as per the definition of the shape in Edselbrunner et al. (1983)

        Examples
        --------
        >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
        >>> alpha = 0.1
        >>> poly = alpha_shape(pts, alpha)
        >>> poly
        0    POLYGON ((0.00000 1.00000, 3.00000 5.00000, 6....
        dtype: geometry
        >>> poly.centroid
        0    POINT (4.69048 3.45238)
        dtype: geometry

        References
        ----------
        Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
            a set of points in the plane. IEEE Transactions on information theory,
            29(4), 551-559.
        """
        if not HAS_JIT:
            warn(NUMBA_WARN, stacklevel=2)
        if xys.shape[0] < 4:
            from shapely import geometry as geom
            from shapely import ops

            return ops.unary_union([geom.Point(xy) for xy in xys]).convex_hull.buffer(0)
        triangulation = spat.Delaunay(xys)
        triangles = xys[triangulation.simplices]
        a_pts = triangles[:, 0, :]
        b_pts = triangles[:, 1, :]
        c_pts = triangles[:, 2, :]
        radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
        del triangles, a_pts, b_pts, c_pts
        geoms = _alpha_geoms(alpha, triangulation.simplices, radii, xys)
        geoms = _filter_holes(geoms, xys)
        return geoms

    def _valid_hull(geoms, points):
        """Sanity check within ``alpha_shape_auto()`` to verify the generated alpha
        shape actually contains the original set of points (xys).

        Parameters
        ----------
        geoms : GeoSeries
            See alpha_geoms()
        points : list
            xys parameter cast as shapely.geometry.Point objects

        Returns
        -------
        flag : bool
            Valid hull for alpha shape [True] or not [False]
        """
        # if there is not exactly one polygon
        if geoms.shape[0] != 1:
            return False
        # if any (xys) points do not intersect the polygon
        return shapely.intersects(geoms[0], points).all()

    @requires("geopandas", "shapely")
    def alpha_shape_auto(
        xys, step=1, verbose=False, return_radius=False, return_circles=False
    ):
        """Computation of alpha-shape delineation with automated selection of alpha.

        This method uses the algorithm proposed by  Edelsbrunner, Kirkpatrick &
        Seidel (1983) to return the tightest polygon that contains all points in
        `xys`. The algorithm ranks every point based on its radius and iterates
        over each point, checking whether the maximum alpha that would keep the
        point and all the other ones in the set with smaller radii results in a
        single polygon. If that is the case, it moves to the next point;
        otherwise, it retains the previous alpha value and returns the polygon
        as `shapely` geometry. Note that this geometry may have holes.

        Parameters
        ----------

        xys : ndarray
            Nx2 array with one point per row and coordinates structured as X and Y
        step : int
            [Optional. Default=1] Number of points in `xys` to jump ahead after
            checking whether the largest possible alpha that includes the point and
            all the other ones with smaller radii
        verbose : Boolean
            [Optional. Default=False] If True, it prints alpha values
            being tried at every step.

        Returns
        -------
        poly : shapely.Polygon
             Tightest alpha-shape polygon containing all points in `xys`

        Examples
        --------
        >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
        >>> poly = alpha_shape_auto(pts)
        >>> poly.bounds
        (0.0, 1.0, 9.0, 7.0)
        >>> poly.centroid.x, poly.centroid.y
        (4.690476190476191, 3.4523809523809526)

        References
        ----------
        Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
            a set of points in the plane. IEEE Transactions on information theory,
            29(4), 551-559.

        """
        if not HAS_JIT:
            warn(NUMBA_WARN, stacklevel=2)
        from shapely import geometry as geom

        if return_circles:
            return_radius = True
        if xys.shape[0] < 4:
            if xys.shape[0] == 3:
                multipoint = geom.MultiPoint(xys)
                alpha_shape = multipoint.convex_hull
            else:
                alpha_shape = geom.Polygon([])
            if xys.shape[0] == 1:
                if return_radius:
                    if return_circles:
                        out = [alpha_shape, 0, alpha_shape]
                    return alpha_shape, 0
                return alpha_shape
            elif xys.shape[0] == 2:
                if return_radius:
                    r = spat.distance.euclidean(xys[0], xys[1]) / 2
                    if return_circles:
                        circle = _construct_centers(xys[0], xys[1], r)
                        return [alpha_shape, r, circle]
                    return [alpha_shape, r]
                return alpha_shape
            elif return_radius:  # this handles xys.shape[0] == 3
                radius = r_circumcircle_triangle_single(xys[0], xys[1], xys[2])
                if return_circles:
                    circles = construct_bounding_circles(alpha_shape, radius)
                    return [alpha_shape, radius, circles]
                return [alpha_shape, radius]
            return alpha_shape
        triangulation = spat.Delaunay(xys)
        triangles = xys[triangulation.simplices]
        a_pts = triangles[:, 0, :]
        b_pts = triangles[:, 1, :]
        c_pts = triangles[:, 2, :]
        radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
        radii[np.isnan(radii)] = 0  # "Line" triangles to be kept for sure
        del triangles, a_pts, b_pts, c_pts
        radii_sorted_i = radii.argsort()
        triangles = triangulation.simplices[radii_sorted_i][::-1]
        radii = radii[radii_sorted_i][::-1]
        geoms_prev = _alpha_geoms((1 / radii.max()) - EPS, triangles, radii, xys)
        points = shapely.points(xys)
        if verbose:
            print(f"Step set to {step}")
        for i in range(0, len(radii), step):
            radi = radii[i]
            alpha = (1 / radi) - EPS
            if verbose:
                print(f"{(i + 1) / radii.shape[0]:.2f}% | Trying a = {alpha:f}")
            geoms = _alpha_geoms(alpha, triangles, radii, xys)
            if _valid_hull(geoms, points):
                geoms_prev = geoms
                radi_prev = radi
            else:
                break
        if verbose:
            print(geoms_prev.shape)
        if return_radius:
            out = [geoms_prev[0], radi_prev]
            if return_circles:
                out.append(construct_bounding_circles(out[0], radi_prev))
            return out
        # Return a shapely polygon
        return geoms_prev[0]

    def construct_bounding_circles(alpha_shape, radius):
        """Construct the bounding circles for an alpha shape, given the radius
        computed from the `alpha_shape_auto` method.

        Parameters
        ----------
        alpha_shape : shapely.Polygon
            An alpha-hull with the input radius.
        radius : float
            The radius of the input alpha_shape.

        Returns
        -------
        center : numpy.ndarray of shape (n,2)
            The centers of the circles defining the alpha_shape.
        """
        coordinates = list(alpha_shape.boundary.coords)
        n_coordinates = len(coordinates)
        centers = []
        for i in range(n_coordinates - 1):
            a, b = coordinates[i], coordinates[i + 1]
            centers.append(_construct_centers(a, b, radius))
        return centers

    @jit(nopython=True)
    def _construct_centers(a, b, radius):
        midpoint_x = (a[0] + b[0]) * 0.5
        midpoint_y = (a[1] + b[1]) * 0.5
        d = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5
        if b[0] - a[0] == 0:
            m = np.inf
            axis_rotation = np.pi / 2
        else:
            m = (b[1] - a[1]) / (b[0] - a[0])
            axis_rotation = np.arctan(m)
        # altitude is perpendicular bisector of AB
        interior_angle = np.arccos(0.5 * d / radius)
        chord = np.sin(interior_angle) * radius

        dx = chord * np.sin(axis_rotation)
        dy = chord * np.cos(axis_rotation)

        up_x = midpoint_x - dx
        up_y = midpoint_y + dy
        down_x = midpoint_x + dx
        down_y = midpoint_y - dy

        # sign gives us direction of point, since
        # shapely shapes are clockwise-defined
        sign = np.sign((b[0] - a[0]) * (up_y - a[1]) - (b[1] - a[1]) * (up_x - a[0]))
        if sign == 1:
            return up_x, up_y
        else:
            return down_x, down_y

    def _filter_holes(geoms, points):  # noqa: ARG001
        """
        Filter hole polygons using a computational geometry solution
        """

        if (geoms.interiors.apply(len) > 0).any():
            from shapely.geometry import Polygon

            # Extract the "shell", or outer ring of the polygon.
            shells = geoms.exterior.apply(Polygon)
            # Compute which original geometries are within each shell, self-inclusive
            inside, outside = shells.sindex.query(geoms, predicate="within")

            # Now, create the sparse matrix relating the inner geom (rows)
            # to the outer shell (cols) and take the sum.
            # A z-order of 1 means the polygon is only inside if its own exterior.
            #   This means it's not a hole.
            # A z-order of 2 means the polygon is inside of exactly one other exterior.
            #   Because the hull generation method is restricted to be planar, this
            #   means the polygon is a hole.
            # In general, an even z-order means that the polygon is always exactly
            #   matched to one exterior, plus some number of intermediate
            #   exterior-hole pairs. Therefore, the polygon is a hole.
            # In general, an odd z-order means that there is an uneven number of exteriors.
            #   This means the polygon is not a hole.
            zorder = sparse.csc_matrix((np.ones_like(inside), (inside, outside))).sum(
                axis=1
            )
            zorder = np.asarray(zorder).flatten()
            # Keep only the odd z-orders
            to_include = (zorder % 2).astype(bool)
            geoms = geoms[to_include]
        return geoms

    if __name__ == "__main__":
        import time

        import geopandas as gpd
        import matplotlib.pyplot as plt

        plt.close("all")
        xys = np.random.random((1000, 2))
        t0 = time.time()
        geoms = alpha_shape_auto(xys, 1)
        t1 = time.time()
        print("%.2f Seconds to run algorithm" % (t1 - t0))
        f, ax = plt.subplots(1)
        gpd.GeoDataFrame({"geometry": [geoms]}).plot(ax=ax, color="orange", alpha=0.5)
        ax.scatter(xys[:, 0], xys[:, 1], s=0.1)
        plt.show()